Finite Element Discretization

Simcenter STAR-CCM+ calculates the displacement of a solid based on the principle of virtual work, which is discretized using the Finite Element method.

This approach follows the total Lagrangian displacement Finite Element formulation by Zienkiewicz and Taylor [866]. For general information on the Finite Element method, see Finite Element Method.

The continuous space domain is discretized into a finite number of elements, which are interconnected at the vertices. In each element, the nodal positions and displacements are interpolated with nodal shape functions $N_{M}$ :

Figure 1. EQUATION_DISPLAY

\begin{matrix} x & = & X + u \\ X & = & N_{M} X_{M} \\ u & = & N_{M} u_{M} \end{matrix}

(4550)

where $x$ and $X$ denotes the position vectors in the current and initial configuration, respectively. $X_{M}$ and $u_{M}$ are the position and displacement at the node M, and $N_{M}$ is a node-oriented Lagrange shape function (see H1 Lagrange Shape Functions).

The discretized form of the variation of the Green-Lagrange strain is:

Figure 2. EQUATION_DISPLAY

δ E = {\hat{B}}_{M} δ u_{M}

(4551)

where ${\hat{B}}_{M}$ is the strain-displacement matrix:

Figure 3. EQUATION_DISPLAY

{\hat{B}}_{M} = [\begin{matrix} F_{11} N_{M, 1} & F_{21} N_{M, 1} & F_{31} N_{M, 1} \\ F_{12} N_{M, 2} & F_{22} N_{M, 2} & F_{32} N_{M, 2} \\ F_{13} N_{M, 3} & F_{23} N_{M, 3} & F_{33} N_{M, 3} \\ F_{11} N_{M, 2} + F_{12} N_{M, 1} & F_{21} N_{M, 2} + F_{22} N_{M, 1} & F_{31} N_{M, 2} + F_{32} N_{M, 1} \\ F_{12} N_{M, 3} + F_{13} N_{M, 2} & F_{22} N_{M, 3} + F_{23} N_{M, 2} & F_{32} N_{M, 3} + F_{33} N_{M, 2} \\ F_{13} N_{M, 1} + F_{11} N_{M, 3} & F_{23} N_{M, 1} + F_{21} N_{M, 3} & F_{33} N_{M, 1} + F_{31} N_{M, 3} \end{matrix}]

(4552)

The deformation gradient in discretized form is:

Figure 4. EQUATION_DISPLAY

F_{i J} = δ_{i J} + u_{M i} N_{M, J} (X)

(4553)

where $N_{M, J}$ is a short hand notation for the gradients of the shape functions with respect to the initial coordinates:

Figure 5. EQUATION_DISPLAY

N_{M, J} \equiv \frac{\partial N_{M} (X)}{\partial X_{J}}

(4554)

In the infinitesimal strain approximation, the strain-displacement matrix ${\hat{B}}_{M}$ reduces to:

Figure 6. EQUATION_DISPLAY

{\hat{B}}_{M} = [\begin{matrix} N_{M, 1} & 0 & 0 \\ 0 & N_{M, 2} & 0 \\ 0 & 0 & N_{M, 3} \\ N_{M, 2} & N_{M, 1} & 0 \\ 0 & N_{M, 3} & N_{M, 2} \\ N_{M, 3} & 0 & N_{M, 1} \end{matrix}]

(4555)

Substituting the discretized displacements (Eqn. (4550)) and the discretized form of $δ E$ into the weak form (Eqn. (4464)), gives:

Figure 7. EQUATION_DISPLAY

δ Π = δ {u_{M}^{T} [\int}_{V_{0}} {\hat{B}}_{M}^{T} S d V + \int_{V_{0}} N_{M} I N_{N} ρ_{0} {\ddot{u}}_{N} d V {- \int}_{V_{0}} N_{M} b d V - \int_{Γ} N_{M} \bar{τ} d Γ - \int_{C} N_{M} \bar{q} d C - f_{M}^{p}] = 0

(4556)

where $\bar{q}$ and $f_{M}^{p}$ were introduced to take into account prescribed line loads and point forces. Since $δ u_{M}$ is zero at the Dirichlet boundaries, but otherwise arbitrary, the expression within brackets must be zero, leading to the discrete equilibrium equations:

Figure 8. EQUATION_DISPLAY

f_{M}^{i n t} + M_{M N} {\ddot{u}}_{N} = f_{M}^{e x t}

(4557)

where:

$f_{M}^{i n t}$ is the internal force at node $M$ :
Figure 9. EQUATION_DISPLAY

$\begin{array}{l} f_{M}^{i n t} = \int_{V_{0}} {\hat{B}}_{M}^{T} S d V \end{array}$
(4558)
$M_{M N} {\ddot{u}}_{N}$ is the inertial term, with the mass matrix $M_{M N}$ defined as:
Figure 10. EQUATION_DISPLAY

