有限元离散

Simcenter STAR-CCM+ 根据虚功原理计算固体的位移，即使用有限元方法进行离散化。

此方法遵循根据 Zienkiewicz 和 Taylor 的研究 [866] 得出的总拉格朗日位移有限元公式。有关有限元方法的一般信息，请参见有限元方法。

连续空间域将离散化成在节点处互连的有限数量的单元。在每个单元中，使用节点形状函数 $N_{M}$ 对节点位置和位移进行插值：

图 1. EQUATION_DISPLAY

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

(4550)

其中， $x$ 和 $X$ 分别表示当前和初始配置中的位置矢量。 $X_{M}$ 和 $u_{M}$ 为节点 M 处的位置和位移， $N_{M}$ 为面向节点的拉格朗日形状函数（请参见H1 拉格朗日形状函数）。

格林-拉格朗日应变的变体的离散化形式如下：

图 2. EQUATION_DISPLAY

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

(4551)

其中， ${\hat{B}}_{M}$ 为应变-位移矩阵：

图 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)

离散形式的变形梯度如下：

图 4. EQUATION_DISPLAY

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

(4553)

其中， $N_{M, J}$ 为形状函数相对于初始坐标的梯度的简化表示：

图 5. EQUATION_DISPLAY

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

(4554)

在无限小应变近似中，应变-位移矩阵 ${\hat{B}}_{M}$ 简化为：

图 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)

将离散位移 (Eqn. (4550)) 和 $δ E$ 的离散形式代入弱型 (Eqn. (4464)) 得出：

图 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)

其中和是为了考虑规定的线载荷和点力而引入的。 $\bar{q}$ $f_{M}^{p}$ 由于 $δ u_{M}$ 在狄利克雷边界处为零（否则为任意值），因此括号内的表达式必须等于零，从而获得离散平衡方程：

图 8. EQUATION_DISPLAY

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

(4557)

其中：

$f_{M}^{i n t}$ 为节点处的内力： $M$
图 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}$ 为惯性项，质量矩阵定义为： $M_{M N}$
图 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}$ 为施加于节点的外力： $M$
图 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)

包括由规定的体力、表面拉力、线载荷和点力分别产生的节点力。

刚度矩阵

仅对于无限小应变且满足线性本构定律的特殊情况，Eqn. (4558) 才为线性。对于大位移，内力是位移的非线性函数。内力相对于位移的灵敏度由刚度矩阵定义：

图 12. EQUATION_DISPLAY

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

(4561)

刚度矩阵可表示为材料刚度和几何刚度这两项的总和：

图 13. EQUATION_DISPLAY

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

(4562)

材料刚度如下：

图 14. EQUATION_DISPLAY

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

(4563)

其中， $D$ 为第二 Piola-Kirchoff 应力张量相对于格林-拉格朗日应变张量的灵敏度，即采用 Voigt 标记法的 6x6 矩阵：

图 15. EQUATION_DISPLAY

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

(4564)

几何刚度如下：

图 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)

对于线性几何（无限小应变假设），将忽略几何刚度。如果应力应变关系也为线性，则内部节点力将变为节点位移的线性函数：

图 17. EQUATION_DISPLAY

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

(4566)

体负载导数

一般情况下，体负载

b

可以是位移、速度和加速度的函数

b (x, u, v, a)

。体负载

b (x, u, v, a)

对位移、速度和加速度的导数可以表示为 3 x 3 个二阶张量，其分量如下：

图 18. EQUATION_DISPLAY

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

(4567)

图 19. EQUATION_DISPLAY

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

(4568)

图 20. EQUATION_DISPLAY

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

(4569)

对于解算方案相关体负载

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

，节点

M

处生成的体负载矢量

f_{M} \in ℝ^{n}

在域

V_{0}

中可以表示为：

图 21. EQUATION_DISPLAY

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

(4570)

其中

N \in ℝ^{3 \times n}

是位移场的全局形状函数。关联的体负载位移导数通过以下项对刚度矩阵计算产生影响：

图 22. EQUATION_DISPLAY

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

(4571)

此刚度矩阵项的对称性由位移导数 $D_{u} b (x, u) \in ℝ^{3 \times 3}$ 的对称性确定。

对于瞬态问题，速度和加速度导数分别对阻尼和质量矩阵的计算产生影响：

图 23. EQUATION_DISPLAY

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

(4572)

图 24. EQUATION_DISPLAY

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

(4573)