3. Energy balance equation setup for virtual test models

In a virtual test, the energy balance equation is the accounting ledger that tells you whether the solver, the contact definitions, and the boundary conditions are all behaving consistently. If the ledger does not close, the model is not trustworthy — regardless of how plausible the animation looks. This section develops the energy balance for a general explicit or implicit transient model, defines each term with a consistent sign convention, and shows how to assemble and audit the equation in practice.

Sign convention and the master equation

Before writing any term, fix the reference frame and the sign convention, because every downstream check depends on them. Adopt the following convention throughout this course:

  • Work done on the system is positive; work done by the system is negative.
  • Energy stored in the system is positive; energy leaving the system is negative.
  • Gravity acts along the negative global Z axis; potential energy increases with height.
  • Contact forces are internal to the assembly but external to each individual body, so contact work appears with opposite signs on the two bodies.

With that convention, the discrete energy balance at any time t is:

E_kin(t) + E_pot(t) + E_int(t) + E_con(t)
  = W_ext(t) + E_kin(0) + E_pot(0) + E_int(0) + E_con(0) - E_diss(t)

Here E_kin is kinetic energy, E_pot gravitational potential energy, E_int internal (strain) energy, and E_con the contact energy stored in penalty springs or contact compliance. W_ext is work by external loads and prescribed motion, and E_diss is energy removed by damping, friction, plasticity, and hourglass control. The residual of this equation, normalized by the initial or peak energy, is the single most useful health metric for a virtual test.

Kinetic energy term

For a lumped-mass explicit model, kinetic energy is summed over nodes:

E_kin = 0.5 * Σ ( m_i * |v_i|^2 )

Two practical points matter. First, use the current nodal velocities at the same output frame as the other terms; mixing frames introduces a spurious residual that looks like a solver problem but is a post-processing error. Second, for models with rigid bodies, include rotational kinetic energy 0.5 * ω^T I ω; omitting it is a common cause of a balance that drifts upward during spin-up. In a typical crash or drop model, E_kin starts at its maximum and decays as deformation and contact absorb it.

Potential energy term

Gravitational potential energy is referenced to a datum plane z_ref:

E_pot = Σ ( m_i * g * (z_i - z_ref) )

Choose z_ref once and record it in the model metadata. If the datum is changed between runs, the absolute values shift but the change in potential energy must still match the work done by gravity. For a drop test, the initial potential energy above the impact plane is the energy input you expect to see converted into kinetic energy just before contact — a useful independent check on the mass and gravity settings.

Internal energy term

Internal energy is the recoverable strain energy plus the irrecoverable dissipated energy in the material:

E_int = Σ_elements [ ∫ σ : dε  dV ]  =  E_elastic + E_plastic + E_damage

In practice you extract E_int from the solver's element energy output (often labeled IE or ALLIE). For an elastic material it oscillates and returns to near zero after unloading; for an elastic-plastic material it grows monotonically and never returns to zero. That distinction is itself a validation signal: if a nominally elastic model shows a permanently rising E_int, check for unintended plasticity, excessive hourglassing, or a misassigned material.

Contact energy term

Penalty-based contact stores energy in the artificial spring between penetrating surfaces:

E_con = Σ_contacts 0.5 * k_p * δ^2

where k_p is the penalty stiffness and δ the penetration depth. E_con should remain a small fraction of the total energy — a common acceptance threshold is below 5% of peak internal energy, and below 1% for a well-tuned model. A large or spiking E_con indicates penalty stiffness that is too low (deep penetration), too high (instability and tiny time steps), or contact definitions that miss initial penetrations. Contact friction dissipates energy rather than storing it, so friction work belongs in E_diss, not E_con.

Report contact energy as a percentage of the peak internal energy at every output frame, not just at the end. A contact energy spike that appears and then vanishes often marks a transient penetration event that a final-frame check would miss entirely.

Assembling and auditing the balance in practice

The workflow below assumes a transient explicit model with a rigid impactor and a deformable target, exported at a fixed output interval. The steps are generic and apply to most explicit solvers.

  1. Request all energy components at every output frame. Enable kinetic, internal, contact, external work, and damping/friction energy in the output request. Do not reconstruct them from stresses afterward; use the solver's consistent set.
  2. Fix the datum and gravity vector in the model header and confirm the potential energy reference matches the post-processor's assumption.
  3. Compute the residual per frame using the master equation, then normalize by the maximum total energy observed over the run.
  4. Plot the normalized residual against time. A healthy explicit run keeps it below roughly 1–2%; a steadily climbing residual points to missing dissipation or a post-processing mismatch.
  5. Cross-check the drop phase. Before first contact, the loss in potential energy should equal the gain in kinetic energy to within numerical damping.

As a concrete example, consider a 2.0 kg impactor dropped from 0.5 m onto a deformable plate, with gravity 9.81 m/s². The initial potential energy relative to the impact plane is 2.0 × 9.81 × 0.5 ≈ 9.81 J. Just before impact, kinetic energy should read approximately 9.81 J minus any damping losses. If the post-processor reports 8.4 J, the missing 1.4 J is either real damping or a mass/gravity error — and the residual plot will tell you which, because a mass error shifts the whole curve while damping produces a gradual divergence. These numbers are illustrative and should be recomputed for your own model.

The figure below shows how the four energy terms exchange during a typical impact sequence and where the residual is evaluated.

Energy Exchange During an Impact Virtual Test Time → Energy (J) → Free fall Impact / deformation Rebound / settle E_pot E_kin E_int E_con Residual check R = (E_kin+E_pot+E_int+E_con − W_ext − E_0 + E_diss) / E_max Potential Kinetic Internal Contact

Common failure symptoms and how to read them

The residual curve is a diagnostic, not just a pass/fail number. A few patterns recur often enough to memorize:

Residual patternLikely causeFirst check
Steady upward driftMissing dissipation term (friction or damping not requested)Output request for friction/damping energy
Step jump at first contactInitial penetration or penalty stiffness mismatchContact gap and k_p settings
Oscillating residualOutput frames not synchronized across energy componentsOutput interval consistency
Large constant offsetWrong potential energy datum or missing rotational kinetic energyDatum plane and rigid-body inertia

When the residual stays within tolerance and the drop-phase cross-check closes, the energy ledger is sound and the model is ready for the contact artifact review that follows. If it does not close, resolve the balance before interpreting any contact result — an unbalanced ledger will contaminate every downstream conclusion about penetration, force, and deformation.