社区首页 >专栏 >几何非线性| 桁架单元(一)

几何非线性| 桁架单元(一)

fem178

发布于 2024-05-10 10:43:47

1790

文章被收录于专栏：数值分析与有限元编程数值分析与有限元编程

在上篇几何非线性| 应变张量，得到拉格朗日应变表达式为

$\mathbf E = \frac{1}{2}(\mathbf H + \mathbf H^T + \mathbf H^T \mathbf H)$

用指标记法

$E_{ij}=\frac{1}{2}(\frac{\partial u_i }{\partial X_j}+\frac{\partial u_j }{\partial X_i}+ \frac{\partial u_k }{\partial X_j}\frac{\partial u_k }{\partial X_i}) \quad (1)$

对于杆系结构，有

$\epsilon_{x}=\frac{\partial u }{\partial x}+\frac{1}{2}(\frac{\partial u }{\partial x})^2+ \frac{1}{2}(\frac{\partial v }{\partial x})^2+\frac{1}{2}(\frac{\partial w }{\partial x})^2 \quad (2)$

拉格朗日应变适用于描述几何非线性。

▲图1

如图1所示的桁架单元，局部坐标下的位移插值

$\begin{split} u(x) &=[1- \frac{x}{l},0,\frac{x}{l},0]\begin{Bmatrix} u_1 \\ v_1 \\ u_2 \\ v_2 \\ \end{Bmatrix} \\ v(x) &=[0,1- \frac{x}{l},0,\frac{x}{l}]\begin{Bmatrix} u_1 \\ v_1 \\ u_2 \\ v_2 \\ \end{Bmatrix} \end{split}$

$\mathbf u = \mathbf N\mathbf q^e \quad (3)$

其中，

$\mathbf q^e$

是单元节点位移矩阵。