$\begin{matrix} M_{M N} = \int_{V_{0}} N_{M} I N_{N} ρ_{0} d V \end{matrix}$
(4559)
$\begin{matrix} f_{M}^{e x t} \end{matrix}$ is the external force applied at node $M$ :
Figure 11. EQUATION_DISPLAY

$\begin{array}{l} f_{M}^{e x t} = \int_{V_{0}} N_{M} b d V \end{array} + \int_{Γ} N_{M} \bar{τ} d Γ + \int_{C} N_{M} \bar{q} d C + f_{M}^{p}$
(4560)

which includes nodal forces resulting from prescribed body forces, surface tractions, line loads, and point forces, respectively.

Stiffness Matrix

Eqn. (4558) is linear only for the special case of infinitesimal strain, and a linear constitutive law. For large displacements, the internal forces are a nonlinear function of the displacement. The sensitivity of the internal forces with respect to the displacement is defined by the stiffness matrix:

Figure 12. EQUATION_DISPLAY

K_{M N} = \frac{\partial f_{M}^{i n t}}{\partial u_{N}}

(4561)

The stiffness matrix can be expressed as the sum of two terms, the material stiffness and the geometric stiffness:

Figure 13. EQUATION_DISPLAY

K_{M N} = K_{M N}^{m a t} + K_{M N}^{g e o m} \forall M, N \in ℕ

(4562)

The material stiffness is:

Figure 14. EQUATION_DISPLAY

K_{M N}^{m a t} = \int_{Ω} {\hat{B}}_{M}^{T} D {\hat{B}}_{N} d V

(4563)

where $D$ is the sensitivity of the 2nd Piola-Kirchoff stress tensor to the Green-Lagrange strain tensor, which is a 6x6 matrix in Voigt notation:

Figure 15. EQUATION_DISPLAY

D = \frac{\partial S}{\partial E}

(4564)

The geometric stiffness is:

Figure 16. EQUATION_DISPLAY

K_{M N}^{g e o m} = I \int_{Ω} {\frac{{\partial N}_{M}}{\partial X}}^{T} S \frac{{\partial N}_{N}}{\partial X} d V

(4565)

For linear geometry (infinitesimal strain assumption), the geometric stiffness is neglected. If the stress-strain relationship is also linear, the internal nodal forces become linear functions of the nodal displacements:

Figure 17. EQUATION_DISPLAY

f_{M}^{int} = K_{M N}^{m a t} u_{N}

(4566)

Body Load Derivatives

In general, a body load

b

can be a function of displacement, velocity, and acceleration,

b (x, u, v, a)

. The derivatives of the body load,

b (x, u, v, a)

, with respect to displacement, velocity, and acceleration can be expressed as 3 x 3 second-order tensors with components:

Figure 18. EQUATION_DISPLAY

{[D_{u} b (x, u)]}_{i j} = \frac{\partial b_{i}}{\partial u_{j}}

(4567)

Figure 19. EQUATION_DISPLAY

{[D_{v} b (x, v)]}_{i j} = \frac{\partial b_{i}}{\partial v_{j}}

(4568)

Figure 20. EQUATION_DISPLAY

{[D_{a} b (x, a)]}_{i j} = \frac{\partial b_{i}}{\partial a_{j}}

(4569)

The resultant body load vector at node

M

f_{M} \in ℝ^{n}

, for the solution-dependent body load,

b (x, u, v, a) \in ℝ^{3}

, on the domain,

V_{0}

, is assembled as:

Figure 21. EQUATION_DISPLAY

f_{M} = \int_{V_{0}}^{} N_{M}^{T} b d V

(4570)

where

N \in ℝ^{3 \times n}

are the global shape functions of the displacement field. The associated body load displacement derivative contributes to the stiffness matrix with the term:

Figure 22. EQUATION_DISPLAY

K_{M N} = \int_{V_{0}}^{} N_{M}^{T} [D_{u} b] N_{N} d V

(4571)

The symmetry of this stiffness matrix term is determined by the symmetry of the displacement derivative $D_{u} b (x, u) \in ℝ^{3 \times 3}$ .

For transient problems, the velocity and acceleration derivatives contributes to the damping and mass matrices, respectively:

Figure 23. EQUATION_DISPLAY

C_{M N} = \int_{V_{0}}^{} N_{M}^{T} [D_{v} b] N_{N} d V

(4572)

Figure 24. EQUATION_DISPLAY

M_{M N} = \int_{V_{0}}^{} N_{M}^{T} [D_{a} b] N_{N} d V

(4573)