1. Calculation of the energy balance of a calculation#
The energy balance of a calculation can be very useful information. On the one hand, it makes it possible to take a critical look at the quality of the results, in the same way as an equilibrium residue in Newton’s algorithm. On the other hand, it provides information on the physical phenomena involved and provides answers to the following questions: what is the energy transmitted to my system? What is she becoming?
To do this, the energy balance for a calculation with Code_Aster is established below on the complete model. We distinguish the nature of each energy term and strive to define where it comes from.
1.1. Work of a force#
Let us consider the evolution at different times \({t}_{i}\) of a deformable solid subjected to a variable external force \(\mathrm{F}(t)\). Let \({\mathrm{U}}_{i}\) be the solution displacement field for each moment \({t}_{i}\) and \({\mathrm{F}}_{i}\mathrm{=}\mathrm{F}({t}_{i})\). The curve of the external force as a function of the displacement is given on the. From an energy point of view, the increment of external work provided to the system between two consecutive moments \(({t}_{0},{t}_{1})\) represents the area under the curve, which is calculated by the trapezius method:
where:
\({\stackrel{ˉ}{\mathrm{F}}}_{1}\) represents the average force on the time increment (the index is chosen to be that of the end time),
\({\Delta \mathrm{U}}_{1}\) represents the displacement increment between \({t}_{0}\) and \({t}_{1}\).
Figure 1.1-1: Force-displacement curve.
The work increment of any force can therefore be considered as the dot product of two vectors: the mean force over the time step and the displacement increment.
This method for calculating the work increment of a force is practical in the context of numerical simulation, since the solution is known only at different times, corresponding to the time discretization chosen.
1.2. Energy balance#
Let us consider the dynamic equilibrium equation of a system in its spatially discretized form:
where:
\(\mathrm{U}\), \(\dot{\mathrm{U}}\) and \(\ddot{\mathrm{U}}\) are the vectors of movements, speeds and accelerations at the nodes,
\(\mathrm{M}\) and \(\mathrm{C}\) are the mass and damping matrices,
\(\mathrm{R}\) represents internal forces (in linear elasticity, this is the term \(\mathrm{K}\mathrm{U}\)), from which what is considered to be damping forces has been removed,
\(\mathrm{L}\) is the vector of external efforts.
This writing, inspired by [R5.05.05], is here simplified compared to the original script which involves Lagrange multipliers for Dirichlet boundary conditions and unilateral conditions. Note that damping matrix \(\mathrm{C}\left(t\right)\) may depend on time if its construction is based on the tangent stiffness matrix and not elastic (option AMOR_RAYL_RIGI in operator DYNA_NON_LINE).
Each term in equation () is homogeneous to a force. The energy evolution of the system between two consecutive calculation moments \({t}_{0}\) and \({t}_{1}\) is directly reflected in the work of each of these forces. Using the mean magnitude and magnitude increment notations used in equation (), we obtain:
the kinetic energy increment:
An expression identical to equation () is found by considering \(\mathrm{F}=\mathrm{M}\ddot{\mathrm{U}}\) as force, with the \(\mathrm{M}\) operator being assumed to be constant.
the increment of energy dissipated by damping:
: label: eq-4
mathrm {Delta} {W} _ {mathit {amor}}} = {overline {mathrm {C}left (mathrm {t}right)dot {mathrm {U}}}} _ {1}} ^ {mathrm {T}}mathrm {.} mathrm {Delta} {mathrm {U}} _ {1}
Since the damping matrix can be time-dependent, the quantity \(\mathrm{C}\left(\mathrm{t}\right)\dot{\mathrm{U}}\) must be evaluated at \({t}_{0}\) and \({t}_{1}\) must be evaluated and then averaged. This step is neglected, and the following quantity is in fact calculated:
: label: eq-5
mathrm {Delta} {W} _ {mathit {amor}}} =overline {{dot {mathrm {U}}}} _ {1} ^ {mathrm {T}}}}mathrm {T}}}}mathrm {C}}}left ({t} _ {1}right)mathrm {delta} {mathrm {U}}}}}mathrm {C}}left ({t} _ {1}right)mathrm {delta} {mathrm {U}}}}}mathrm {C}}left ({t} _ {1}right)mathrm {delta} {mathrm {U}}}}}
The energy dissipated by damping is therefore calculated in an approximate manner.
The total deformation energy increment:
The work increment of external efforts:
In the presence of rubbing contact, an additional effort term appears in the equilibrium equation, corresponding to the contact force. From an energy point of view, its contribution is easily calculated in a manner similar to (), but it cannot be taken into account in any of the energies defined above. We therefore create a new energy that represents the energy stored and/or dissipated by the bonding forces, which we call \(\Delta {W}_{\mathit{liai}}\). The mechanical energy balance of the system is then written as:
Assuming that the resolution is strictly accurate, the term on the left is equal to the sum of the terms on the right. In reality, the numerical resolution scheme induces differences for various reasons:
Newton’s algorithm is an iterative resolution algorithm, whose convergence is detected when the residual relative balance of forces is less than a certain criterion. By definition, the solution found is therefore correct within a tolerance, which inevitably creates a balance gap in the energy balance. It is then possible to improve the energy balance by tightening the convergence criterion, but the precision achieved with the criterion by default is generally satisfactory.
Among the numerical schemes of temporal integration, only some are not dispative. Changing the default values of these diagrams, such as the Newmark schema, causes numerical dissipation. The user can refer to the [R5.05.05] documentation for more information on the dissipation of numerical temporal integration schemes.
The energy balance is of first importance here: it makes it possible to quantify the energy dissipated by the numerical diagram. In shock calculations, for example, it is common to use dissipative numerical schemes to filter high frequencies. The energy dissipated by the numerical diagram is defined as follows:
1.3. Construction of the different energy terms#
Energy is considered to be the work of a force. Numerically, calculating an energy term consists first of all in constructing a global vector corresponding to the sum of all the forces whose contribution to the energy balance is stored in the same term. There are many force vectors and in this section we define which ones feed each of the energy terms.
1.3.1. Work of external forces#
External work is defined as work:
of all forces associated with commands AFFE_CHAR_MECA and AFFE_CHAR_MECA_F, with the exception of FORCE_SOLcomptabilisée in the work of the liaison forces,
the forces resulting from the calculation by substructuration (use of macro-elements),
of all forces associated with commands AFFE_CHAR_CINEetAFFE_CHAR_CINE_F. In the strict sense of the word, the AFFE_CHAR_CINE command makes it possible to define conditions at the limits of imposed displacement, whose work should be counted as binding energy. For practical reasons, its contribution is added to external work. In general, it is possible to recover the force corresponding to this imposed displacement. However, some cases may cause problems, such as the ssnop155a test case (see § 2.1.2.1).
1.3.2. Liaison forces work#
The work of liaison forces is defined as work:
contact and friction forces from command DEFI_CONTACT,
the force coming from the absorbing border elements,
of the force from option FORCE_SOLdans the AFFE_CHAR_MECA command.
1.3.3. Work of damping forces#
Work of amortization forces is defined as work:
the damping force from the damping matrix (term \(\mathrm{C}\dot{\mathrm{U}}=(\alpha \mathrm{K}+\beta \mathrm{M})\dot{\mathrm{U}}\)),
the modal damping force (keyword AMOR_MODAL in calculation operators).
1.3.4. Total deformation energy#
The total deformation energy includes all the contributions that require the integration of a law of behavior. This can lead to arbitrary choices. For example, the contribution of joint elements or cohesive elements is included in the total deformation energy, although it would be more natural to include them in the calculation of the energy linked to the bonds. Likewise, when laws of behavior are used to represent contact with discrete elements, the contribution is integrated into the total deformation energy whereas it would be more appropriate in the energy linked to the bonds. The rule is that, when it comes to behavior carried by finite elements defined in the model, then the contribution is added to \({E}_{\mathit{tot}}\).
1.4. Taking into account an initial state#
Let us consider the example of a mass-spring system subjected to an initial shock. In Code_Aster, the user can simply give an initial acceleration. This corresponds to the following equation:
\(\begin{array}{c}m\ddot{u}\left(t\right)+ku\left(t\right)=F\left(t\right)\\ u=0\\ F=0\\ \ddot{u}\ne 0\end{array}\) with as initial conditions: \(\begin{array}{c}u(0)=0\\ \dot{u}(0)=0\\ F(0)=0\\ \ddot{u}(0)\ne 0\end{array}\) |
It is noted that there is an incompatibility: the user provides an unbalanced initial state. In all seriousness, since the equations of dynamics are of the second order in terms of time, there can only be a displacement or a speed imposed at the initial moment. If we want to impose an initial acceleration, we should in principle give the force that balances the equation of the system. The numerical diagram does not need it because it does not affect the resolution. However, this affects the calculation of the energy balance because without this initial force, the balance is incorrect. In fact, the energy balance over the first step of time is then written as:
The contribution of external force that corresponds to this initial acceleration is missing. Taking it into account, we then get a term \({W}_{\mathit{ext}}\) that is not zero.
In the general case, the energy balance can only be correct if the calculation starts from a balanced state. In the absence of an initial state, this condition is implicitly met. Particular attention must therefore be paid to the various possibilities of introducing a non-virgin initial state.
1.4.1. Case of returns based on a previously calculated statement#
In this part, we are interested in repeats based on a previously calculated state, namely the use of keywords in ETAT_INIT:
EVOL_NOLI for STAT_NON_LINEet DYNA_NON_LINE,
RESULTATpour DYNA_VIBRA.
In order to calculate each of the energy terms in equation (), it is therefore necessary to have at the recovery step:
nodal vectors of movement, speed and accelerations,
mass and damping matrices,
internal forces resulting from the stress state,
force vectors \({\mathrm{F}}_{\mathit{ext}}\), \({\mathrm{F}}_{\mathit{amor}}\) and \({F}_{\mathit{liai}}\).
The nodal vectors of movement, speed, and acceleration, as well as the mass matrix and internal forces, are still available. The damping matrix is in fact not necessary since equation () is used to evaluate the energy dissipated by damping, where this matrix is evaluated at time \({t}_{1}\). Only the last three force vectors require special attention.
Vectors \({\mathrm{F}}_{\mathit{amor}}\) and \({\mathrm{F}}_{\mathit{liai}}\) must have been archived as it is impossible to recreate them. Thus, when the calculation of the energy balance is activated, the archiving of two additional nodal fields called FORC_AMOR and FORC_LIAI is automatically activated. These fields are intended to be re-read in case of resumption in order to correctly initialize the calculation of the energy balance.
Finally, concerning the vector \({\mathrm{F}}_{\mathit{ext}}\) of external forces, for STAT_NON_LINE and DYNA_NON_LINE, it is recalculated at the recovery time \({t}_{0}\) in accordance with the system equilibrium equation ():
In the case of a calculation with DYNA_VIBRA on a physical basis, the external force FORC_EXTE is saved directly.
1.4.2. Case of a calculation with an explicitly given initial state#
The user may be required to provide a non-blank initial state in the form of a combination of the following fields:
a field for movement (keyword DEPL)
a speed field (keyword VITE)
an acceleration field (keyword ACCE)
a constraint field (keyword SIGM)
a field of internal variables (keyword VARI)
The data for the displacement fields, speed and acceleration is automatically managed as in the case of an initial state resulting from a previous calculation, using equation (). Note that the time diagrams of the dynamics operators DYNA_NON_LINEet DYNA_VIBRA are required to calculate an initial acceleration if the user does not provide one. This case is handled in the same way. As an example, in the case of the mass-spring system subject to initial acceleration, we construct the following initial vector \({F}_{\mathit{ext}}\):
The data of a constraint field or internal variables is also automatically managed, since these quantities are used in the calculation of the term \(\mathrm{R}(\mathrm{U}({t}_{0}),\dot{\mathrm{U}}({t}_{0}))\) of equation ().
1.4.3. Initial energy state#
It is important to note that the energy balance is provided in incremental form. Thus, it is only possible to calculate a variation of the various energy terms between two calculation moments. Moreover, it is this information that is useful in order to know the different energy transfers. When the initial state of a calculation is provided from a result (see 1.4.1), the initial energy state is recovered. On the other hand, when the initial state is not given by a « result » data structure, the initial energy state cannot be recovered or calculated. As a result, the various terms of the energy balance are initialized to 0.