$\begin{split} \frac{\partial u }{\partial x} &=u^{'}\\ &= \frac{1}{l}\begin{bmatrix} -1 & 0 & 1 & 0 \\ \end{bmatrix} \begin{Bmatrix} u_1 \\ v_1 \\ u_2 \\ v_2 \\ \end{Bmatrix} \\ &= \mathbf C \mathbf q^e \end{split}$

$\begin{split} \frac{\partial v }{\partial x} &=v^{'}\\ &= \frac{1}{l}\begin{bmatrix} 0 & -1 & 0 & 1 \\ \end{bmatrix} \begin{Bmatrix} u_1 \\ v_1 \\ u_2 \\ v_2 \\ \end{Bmatrix} \\ &= \mathbf D \mathbf q^e \end{split}$

拉格朗日应变

$\epsilon = \mathbf C \mathbf q^e + \frac{1}{2} {\mathbf q^e}^T {\mathbf C}^T \mathbf C \mathbf q^e +\frac{1}{2} {\mathbf q^e}^T {\mathbf D}^T \mathbf D \mathbf q^e \quad (4)$

虚位移

$\delta \mathbf u = \mathbf N \delta \mathbf q^e \quad (5)$

虚应变

$\begin{split} \delta \epsilon &= \frac{\partial \epsilon }{\partial \mathbf q^e }\delta \mathbf q^e \\ &= \mathbf B(\mathbf q^e) \delta \mathbf q^e \end{split} \quad (6)$

这里，(6)用到了变分运算公式

$\delta f = \frac{\partial f }{\partial x }\delta x + \frac{\partial f }{\partial y }\delta y$

内力虚功为

$\begin{split} \delta W_i & = \int_V \delta \boldsymbol {\epsilon}^T \boldsymbol{\sigma} dV \\ & = \delta \mathbf q^{eT} \int_V \mathbf B^T \boldsymbol {\sigma} dV \end{split} \quad (7)$

记

$f_i= \int_V \mathbf B^T \boldsymbol {\sigma} dV \quad (8)$

则

$\delta W_i = \delta \mathbf q^{eT} f_i \quad (9)$

由(6)可得

$\begin{split} \bf B &= \frac{\partial \epsilon }{\partial \mathbf q^{e}}\\ &= \frac{\partial \epsilon }{\partial u^{'}}\frac{\partial u^{'} }{\partial \mathbf q^{e}}+ \frac{\partial \epsilon }{\partial v^{'}}\frac{\partial v^{'} }{\partial \mathbf q^{e}}\\ &= (1+u^{'}) \mathbf {C} + v^{'} \mathbf D \\ &= \mathbf C + \mathbf q^{eT} \mathbf C^T\mathbf C + \mathbf q^{eT} \mathbf D^T\mathbf D \end{split} \quad (10)$

应力

$\begin{split} \sigma &= E(u^{'}+\frac{1}{2}u^{'2} + \frac{1}{2}v^{'2})\\ &= E(\mathbf C \mathbf q^e + \frac{1}{2} {\mathbf q^e}^T {\mathbf C}^T \mathbf C \mathbf q^e +\frac{1}{2} {\mathbf q^e}^T {\mathbf D}^T \mathbf D \mathbf q^e) \end{split} \quad (11)$

由(8)(10)(11)可得

$f_i= ((1+u^{'}) \mathbf C^T + v^{'} \mathbf D^T )\sigma Al \quad (12)$

▲图2

如图2所示的非线性迭代过程，当某一迭代步

$i$

达到收敛标准时，可以认为处于平衡状态，即

$f_i(\mathbf q^i) = f_e(\mathbf q^i) \quad (13)$

式中

$f_i$

是结构内力，

$f_e$

是外荷载，

$\mathbf q^i$

是

$i$

迭代步时的节点位移。

$i+1$

迭代步时的内力用一阶泰勒展开

$f_i(\mathbf q^{i+1}) \approx f_i(\mathbf q^i) + \frac{\partial f_i }{\partial \mathbf q} \Delta \mathbf q^{i} \quad (14)$

由(13)(14)得

$\frac{\partial f_i }{\partial \mathbf q} \Delta \mathbf q^{i} = f_e(\mathbf q^{i+1})-f_i(\mathbf u^i) \quad (15)$

记

$\mathbf K_T = \frac{\partial f_i }{\partial \mathbf q} \quad (16)$

其中，

$\mathbf K_T$

叫做切线刚度矩阵，(15)可写成

$\mathbf K_T \Delta \mathbf q^{i} = f_e(\mathbf q^{i+1})-f_i(\mathbf q^i) \quad (17)$

$\mathbf K_T$

是内力的导数，

$f_e(\mathbf q^{i+1})$

是新的荷载步下的外荷载。

$\begin{split} \mathbf K_T &= \int_V \mathbf B^T \frac{\partial \boldsymbol {\sigma} }{\partial \boldsymbol {\epsilon}} \frac{\partial \boldsymbol {\epsilon}} {\partial \mathbf q} dV + \int_V \frac{\partial \mathbf B^T }{\partial \mathbf q} \boldsymbol {\sigma}dV\\ &= \int_V \mathbf B^T E \mathbf B dV + \int_V \frac{\partial \mathbf B^T }{\partial \mathbf q} \boldsymbol {\sigma}dV\\ &= \mathbf K_{\mathbf q} + \mathbf K_{\boldsymbol {\sigma}}\\ \end{split} \quad (17)$

其中

$\mathbf K_{\mathbf q}$

叫做初始刚度矩阵，

$\mathbf K_{\sigma}$

叫做几何刚度矩阵。对于桁架单元

$\begin{split} \mathbf K_{\mathbf q} &= [(1+u^{'}) \mathbf {C}^T + v^{'} \mathbf D^T][(1+u^{'}) \mathbf {C} + v^{'} \mathbf D]EAl\\ &= [(1+u^{'})^2\mathbf {C}^T\mathbf {C}+ v^{'}(1+u^{'})(\mathbf {D}^T\mathbf {C}+\mathbf {C}^T\mathbf {D} ) + v^{'2}\mathbf {D}^T\mathbf {D}] EAl\\ \end{split} \quad (18)$

几何刚度矩阵

$\begin{split} \mathbf K_{\boldsymbol {\sigma}} &= \frac{\partial }{\partial u^{'} }(1+u^{'}) \mathbf {C}^T \frac{\partial u^{'} }{\partial \mathbf q } +\frac{\partial }{\partial v^{'} }(v^{'}) \mathbf {D}^T \frac{\partial v^{'} }{\partial \mathbf q } \boldsymbol {\sigma}Al \\ &= (\mathbf {C}^T\mathbf {C}+\mathbf {D}^T\mathbf {D}) \boldsymbol {\sigma}Al \end{split} \quad (19)$

其中

$\begin{split} \mathbf {C}^T\mathbf {C} &= \frac{1}{l^2} \begin{bmatrix} -1 \\ 0 \\ 1 \\ 0 \\ \end{bmatrix} \begin{bmatrix} -1 & 0 & 1 & 0 \\ \end{bmatrix} \\ &= \frac{1}{l^2} \begin{bmatrix} 1 & 0 & -1 & 0 \\ 0 & 0 & 0 & 0 \\ -1 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 \\ \end{bmatrix} \\ \end{split}$

同理

$\begin{split} \mathbf {D}^T\mathbf {D} &= \frac{1}{l^2} \begin{bmatrix} 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & -1 \\ 0 & 0 & 0 & 0 \\ 0 & -1 & 0 & 1 \\ \end{bmatrix} \\ \end{split}$

最终得到局部坐标下的切线刚度矩阵为

$\begin{split} \mathbf K_T &= \frac{EA}{l} \begin{bmatrix} (1+u^{'})^2 & (v^{'}+v^{'}u^{'}) & -(1+u^{'})^2 & -(v^{'}+v^{'}u^{'}) \\ (v^{'}+v^{'}u^{'}) & v^{'2} & -(v^{'}+v^{'}u^{'}) & -v^{'2} \\ -(1+u^{'})^2 & -(v^{'}+v^{'}u^{'}) & (1+u^{'})^2 & (v^{'}+v^{'}u^{'}) \\ -(v^{'}+v^{'}u^{'}) & -v^{'2} & (v^{'}+v^{'}u^{'}) & v^{'2} \\ \end{bmatrix} \\ &\quad + \frac{\sigma A}{l} \begin{bmatrix} 1 & 0 & -1 & 0 \\ 0 & 1 & 0 & -1 \\ -1 & 0 & 1 & 0 \\ 0 & -1 & 0 & 1 \\ \end{bmatrix} \\ \end{split} \quad (20)$

局部坐标和整体坐标下的节点位移转换关系

$\mathbf q^e = \mathbf T \mathbf q^g$

整体坐标下的切线刚度矩阵

$\mathbf K_T^g = \mathbf T^T \mathbf K_T\mathbf T$

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2024-05-06，如有侵权请联系 cloudcommunity@tencent.com 删除

partial