Solution Methods

Displacement Field Solution

For infinitesimal strain (linear geometry), the internal forces are linear functions of the nodal displacements (see Eqn. (4566)). For large deformations, the internal forces are nonlinear in the displacements and Simcenter STAR-CCM+ solves the governing equations with Newton iterations:

where $\begin{array}{c}{r}_M^i\end{array}$ are the residual forces at node $M$.

Statics

The nomenclature that is used in solid mechanics is different from that used in fluid mechanics. In fluid mechanics, for transient simulations, the inertial terms are included and the solution depends on the initial conditions and time. In solid mechanics, a load can be time-dependent but the solution can still be considered static because the typical time scales to load or unload a part are very long compared to the time scales for flexural or sound waves to travel through the part. However, in truly dynamic problems, the transients due to the low-frequency flexural waves or high-frequency sound waves are important and the inertial terms must be included.

The static solution seeks the displacement field $u_M$ such that the internal forces are in equilibrium with the external forces. In static problems, the inertial terms are neglected and Eqn. (4557) reads:

$$f_M^{\mathrm{int}}=f_M^{ext}\forall M\in \mathbb{N}$$

(4595)

where $\mathbb{N}$ is the set of nodes of an element.

The residual forces are then:

$$\begin{array}{c}{r}_M^i=f_M^{ext}-f_M^{int}\end{array}$$

(4596)

Simcenter STAR-CCM+ solves Eqn. (4594) for the displacement increments $\delta u_N^i$ and updates the displacements as:

$$u_N^{i+1}=u_N^i+\delta u_N^i$$

(4597)

The iteration starts with a given initial condition $u_N^0$. For a linear problem, the solution is independent of the initial conditions. In addition, a direct solver computes the solution in one iteration.

Dynamics

The dynamic solution seeks the displacement field $u_M$ that satisfies the equation:

$$M_{MN}{\ddot{u}}_N+C_{MN}{\dot{u}}_N=f_M^{ext}-f_M^{\mathrm{int}},\forall M\in \mathbb{N}$$

(4598)

where $M_{MN}$ and $C_{MN}$ are the mass and damping matrices. The residual forces are then:

$$r_M=-M_{MN}{\ddot{u}}_N-C_{MN}{\dot{u}}_N-f_M^{\mathrm{int}}+f_M^{ext},\forall M\in \mathbb{N}$$

(4599)

Simcenter STAR-CCM+ provides two different approximations of the accelerations and velocities:

1st Order Backward Euler Method

The 1st order Backward Euler method approximates the acceleration and velocity at time-step n as:

$$\begin{array}{l}{\ddot{u}}_N^n=\frac{{\dot{u}}_N^n-{\dot{u}}_N^{n-1}}{\Delta t}\\ {\dot{u}}_N^n=\frac{u_N^n-u_N^{n-1}}{\Delta t}\end{array}$$

(4600)

This first order approximation is not recommended for high-resolution structural dynamics, as it can introduce a large amount of numerical damping. However, the numerical damping can be used to remove unwanted initial transients, or when the goal is to reach a quasi-static solution.

2nd Order Newmark Method

The 2nd order Newmark method approximates the velocity and position at the time-step n as:

$$\begin{array}{l}{\dot{u}}_N^n={\dot{u}}_N^{n-1}+\left(\gamma {\ddot{u}}_N^n+\left(1-\gamma \right){\ddot{u}}_N^{n-1}\right)\Delta t\\ u_N^n=u_N^{n-1}+{\dot{u}}_N^{n-1}\Delta t+\left(\beta {\ddot{u}}_N^n+\left(\frac{1}{2}-\beta \right){\ddot{u}}_N^{n-1}\right)\Delta t^2\end{array}$$

(4601)

On the initial time-step, the acceleration ${\ddot{u}}_N^0$ is obtained by solving Eqn. (4598) at $t=0$ :

$${\mathbf{M}}_{MN}{\ddot{\mathbf{u}}}_N^0={\mathbf{f}}_M^{ext}-{\mathbf{f}}_M^{\mathrm{int}}$$

(4602)

where the damping term ${\mathbf{C}}_{MN}{\dot{\mathbf{u}}}_N^0$ is neglected.

The method is 2nd order accurate when $\gamma =\frac{1}{2}$ and $\beta =\frac{1}{4}$. The method is absolutely stable when:

$$\gamma \ge \frac{1}{2},\beta =\frac{1}{4}{\left(\frac{1}{2}+\gamma \right)}^2$$

(4603)

Values of $\gamma >\frac{1}{2}$ introduce numerical damping, but reduces time integration to 1st order accuracy.

The effective stiffness matrix for the Newmark method is:

$${\tilde{K}}_{MN}=-\frac{\partial r_M}{\partial u_N}=+K_{MN}+\frac{M_{MN}}{\beta \Delta t^2}+\frac{\gamma C_{MN}}{\beta \Delta t}$$

(4604)

Generalized- $\alpha$ Method

The Generalized- $\alpha$ Method introduces the terms ${\alpha }_f$ and ${\alpha }_m$ , which denote an offset in the time-step between $n-1$ and $n$ .

The implementation of the Generalized- $\alpha$ time integration method in Simcenter STAR-CCM+ is based on the work of Arnold and Brüls [872]. This integration method introduces the pseudo-acceleration variable $a_N$ related to the physical acceleration ${\ddot{u}}_N$ of the system by means of the recurrence relation:

$$(1-{\alpha }_m)a_N^n+{\alpha }_m{a}_N^{n-1}=(1-{\alpha }_f){\ddot{u}}_N^n+{\alpha }_f{\ddot{u}}_N^{n-1}$$

(4605)

Where $a_N^0={\ddot{u}}_N^0$

