壳单元列式
壳体插值函数
三维空间中壳单元的插值函数为 \[ \begin{aligned} \boldsymbol{x}(\xi, \eta, \zeta) &= \overline{\boldsymbol{x}}(\xi, \eta) + \boldsymbol{X}(\xi, \eta, \zeta) \\ \overline{\boldsymbol{x}}(\xi, \eta) &= \sum_{a}^{\mathrm{nen}}N^{a}(\xi, \eta) \overline{\boldsymbol{x}}^a \\ \boldsymbol{X}(\xi, \eta, \zeta) &= \sum_{a}^{\mathrm{nen}} N^{a}(\xi, \eta) \boldsymbol{X}^a(\zeta) \\ \boldsymbol{X}^a(\zeta) &= z^{a}(\zeta) \hat{\boldsymbol{X}}^a \quad(\text{不求和})\\ z^a(\zeta) &= N_{+}(\zeta) z_{+}^{a} + N_{-}(\zeta) z_{-}^{a} \\ N_{+}(\zeta) &= \frac{1}{2}(1 + \zeta), \quad N_{-}(\zeta) = \frac{1}{2}(1 - \zeta) \end{aligned} \tag{1} \] 式中, \(\overline{\boldsymbol{x}}\) 为壳单元中参考面中的点; \(\boldsymbol{X}\) 是参考平面中指向 “纤维” 方向的向量, \(\hat{\boldsymbol{X}}^a\) 是节点处纤维法向量 (注意, 纤维方向并不一定垂直于参考平面); \(z_a\) 是厚度方程. 壳单元需要输入的参数为 \[ \boxed{ \overline{\boldsymbol{x}}^a, \quad \hat{\boldsymbol{X}}^a, \quad z_{+}^{a}, \quad z_{-}^{a} } \tag{2} \]
- 纤维: 固定坐标 \((\xi, \eta)\) 得到的关于 \(\zeta\) 的参数曲线;
- lamina: 固定坐标 \(\zeta\) 得到的关于 \((\xi, \eta)\) 的参数曲面
壳单元的另一类参数: 如果给出节点纤维的上下端点 \(\boldsymbol{x}_{+}^a\), \(\boldsymbol{x}_{-}^a\), 以及参考面对应的参数 \(\overline{\zeta}\), 那么由此得到式 (2) 所需要的参数为 \[ \overline{\boldsymbol{x}}^a = \frac{1}{2}(1 + \overline{\zeta}) \boldsymbol{x}_{+}^a + \frac{1}{2}(1 - \overline{\zeta}) \boldsymbol{x}_{-}^a \\ \hat{\boldsymbol{X}}^a = \frac{\boldsymbol{x}_{+}^a - \boldsymbol{x}_{-}^a} {\| \boldsymbol{x}_{+}^a - \boldsymbol{x}_{-}^a \|} \\ \quad z_{+}^{a} = \frac{1}{2}(1 - \overline{\zeta}) \| \boldsymbol{x}_{+}^a - \boldsymbol{x}_{-}^a \| \\ \quad z_{-}^{a} = -\frac{1}{2}(1 + \overline{\zeta}) \| \boldsymbol{x}_{+}^a - \boldsymbol{x}_{-}^a \| \]
壳单元坐标系
单元积分点处 Lamina 坐标系
定义在单元积分点上的笛卡尔坐标系 \(\{ \boldsymbol{e}_1^l, ~\boldsymbol{e}_2^l, ~\boldsymbol{e}_3^l \}\), 其中 \(\boldsymbol{e}_3^l\) 垂直于过该点的 Lamina. 该坐标系需要进行构造. 首先根据单元的插值函数得到: \[ \boldsymbol{e}_{\xi} = \frac{\boldsymbol{x}_{,\xi}}{\|\boldsymbol{x}_{,\xi}\|}, \quad \boldsymbol{e}_{\eta} = \frac{\boldsymbol{x}_{,\eta}}{\|\boldsymbol{x}_{,\eta}\|} \] 这两个单位向量不一定相互垂直, 但是切于 Lamina. 因此 \(e_{3}^l\) 可通过这两个向量叉乘得到 \[ \boldsymbol{e}_{3}^{l}\triangleq \frac{\boldsymbol{e}_{\xi} \times \boldsymbol{e}_{\eta}} {\| \boldsymbol{e}_{\xi} \times \boldsymbol{e}_{\eta} \|} \] \(\boldsymbol{e}_\xi\), \(\boldsymbol{e}_\eta\) 的角平分单位向量为 \[ \boldsymbol{e}_{\alpha} = \frac{\boldsymbol{e}_\xi + \boldsymbol{e}_\eta} {\|\boldsymbol{e}_\xi + \boldsymbol{e}_\eta\|}, \quad \boldsymbol{e}_{\beta} = \frac{\boldsymbol{e}_3^l \times \boldsymbol{e}_\alpha} {\|\boldsymbol{e}_3^l \times \boldsymbol{e}_\alpha\|} \] 绕 \(\boldsymbol{e}_3^l\) 顺时针旋转 \(45\deg\) 得到 \[ \boldsymbol{e}_1^l = \frac{\sqrt{2}}{2} (\boldsymbol{e}_{\alpha} - \boldsymbol{e}_{\beta}), \quad \boldsymbol{e}_2^l = \frac{\sqrt{2}}{2} (\boldsymbol{e}_{\alpha} + \boldsymbol{e}_{\beta}) \]
节点处 fiber direction 坐标系
节点处坐标系要求 \(\boldsymbol{e}_3^f\) 与 \(\hat{\boldsymbol{X}}^a\) 重合. 其它两个方向不定, Hughes 给出一个算法, 能够在 \(\hat{\boldsymbol{X}}^a\) 接近 \(e_3\) 时 保持与全局坐标系近乎一致. 不过这类算法, 应该不是唯一的, 只要能避免 \(\hat{\boldsymbol{X}}^a\) 与一个接近平行的向量叉乘即可.
坐标变换
宏观坐标系与上述定义的坐标系之间的坐标变换关系为 \[ \boldsymbol{e}_{i}^{\cdot} = Q_{\alpha i}^{\cdot} \boldsymbol{e}_{\alpha} \] 式中, \(\cdot\) 表示积分点上的 lamina 坐标系或节点处的 fiber 坐标系. 因此, 向量 \(\boldsymbol{e}_{i}^{\cdot}\) 在宏观坐标系中的坐标列向量排布得到旋转矩阵 \([\mathrm{Q}]\): \[ [\mathrm{Q}^{\cdot}] = [\boldsymbol{e}_{1}^{\cdot}, ~\boldsymbol{e}_{2}^{\cdot}, ~\boldsymbol{e}_{3}^{\cdot}] \] 由此得到坐标变换矩阵为 \[ [\mathrm{x}^{\cdot}] = [\mathrm{Q}^{\cdot}]^{\top} [\mathrm{x}] \] lamina 和 fiber 坐标系之间的关系为 \[ \boldsymbol{e}_{i}^{f} = Q_{\alpha i}^{f} \boldsymbol{e}_{\alpha} \\ \boldsymbol{e}_{i}^{l} = Q_{\alpha i}^{l} \boldsymbol{e}_{\alpha} \\ \]
单元的运动学描述
\[ \boldsymbol{u}(\xi, \eta, \zeta) = \overline{\boldsymbol{u}}(\xi, \eta) + \boldsymbol{U}(\xi, \eta, \zeta) \\ \overline{\boldsymbol{u}}(\xi, \eta) = \sum_{a}^{\mathrm{nen}} N^{a}(\xi, \eta) \overline{\boldsymbol{u}}^a \\ \boldsymbol{U}(\xi, \eta, \zeta) = \sum_{a}^{\mathrm{nen}} N^{a}(\xi, \eta) \boldsymbol{U}^a(\zeta) \\ \boldsymbol{U}^a(\zeta) = z^{a}(\zeta) \hat{\boldsymbol{U}}^a \quad(\text{不求和}) \]
壳单元假设纤维方向只有旋转没有拉伸, 因此在小转动假设下, 节点处的向量为 \(\hat{\boldsymbol{U}}^a\) \[ \hat{\boldsymbol{U}}^a = \theta_2^a \boldsymbol{e}_{a1}^f - \theta_1^a \boldsymbol{e}_{a2}^f \]
本构方程
在壳体 Lamina 坐标系下的线弹性本构方程为 \[ \sigma^{l} = \mathbb{L}^{l} : \varepsilon^{l} \] 表示成 Voigt 记法后的矩阵形式为 \[ [\sigma^{l}] = [\mathrm{L}^{l}] : [\varepsilon^{l}] \] 壳体的 zero normal stress 假设要求应力分量 \(\sigma_{33}^{l}\) 等于 0, 由此得到关于 \(\varepsilon_{3}^{l}\) 的方程: \[ \varepsilon_{3}^{l} = -\frac{1}{L_{33}^{l}} (L_{3i'}^{l} \varepsilon_{i'}^{l}) \]
| \(i'\) | \(1\) | \(2\) | \(3\) | \(4\) | \(5\) |
|---|---|---|---|---|---|
| \(ij\) | \(11\) | \(22\) | \(12\) | \(23\) | \(13\) |
本构方程 \[ \sigma_{i'}^{l} = L_{i'j}^{l} \varepsilon_{j}^{l} = L_{i'j'}^{l} \varepsilon_{j'}^{l} + L_{i'3}^{l} \varepsilon_{3}^{l} = \left( L_{i'j'}^{l} - \frac{1}{L_{33}}^{l} L_{i'3}L_{3j'}^{l} \right) \varepsilon_{j'}^{l} \] 得到各向同性线弹性本构矩阵 \[ [\mathrm{L}^{l}] = \frac{E}{(1-\nu^2)} \begin{bmatrix} 1 & \nu & & & \\ & 1 & & & \\ && \dfrac{1-\nu}{2} && \\ &&& \kappa\dfrac{1-\nu}{2} & \\ &&&& \kappa \dfrac{1-\nu}{2} \end{bmatrix} \] 再加上 bending 修正系数 \(\kappa = 5/6\)
B 矩阵的构造
节点的自由度为 3 个位移自由度和 2 个转角自由度. 由此得到 B
矩阵的列式为 \[
[\varepsilon^{l}] = \sum_{a}^{\mathrm{nen}} [\mathrm{B}^a]
\begin{bmatrix}
\boldsymbol{u}^a \\
\theta_{a1} \\
\theta_{a2}
\end{bmatrix}
\] 式中, 矩阵 \([\mathrm{B}^a]\)
定义为 \[
[\mathrm{B}^a] = \left[\begin{array}{ll}
b_I^u & b_I^\theta
\end{array}\right]=\left[\begin{array}{cc}
b_1^u & b_1^\theta \\
b_2^u & b_2^\theta \\
b_3^u & b_3^\theta \\
\hdashline
b_4^u & b_4^\theta \\
b_5^u & b_5^\theta
\end{array}\right]
\] 式中的 \(b_{I}^u\) 和 \(b_{I}^{\theta}\) 是行向量, \[
\begin{aligned}
& b_I^u=\left\langle b_{I
m}^u\right\rangle=\left\langle\begin{array}{lll}
b_{I1}^u & b_{I2}^u & b_{I3}^u
\end{array}\right\rangle \\
& b_I^\theta=\left\langle b_{I
a}^\theta\right\rangle=\left\langle\begin{array}{lll}
b_{I 1}^\theta & b_{I2}^\theta
\end{array}\right\rangle
\end{aligned}
\] 注意到,
节点上的位移自由度是在全局坐标系下的分量, 应变是 Lamina
坐标系下的表达式, 因此在求位移导数时需要进行坐标变换.
根据壳单元的运动学方程, 得到 Lamina 坐标系下的位移导数为 \[
\frac{ \partial u_{i}^{l} }{ \partial x_{j}^{l}}
= Q_{im}^{l} \frac{ \partial u_{m} }{ \partial x_{j}^{l}}
= Q_{im}^{l} \sum_{a}^{\mathrm{nen}}\left(
\frac{ \partial N^{a} }{ \partial x_{j}^{l}} \overline{u}_m^a
+ \frac{ \partial (N^{a} z^a)}{ \partial x_{j}^{l}} (\theta_2^a
\boldsymbol{e}_{a1}^f - \theta_1^a \boldsymbol{e}_{a2}^f) \cdot
\boldsymbol{e}_{m}
\right)\\
= Q_{im}^{l} \sum_{a}^{\mathrm{nen}}\left(
\frac{ \partial N^{a} }{ \partial x_{j}^{l}} \overline{u}_m^a
- \frac{ \partial (N^{a} z^a)}{ \partial x_{j}^{l}}
\boldsymbol{e}_{am2}^f ~\theta_1^a
+ \frac{ \partial (N^{a} z^a)}{ \partial x_{j}^{l}}
\boldsymbol{e}_{am1}^f ~\theta_2^a
\right)
\] 
式中出现的关于 lamina 坐标的导数, 可以通过坐标变换转换为宏观坐标导数的计算: \[ \frac{ \partial ~}{ \partial x_{i}^{l}} = Q_{im} \frac{ \partial ~}{ \partial x_{m}} \]
刚度矩阵
\[ \underbrace{[\mathrm{k}]}_{\rm 5 nen \times 5 nen} = [[\mathrm{k}^{ab}]], \\ \underbrace{[\mathrm{k}^{ab}]}_{\rm 5 \times 5} = \int_{\square} \int_{-1}^{+1} [\mathrm{B}^a]^{\top} [\mathrm{L}^{l}] [\mathrm{B}] J ~\mathrm{d} \zeta ~\mathrm{d} \square \]
外力向量
体积力
\[ \underbrace{[f^{\mathrm{body}}]}_{\rm 5 nen \times 1} = [[f_a^{\mathrm{body}}]] \\ \underbrace{[f_a^{\mathrm{body}}]}_{\rm 5 \times 1} = \int_{\square} \int_{-1}^{+1} [\mathrm{N}^a]^{\top} [\mathrm{b}] J ~\mathrm{d} \zeta ~\mathrm{d} \square \]
式中, 矩阵 \([\mathrm{N}^a]\) 为 \[ [\mathrm{N}^a] = \left[\begin{array}{ccc:cc} N^a & 0 & 0 & -N^a z^a e_{a 12}^f & N^a z^a e_{a 11}^f \\ 0 & N^a & 0 & -N^a z^a e_{a 22}^f & N^a z^a e_{a 21}^f \\ 0 & 0 & N^a & -N^a z^a e_{a 32}^f & N^a z^a e_{a 31}^f \end{array}\right] \]
面力
\[ \underbrace{[f^{\mathrm{surf}}]}_{\rm 5 nen \times 1} = [[f_a^{\mathrm{surf}}]] \\ \underbrace{[f_a^{\mathrm{surf}}]}_{\rm 5 \times 1} = \int_{\square} [\mathrm{N}^a]^{\top} [\mathrm{t}] J_s ~\mathrm{d} \square, \quad \zeta = \begin{cases} +1 & \text{top surface} \\ -1 & \text{bottom surface} \end{cases} \]
式中, \(J_s = \| \boldsymbol{x}_{,\xi} \times \boldsymbol{x}_{,\eta} \|\).
当面力为压力时, 对应的面力为 \[ t = -\zeta p \boldsymbol{n}, \quad \boldsymbol{n} = \frac{\boldsymbol{x}_{,\xi} \times \boldsymbol{x}_{,\eta}} {\| \boldsymbol{x}_{,\xi} \times \boldsymbol{x}_{,\eta} \|} \] 当面力为壳体面内的剪切力时, 假设已知沿参数坐标线的面力分量, 由此得到 \[ t = h_{\xi} \boldsymbol{e}_{\xi} + h_{\eta} \boldsymbol{e}_{\eta} \]
边力
\[ [f^{\mathrm{edge}}] = [[f_a^{\mathrm{edge}}]] \\ [f_a^{\mathrm{edge}}] = \int_{-1}^{+1} \int_{-1}^{+1} [\mathrm{N}^a]^{\top} [\mathrm{t}] J_e ~\mathrm{d} \xi ~\mathrm{d} \zeta, \quad \eta = \{+1, -1\} \]
式中, \(J_e = \| \boldsymbol{x}_{,\xi} \times \boldsymbol{x}_{,\zeta} \|\).
应力结果
弯矩, 膜力以及横向剪力可以通过对应力分量 在 \(\zeta\) 方向积分得到