Solid Mechanics Solver

List of Symbols

$i,j,k &\equiv \text {indices over spatial dimensions} \notag \\ a,b,c &\equiv \text {indices over nodes} \notag \\ l &\equiv \text {indices over volumetric elements} \notag \\ q,r,s &\equiv \text {indices over faces} \notag \\ n &\equiv \text {indices over time} \notag \\ kiter &\equiv \text {iteration count for non-linear solution scheme} \notag \\ \Omega &\equiv \text {Volume of continuum body} \notag \\ \Omega_{crack} &\equiv \text {Volume of open crack} \notag \\ \Gamma &\equiv \text {External surface of } \Omega \notag \\ \Gamma_t &\equiv \text {External surface where tractions are applied} \notag \\ \Gamma_u &\equiv \text {External surface where kinematics are specified} \notag \\ \Gamma_{crack} &\equiv \text {entire surface of crack} \notag \\ \Gamma_{cohesive} &\equiv \text {surface of crack subject to cohesive tractions} \notag \\ \eta_0 &\equiv \text {set of all nodes} \notag \\ \eta_f & \equiv \text {set of all nodes on flow mesh} \notag \\ m &\equiv \text{mass} \notag \\ \kappa_k &\equiv \text{all elements connected to element k} \notag \\ \phi & \equiv \text { porosity} \notag \\ p_f & \equiv \text { fluid pressure} \notag \\ \mathbf{u} & \equiv \text { displacement} \notag \\ \mathbf{q} & \equiv \text { volumetric flow rate} \notag \\ \mathbf{T} & \equiv \text { Cauchy stress} \notag \\ \rho & \equiv \text { density in the current configuration} \notag \\ \mathbf{x}& \equiv \text { current position} \notag \\ \mathbf{w}&\equiv \text { aperture, or gap vector} \notag$

Introduction

The SolidMechanics_LagrangianFEM solver applies a Continuous Galerkin finite element method to solve the linear momentum balance equation. The primary variable is the displacement field which is discretized at the nodes.

Theory

Governing Equations

The SolidMechanics_LagrangianFEM solves the equations of motion as given by

$T_{ij,j} + \rho(b_{i}-\ddot{x}_{i}) = 0,$

which is a 3-dimensional expression for the well known expression of Newtons Second Law ( $F = m a$ ). These equations of motion are discretized using the Finite Element Method, which leads to a discrete set of residual equations:

$(R_{solid})_{ai}=\int\limits_{\Gamma_t} \Phi_a t_i dA - \int\limits_\Omega \Phi_{a,j} T_{ij} dV +\int\limits_\Omega \Phi_a \rho(b_{i}-\Phi_b\ddot{x}_{ib}) dV = 0$

Quasi-Static Time Integration

The Quasi-Static time integration option solves the equation of motion after removing the inertial term, which is expressed by

$T_{ij,j} + \rho b_{i} = 0,$

which is essentially a way to express the equation for static equilibrium ( $\Sigma F=0$ ). Thus, selection of the Quasi-Static option will yield a solution where the sum of all forces at a given node is equal to zero. The resulting finite element discretized set of residual equations are expressed as

$(R_{solid})_{ai}=\int\limits_{\Gamma_t} \Phi_a t_i dA - \int\limits_\Omega \Phi_{a,j} T_{ij} dV + \int\limits_\Omega \Phi_a \rho b_{i} dV = 0,$

Taking the derivative of these residual equations wrt. the primary variable (displacement) yields

$\pderiv{(R_{solid}^e)_{ai}}{u_{bj}} &= - \int\limits_{\Omega^e} \Phi_{a,k} \frac{\partial T_{ik}}{\partial u_{bj}} dV,$

And finally, the expression for the residual equation and derivative are used to express a non-linear system of equations