Then, the generalized- $\alpha$ expressions for the configuration of the system $u_N$ and the velocity ${\dot{u}}_N$ is obtained by using $a_N$ in the Newmark difference formulae:

$$u_N^n=u_N^{n-1}+\mathrm{\Delta }t{\dot{u}}_N^{n-1}+{\mathrm{\Delta }t}^2(0.5-\beta )a_N^{n-1}+{\mathrm{\Delta }t}^2\beta a_N^n$$

(4606)

$${\dot{u}}_N^n={\dot{u}}_N^{n-1}+\mathrm{\Delta }t(1-\gamma )a_N^{n-1}+\gamma \mathrm{\Delta }t{a}_N^n$$

(4607)

Where $\mathrm{\Delta }t$ is the time-step.

The residual characterizing the problem is evaluated in an implicit manner in terms of the values of $u_N$ , ${\dot{u}}_N$ , and ${\ddot{u}}_N$ at the current time-step. By writing the integrator in terms of corrections to the position variable, the following parameters are introduced:

$${\beta }^{\prime }=\frac{1-{\alpha }_m}{{\mathrm{\Delta }t}^2\beta (1-{\alpha }_f)}$$

(4608)

$${\gamma }^{\prime }=\frac{\gamma }{\mathrm{\Delta }t\beta }$$

(4609)

and the tangent matrix of the system, ${\tilde{K}}_{MN}$ , takes the expression:

$${\tilde{K}}_{MN}={\beta }^{\prime }{M}_{MN}+{\gamma }^{\prime }{C}_{MN}+K_{MN}$$

(4610)

Where $M_{MN}$ is the mass matrix of the system, $C_{MN}$ is the tangent damping matrix, and $K_{MN}$ is the tangent stiffness matrix.

The scheme is unconditionally stable and second-order accurate providing the following conditions are met:

$$\gamma =\frac{1}{2}+{\alpha }_f-{\alpha }_m$$

(4611)

$$\beta =\frac{1}{4}{(1+{\alpha }_f-{\alpha }_m)}^2$$

(4612)

$${\alpha }_m\le {\alpha }_f\le \frac{1}{2}$$

(4613)

The numerical dissipation of the Generalized- $\alpha$ methods can be characterized by the spectral radius at infinity, ${\rho }_{\infty }$ . For ${\rho }_{\infty }=1$ , there is no numerical dissipation. When ${\rho }_{\infty }=0$ , the maximum possible numerical dissipation occurs (known as asymptotic annihilation). The numerical dissipation can be described by the amplitude decay factor, where $\lambda$ , which is defined in terms of the spectral radius ${\rho }_{\infty }\in [0,1]$ :

$$\lambda =\frac{1-{\rho }_{\infty }}{1+{\rho }_{\infty }}$$

(4614)

The following methods are associated to the family of the Generalized- $\alpha$ method:

Optimal Method

For the Optimal (Chung-Hulbert) time integration method, both ${\alpha }_m\ne 0$ and ${\alpha }_f\ne 0$ are non-zero and can be defined in terms of ${\rho }_{\infty }$ using [873]:

$${\alpha }_m=\frac{2{\rho }_{\infty }-1}{{\rho }_{\infty }+1}$$

(4615)

$${\alpha }_f=\frac{{\rho }_{\infty }}{{\rho }_{\infty }+1}$$

(4616)

HHT- $\alpha$ Method

For the HHT- $\alpha$ (Hilber-Hughes-Taylor) time integration method ${\alpha }_m=0$ . ${\alpha }_f$ can be defined in terms of ${\rho }_{\infty }$ using [873]:

$${\alpha }_f=\frac{1-{\rho }_{\infty }}{1+{\rho }_{\infty }}$$

(4617)

WBZ- $\alpha$ Method

For the WBZ- $\alpha$ (Wood-Bossak-Zienkiewicz) time integration method ${\alpha }_f=0$ . ${\alpha }_m$ can be defined in terms of ${\rho }_{\infty }$ using [873]:

$${\alpha }_m=\frac{{\rho }_{\infty }-1}{{\rho }_{\infty }+1}$$

(4618)

Two-Field Solution

The governing equations contain two independent variables, the displacement field $u$ and the pressure $p$ (see Nearly Incompressible Materials).

The linearized system of equations for the static problem is:

where:

$\begin{array}{c}{K}_{uu}\end{array}$ is the displacement-based stiffness matrix, $\begin{array}{c}{K}_{up}\end{array}$ is the stiffness matrix coupling the displacement and pressure fields, $\begin{array}{c}{K}_{pu}={K_{up}}^T\end{array}$, $K_{pp}$ is the pressure-based stiffness matrix, $\begin{array}{c}{\begin{array}{c}f\end{array}}_{ext}\end{array}$ is the external load vector, $\begin{array}{c}{\begin{array}{c}f\end{array}}_{\mathrm{int},u}\end{array}$ is the displacement-based internal load vector, $\begin{array}{c}{\begin{array}{c}f\end{array}}_{\mathrm{int},p}\end{array}$ is the pressure-based internal load vector, $B$ is the displacement-based strain operator, ${ℂ}_{ab}$ is the material stiffness tangent coupling fields $a$ and $b$, $N_u$ and $N_p$ are the vectors of shape functions interpolating the displacement and pressure fields, respectively. The pressure field is approximated with shape functions that are discontinuous between elements. The unknown pressure values are eliminated on the element level by static condensation.

Nonlinear Elasticity using Enhanced Assumed Strains

where $\mathrm{\Psi }$ is the strain energy potential, and $$\begin{gathered} W^{ext \\ } \end{gathered}$$ is the work done by applied body force or surface loads.