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

fem178

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

1240

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

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

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

\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所示的非线性迭代过程，当某一迭代步

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

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

式中

f_i

是结构内力，

f_e

是外荷载，

\mathbf q^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