$\left. \left(\pderiv{(R_{solid}^e)_{ai}}{u_{bj}} \right)\right|^{n+1}_{kiter} \left( \left. \left({u}_{bj} \right) \right|^{n+1}_{{kiter}+1} - \left. \left({u}_{bj} \right) \right|^{n+1}_{kiter} \right) = - (R_{solid})_{ai}|^{n+1}_{kiter} ,$

which are solved via the solver package.

Implicit Dynamics Time Integration (Newmark Method)

For implicit dynamic time integration, we use an implementation of the classical Newmark method. This update method can be posed in terms of a simple SDOF spring/dashpot/mass model. In the following, $M$ represents the mass, $C$ represent the damping of the dashpot, $K$ represents the spring stiffness, and $F$ represents some external load.

$M a^{n+1} + C v^{n+1} + K u^{n+1} &= F_{n+1}, \\$

and a series of update equations for the velocity and displacement at a point:

$u^{n+1} &= u^n + v^{n+1/2} \Delta t, \\ u^{n+1} &= u^n + \left( v^{n} + \inv{2} \left[ (1-2\beta) a^n + 2\beta a^{n+1} \right] \Delta t \right) \Delta t, \\ v^{n+1} &= v^n + \left[(1-\gamma) a^n + \gamma a^{n+1} \right] \Delta t.$

As intermediate quantities we can form an estimate (predictor) for the end of step displacement and midstep velocity by assuming zero end-of-step acceleration.

$\tilde{u}^{n+1} &= u^n + \left( v^{n} + \inv{2} (1-2\beta) a^n \Delta t \right) \Delta t = u^n + \hat{\tilde{u}}\\ \tilde{v}^{n+1} &= v^n + (1-\gamma) a^n \Delta t = v^n + \hat{\tilde{v}}$

This gives the end of step displacement and velocity in terms of the predictor with a correction for the end step acceleration.

$u^{n+1} &= \tilde{u}^{n+1} + \beta a^{n+1} \Delta t^2 \\ v^{n+1} &= \tilde{v}^{n+1} + \gamma a^{n+1} \Delta t$

The acceleration and velocity may now be expressed in terms of displacement, and ultimately in terms of the incremental displacement.

$a^{n+1} &= \frac{1}{\beta \Delta t^2} \left(u^{n+1} - \tilde{u}^{n+1} \right) = \frac{1}{\beta \Delta t^2} \left( \hat{u} - \hat{\tilde{u}} \right) \\ v^{n+1} &= \tilde{v}^{n+1} + \frac{\gamma}{\beta \Delta t} \left(u^{n+1} - \tilde{u}^{n+1} \right) = \tilde{v}^{n+1} + \frac{\gamma}{\beta \Delta t} \left(\hat{u} - \hat{\tilde{u}} \right)$

plugging these into equation of motion for the SDOF system gives:

$M \left(\frac{1}{\beta \Delta t^2} \left(\hat{u} - \hat{\tilde{u}} \right)\right) + C \left( \tilde{v}^{n+1} + \frac{\gamma}{\beta \Delta t} \left(\hat{u} - \hat{\tilde{u}} \right) \right) + K u^{n+1} &= F_{n+1} \\$

Finally, we assume Rayliegh damping for the dashpot.

$C = a_{mass} M + a_{stiff} K$

Of course we know that we intend to model a system of equations with many DOF. Thus the representation for the mass, spring and dashpot can be replaced by our finite element discretized equation of motion. We may express the system in context of a nonlinear residual problem

$(R_{solid}^e)_{ai} &= \int\limits_{\Gamma_t^e} \Phi_a t_i dA \\ &- \int\limits_{\Omega^e} \Phi_{a,j} \left(T_{ij}^{n+1}+ a_{stiff} \left(\pderiv{T_{ij}^{n+1}}{\hat{u}_{bk}} \right)_{elastic} \left( \tilde{v}_{bk}^{n+1} + \frac{\gamma}{\beta \Delta t} \left(\hat{u}_{bk} - \hat{\tilde{u}}_{bk} \right) \right) \right) dV \notag \\ &+\int\limits_{\Omega^e} \Phi_a \rho \left(b_{i}- \Phi_b \left( a_{mass} \left( \tilde{v}_{bi}^{n+1} + \frac{\gamma}{\beta \Delta t} \left(\hat{u}_{bi} - \hat{\tilde{u}}_{bi} \right) \right) + \frac{1}{\beta \Delta t^2} \left( \hat{u}_{bi} - \hat{\tilde{u}}_{bi} \right) \right) \right) dV ,\notag \\ \pderiv{(R_{solid}^e)_{ai}}{\hat{u}_{bj}} &= - \int\limits_{\Omega^e} \Phi_{a,k} \left(\pderiv{T_{ik}^{n+1}}{\hat{u}_{bj}}+ a_{stiff} \frac{\gamma}{\beta \Delta t} \left(\pderiv{T_{ik}^{n+1}}{\hat{u}_{bj}} \right)_{elastic} \right) dV \notag \\ &- \left( \frac{\gamma a_{mass}}{\beta \Delta t} + \frac{1}{\beta \Delta t^2} \right) \int\limits_{\Omega^e} \rho \Phi_a \Phi_c \pderiv{ \hat{u}_{ci} }{\hat{u}_{bj}}dV .$

Again, the expression for the residual equation and derivative are used to express a non-linear system of equations

$\left. \left(\pderiv{(R_{solid}^e)_{ai}}{u_{bj}} \right)\right|^{n+1}_{kiter} \left( \left. \left({u}_{bj} \right) \right|^{n+1}_{{kiter}+1} - \left. \left({u}_{bj} \right) \right|^{n+1}_{kiter} \right) = - (R_{solid})_{ai}|^{n+1}_{kiter} ,$

which are solved via the solver package. Note that the derivatives involving $u$ and $\hat{u}$ are interchangable, as are differences between the non-linear iterations.

Explicit Dynamics Time Integration (Special Implementation of Newmark Method with gamma=0.5, beta=0)

For the Newmark Method, if gamma=0.5, beta=0, and the inertial term contains a diagonalized “mass matrix”, the update equations may be carried out without the solution of a system of equations. In this case, the update equations simplify to a non-iterative update algorithm.

First the mid-step velocity and end-of-step displacements are calculated through the update equations

$\tensor{v}^{n+1/2} &= \tensor{v}^{n} + \tensor{a}^n \left( \frac{\Delta t}{2} \right), \text{ and} \\ \tensor{u}^{n+1} &= \tensor{u}^n + \tensor{v}^{n+1/2} \Delta t.$

Then the residual equation/s are calculated, and acceleration at the end-of-step is calculated via

$\left( \tensor{M} + \frac{\Delta t}{2} \tensor{C} \right) \tensor{a}^{n+1} &= \tensor{F}_{n+1} - \tensor{C} v^{n+1/2} - \tensor{K} u^{n+1} .$

Note that the mass matrix must be diagonal, and damping term may not include the stiffness based damping coefficient for this method, otherwise the above equation will require a system solve. Finally, the end-of-step velocities are calculated from the end of step acceleration:

$\tensor{v}^{n+1} &= \tensor{v}^{n+1/2} + \tensor{a}^{n+1} \left( \frac{\Delta t}{2} \right).$

Note that the velocities may be stored at the midstep, resulting one less kinematic update. This approach is typically referred to as the “Leapfrog” method. However, in GEOS we do not offer this option since it can cause some confusion that results from the storage of state at different points in time.

Parameters

In the preceding XML block, The SolidMechanics_LagrangianFEM is specified by the title of the subblock of the Solvers block. The following attributes are supported in the input block for SolidMechanics_LagrangianFEM:

XML Element: SolidMechanics_LagrangianFEM

Name	Type	Default	Description
cflFactor	real64	0.5	Factor to apply to the CFL condition when calculating the maximum allowable time step. Values should be in the interval (0,1]
contactPenaltyStiffness	real64	0	Value of the penetration penalty stiffness. Units of Pressure/length
contactRelationName	groupNameRef	NOCONTACT	Name of contact relation to enforce constraints on fracture boundary.
discretization	groupNameRef	required	Name of discretization object (defined in the Numerical Methods) to use for this solver. For instance, if this is a Finite Element Solver, the name of a Finite Element Discretization should be specified. If this is a Finite Volume Method, the name of a Finite Volume Discretization discretization should be specified.
initialDt	real64	1e+99	Initial time-step value required by the solver to the event manager.
logLevel	integer	0	Sets the level of information to write in the standard output (the console typically). Level 0 outputs no specific information for this solver. Higher levels require more outputs. 1 - Line search information - Solution information (scaling, maximum changes, quality check) - Convergence information - Time step information - Linear solver information - Nonlinear solver information - Solver timers information 2 - The summary of declared fields and coupling
massDamping	real64	0	Value of mass based damping coefficient.
maxNumResolves	integer	10	Value to indicate how many resolves may be executed after some other event is executed. For example, if a SurfaceGenerator is specified, it will be executed after the mechanics solve. However if a new surface is generated, then the mechanics solve must be executed again due to the change in topology.
name	groupName	required	A name is required for any non-unique nodes
newmarkBeta	real64	0.25	Value of $\beta$ in the Newmark Method for Implicit Dynamic time integration option. This should be pow(newmarkGamma+0.5,2.0)/4.0 unless you know what you are doing.
newmarkGamma	real64	0.5	Value of $\gamma$ in the Newmark Method for Implicit Dynamic time integration option
stiffnessDamping	real64	0	Value of stiffness based damping coefficient.
strainTheory	integer	0	Indicates whether or not to use Infinitesimal Strain Theory, or Finite Strain Theory. Valid Inputs are: 0 - Infinitesimal Strain 1 - Finite Strain
surfaceGeneratorName	string		Name of the surface generator to use
targetRegions	groupNameRef_array	required	Allowable regions that the solver may be applied to. Note that this does not indicate that the solver will be applied to these regions, only that allocation will occur such that the solver may be applied to these regions. The decision about what regions this solver will beapplied to rests in the EventManager.
timeIntegrationOption	geos_SolidMechanicsLagrangianFEM_TimeIntegrationOption	ExplicitDynamic	Time integration method. Options are: * QuasiStatic * ImplicitDynamic * ExplicitDynamic
writeLinearSystem	integer	0	Write matrix, rhs, solution to screen ( = 1) or file ( = 2).
LinearSolverParameters	node	unique	XML Element: LinearSolverParameters
NonlinearSolverParameters	node	unique	XML Element: NonlinearSolverParameters

The following data are allocated and used by the solver:

Datastructure: SolidMechanics_LagrangianFEM

Name	Type	Description
maxForce	real64	The maximum force contribution in the problem domain.
maxStableDt	real64	Value of the Maximum Stable Timestep for this solver.
meshTargets	geos_mapBase<std___1_pair<std___1_basic_string<char, std___1_char_traits<char>, std___1_allocator<char>>, std___1_basic_string<char, std___1_char_traits<char>, std___1_allocator<char>>>, LvArray_Array<std___1_basic_string<char, std___1_char_traits<char>, std___1_allocator<char>>, 1, camp_int_seq<long, 0l>, int, LvArray_ChaiBuffer>, std___1_integral_constant<bool, true>>	MeshBody/Region combinations that the solver will be applied to.
LinearSolverParameters	node	Datastructure: LinearSolverParameters
NonlinearSolverParameters	node	Datastructure: NonlinearSolverParameters
SolverStatistics	node	Datastructure: SolverStatistics

Example

An example of a valid XML block is given here:

  <Solvers>
    <SolidMechanics_LagrangianFEM
      name="lagsolve"
      strainTheory="1"
      cflFactor="0.25"
      discretization="FE1"
      targetRegions="{ Region2 }"
      />
  </Solvers>