采用逆向运动的方法从单元节点全位移中扣除刚体位移（包括刚体平动和刚体转动）部分，以得到节点纯变形位移量。图 1所示是4节点膜单元abcd在t₀到t时段的空间运动图，在t₀时刻单元处于abcd位置并于t时刻运动到a′b′c′d′位置。

图 1

图 1 膜单元的空间运动

由于单元4个节点的实际位置组成的是空间四边形，本文采用投影处理方式将其转换为平面4边形模式（即位置模式Ⅰ），详见2.1节讨论。单元节点i（i = a, b, c, d）在任意时刻t的实际位置向量记为x_i，则单位边向量 $\boldsymbol{e}_{ij}=\frac{\boldsymbol{x}_j－\boldsymbol{x}_i}{\left|\boldsymbol{x}_j－\boldsymbol{x}_i\right|}，$ 单元平均法线向量 $\boldsymbol{n}$ 为：

其中， $\boldsymbol{n}^0_i\begin{array}{l}\left(i=1,2,3,4\right)\end{array}$ 是空间四边形的两两相邻边所在平面的单位法线向量；对应的单元平均单位法线向量为 $\boldsymbol{n}^0=\frac{\boldsymbol{n}}{\left|\boldsymbol{n}\right|}$ 。

考虑单元在垂直于 $\boldsymbol{n}$ 且通过其形心C的平面上的投影位置作为位置模式Ⅰ（如图 2虚线所示位置）。单元的形心位置为 $\boldsymbol{x}_C=\frac{1}{4}∑\limits^4_{i=1}\boldsymbol{x}_i，$ 单元节点对形心的相对位置为 $\text{d}\boldsymbol{x}_i=\boldsymbol{x}_i－\boldsymbol{x}_C$ ，则单元节点i在该投影平面上的位置为：

图 2

图 2 膜单元的投影位置

经过投影处理后，单元从空间四边形模式简化为平面四边形模式。本文采用位置模式Ⅰ（$\bar{\boldsymbol{x}}_i$, $\bar{\boldsymbol{x}}′_i$和 $\boldsymbol{n},\boldsymbol{n}′$ ）求解t₀到t时刻的节点纯变形量和单元内力。

从t₀到t时刻的运动过程中，单元经历的刚体位移包含刚体平移和刚体转动，单元各节点i（i = a, b, c, d）的总位移向量为$\boldsymbol{u}_i=\bar{\boldsymbol{x}}′_i-\bar{\boldsymbol{x}}_i$。选择节点a为参考点，则a点的位移向量u_a为单元的刚体平移，让单元按向量(-u_a)作逆向平移运动，得到各节点扣除刚体平移后的相对位移为：

扣除刚体平移后，计算单元的刚体转动，总的刚体转动向量θ可分为2个部分：

式中，θ₁表示单元平面外转动向量，θ₂表示单元平面内转动向量。

如图 3所示，单元扣除刚体平移后，处于a″b″c″d″位置，此时与t₀时刻的单元abcd并不处于同一平面内，平面外转动角度θ₁（绕转轴单位向量 $\boldsymbol{n}^0_{op}$ 从n转动到n″的角度）为：

图 3

图 3 膜单元的平面外转动

其转轴单位向量 $\boldsymbol{n}^0_{op}=\frac{\boldsymbol{n}×\boldsymbol{n}″}{\left|\boldsymbol{n}×\boldsymbol{n}″\right|}=\{l_{op} \ m_{op} \ n_{op}\}^\text{T}，$ 从而有单元平面外转动向量 $\mathit{\boldsymbol{\theta }}_1=\theta _1\boldsymbol{n}^0_{op}。$

单元节点 $i（i = a,b,c,d）$ 的平面外逆向刚体转动位移可由下式计算得到：

式中，平面外逆向转动矩阵 $\boldsymbol{R}^*_{op}(-\theta _1)=[1-\text{cos}(-\theta _1)] \boldsymbol{A}^2_{op}+\text{sin}(-\theta _1)\boldsymbol{A}_{op}，$ 边向量 $\boldsymbol{r}″_{ai }= \boldsymbol{x}″_{i }-\boldsymbol{x}″_{a }=\mathit{\boldsymbol{\bar{x}}}′_{i }-\mathit{\boldsymbol{\bar{x}}}′_{a}。$

在经历平面外的逆向转动之后，单元回到 $t_0$ 时刻所处的平面上，得到各节点坐标为：

平面内的逆向转动位移计算如图 4所示。

图 4

图 4 膜单元的平面内转动

平面内刚体转动角度 $\theta _2$ 为：

式中， $\Delta \varphi _i（i = a,b,c,d）$ 是节点面内转动的角度，根据节点与形心连线向量的转动夹角计算得到： ${|\Delta }\varphi _i|=\text{cos}^{-1}(\boldsymbol{e}_i·\boldsymbol{e}‴)，(i = a,b,c,d)，|\Delta \varphi _i|∈\left[0,\pi \right]。$ 其中， $\boldsymbol{e}_i=\frac{\boldsymbol{x}_i－\boldsymbol{x}_C}{\left|\boldsymbol{x}_i－\boldsymbol{x}_C\right|}，\boldsymbol{e}_{i‴}=\frac{\boldsymbol{x}_{i‴}－\boldsymbol{x}_C}{\left|\boldsymbol{x}_{i‴}－\boldsymbol{x}_C\right|}，$ 下标 $C$ 表示形心。

$\Delta \varphi _i$ 的正负性可根据下式判断获得：

其转轴单位向量 $\boldsymbol{n}^0_{ip}=\boldsymbol{n}^0=\frac{\boldsymbol{n}}{\left|\boldsymbol{n}\right|}=\{l_{ip}\ m_{ip}\ n_{ip}\}^\text{T}$ 即为膜单元在 $t_0$ 时刻的单位法向向量。

从而有节点 $i（i = a,b,c,d）$ 的面内转动向量为 $\Delta \mathit{\boldsymbol{\varphi }}_i=\Delta \varphi _i\boldsymbol{n}^0,\Delta \varphi _i∈\left[－\pi ,\pi \right]，$ 单元平面内转动向量 $\mathit{\boldsymbol{\theta }}_2=\theta _2\boldsymbol{n}^0,\theta _2∈\left[－\pi ,\pi \right]。$

参考式(7) 可以得到单元节点在平面内的逆向转动位移为：

式中，平面内逆向转动矩阵 $\boldsymbol{R}^*_{ip}(-\theta _2)=[1-\text{cos}(-\theta _2)]\boldsymbol{A}^2_{ip}+\text{sin}(-\theta _2)\boldsymbol{A}_{ip}，$ 边向量 $\boldsymbol{r}‴_{ai }={\boldsymbol{x}‴_{i}-}\boldsymbol{x}‴_{a}=\mathit{\boldsymbol{\bar{x}}}′_{i }-\mathit{\boldsymbol{\bar{x}}}′_{a } +\Delta \mathit{\boldsymbol{\eta }}^r_{i－op}。$

扣除单元的刚体平移、面外刚体转动和面内刚体转动后，得到单元质点的纯变形位移为：

1.2 单元节点内力计算
跳转到：本文顶部跳转到：本文顶部　1 基本公式推导　    1.1 单元节点位移计算　    1.2 单元节点内力计算　    1.3 应力转换计算　2 数值计算　    2.1 4节点膜单元的位置模式处理　    2.2 单元节点内力计算的数值积分　3 程序编制与算例　4 结论　参考文献

单元变形满足的虚功方程为：

$\begin{split} ∑\limits_i\delta (\Delta \mathit{\boldsymbol{\eta }}^d_i)^\text{T}\boldsymbol{f}_i=∫_V\delta (\Delta \epsilon )^\text{T}\mathit{\boldsymbol{\sigma }}\text{d}V \end{split}$

(12)

其中 $\boldsymbol{f}_i$ 是整体坐标系下膜单元节点 $i$ 的节点内力。

求解膜单元节点内力时，采用定义变形坐标系的方法将空间膜元问题转化到平面上，如图 5所示定义变形坐标系的原点在参考节点 $a，$ 这样节点 $a$ 的位移变为零，即可扣除两个未知变量。坐标轴 $\hat{\boldsymbol{x}}$ 的方向则可简单地定义为单元边线 $ab$ 的方向。

图 5

图 5 单元变形坐标系

变形坐标系与整体坐标系之间的转换关系为：

$\begin{split} \hat{\boldsymbol{x}}=\hat{\boldsymbol{Q}}(\boldsymbol{x}－\boldsymbol{x}_a)=\left\{\hat{x}_x\ \hat{x}_y\ \hat{x}_z\right\}^\text{T} \end{split}$

(13)

式中， $\hat{\boldsymbol{Q}}=\left\{\begin{array}{l}\hat{\boldsymbol{e}}_1&\hat{\boldsymbol{e}}_2&\hat{\boldsymbol{e}}_3\end{array}\right\}^\text{T}$ 是坐标系间转换矩阵， ${\hat{\boldsymbol{e}}_1=}\frac{\mathit{\boldsymbol{\bar{x}}}_b－\mathit{\boldsymbol{\bar{x}}}_a}{\left|\mathit{\boldsymbol{\bar{x}}}_b－\mathit{\boldsymbol{\bar{x}}}_a\right|},\hat{\boldsymbol{e}}_3=\boldsymbol{n}^0,\hat{\boldsymbol{e}}_2=\hat{\boldsymbol{e}}_3×\hat{\boldsymbol{e}}_1。$ 变形坐标系下必有 $\hat{x}_z=0，$ 因而可简写成 $\hat{\boldsymbol{x}}=\left\{\hat{x}_x\ \hat{x}_y\right\}^\text{T}。$

在变形坐标系中单元节点的变形位移向量为：

$\begin{split} \hat{\boldsymbol{u}}_i=\hat{\boldsymbol{Q}}\Delta \mathit{\boldsymbol{\eta }}^d_i=\left\{\hat{u}_i\ \hat{v}_i\ \hat{w}_i\right\}^\text{T} \left(i=a,b,c,d\right) \end{split}$

(14)

由于变形坐标系下必有 $\hat{w}_i=0\left(i=a,b,c,d\right)，$ 因而可简写成 $\hat{\boldsymbol{u}}_i=\left\{\hat{u}_i\ \hat{v}_i\right\}^\text{T}。$

引入传统有限元的形函数得到变形系下单元内任意点的变形位移 $\hat{\boldsymbol{u}}=\left[\hat{u}\ \hat{v}\right]^\text{T}$ 为：

$\begin{split} \hat{\boldsymbol{u}}=N_a\bar{u}_a+N_b\bar{u}_b+N_c\bar{u}_c+N_d\bar{u}_d \end{split}$

(15)

式中， $N_i(\hat{x},\hat{y})$ 是变形系下的单元形函数。

对于四边形等参单元，引入局部坐标系下的正方形单元坐标 $\left(\xi ,\eta \right)$ 来求解变形系下的相关变量值。单元形函数的表达式为：

$\begin{split} N_i(\xi ,\eta )=\frac{1}{4}\left(1+\xi _i\xi \right)\left(1+\eta _i\eta \right) \end{split}$

(16)

式中， $\xi _i$ 和 $\eta _i$ 是4个节点的单元坐标，值为±1。

单元坐标系和变形坐标系下的转换通过雅可比矩阵 $\boldsymbol{J}$ 来完成，其表达式详见文献[18]。

单元节点位移的组合向量可记为：

$\begin{split} \hat{\boldsymbol{u}}^0=\left\{\hat{u}_a\ \hat{v}_a\ \hat{u}_b\ \hat{v}_b\ \hat{u}_c\ \hat{v}_c\ \hat{u}_d\ \hat{v}_d\right\}^\text{T} \end{split}$

(17)

由于变形坐标系下有 $\hat{u}_a=\hat{v}_a=0，$ 因而可简写成 $\hat{\boldsymbol{u}}^*=\left\{\hat{u}_b\ \hat{v}_b\ \hat{u}_c\ \hat{v}_c\ \hat{u}_d\ \hat{v}_d\right\}^\text{T}。$

式(15) 经简化后的单元任意点变形位移为：

$\begin{split} \hat{\boldsymbol{u}}=N_b\hat{u}_b+N_c\hat{u}_c+N_d\hat{u}_d \end{split}$

(18)

由单元形函数描述的整个膜单元变形分布，可获得单元的应变分布为：

$\begin{split} \Delta \hat{\mathit{\boldsymbol{\epsilon }}}=\mathit{\boldsymbol{B\hat{u}}}=\boldsymbol{B}^*\hat{\boldsymbol{u}}^* \end{split}$

(19)

式中， $\Delta \hat{\mathit{\boldsymbol{\epsilon }}}=\left\{\Delta \hat{\epsilon }_x\ \Delta \hat{\epsilon }_y\ \Delta \hat{\gamma }_{xy}\right\}^\text{T}，$ 应变矩阵 $\boldsymbol{B}^*$ 的表达式见文献[18]。

假设材料为线弹性，在得到单元的应变矩阵后，可得到单元的应力矩阵：

$\begin{split} \Delta \hat{\mathit{\boldsymbol{\sigma }}}=\boldsymbol{D}\Delta \hat{\epsilon }=\mathit{\boldsymbol{DB}}^*\hat{\boldsymbol{u}}^* \end{split}$

(20)

得到单元的应力应变后计算单元的变形虚功为：

$\begin{split} &\ \ \ \ \ \ \ \ \delta U=∫_V\delta (\Delta \hat{\epsilon })^\text{T}(\hat{\mathit{\boldsymbol{\sigma }}}_0+\Delta \hat{\mathit{\boldsymbol{\sigma }}})\text{d}V=\\ &(\delta \hat{\boldsymbol{u}}^*)^\text{T}[∫_V(\boldsymbol{B}^*)^\text{T}\hat{\sigma }_0\text{d}V+(∫_V(\boldsymbol{B}^*)^\text{T}\mathit{\boldsymbol{DB}}^*\text{d}V)\hat{\boldsymbol{u}}^*] \end{split}$

(21)

式中， $\hat{\mathit{\boldsymbol{\sigma }}}_0$ 是单元的初始应力。

因而有节点内力向量：

$\begin{split} &\ \ \ \ \ \ \ \ \hat{\boldsymbol{f}}^*=∫_V(\boldsymbol{B}^*)^\text{T}\hat{\sigma }_0\text{d}V+(∫_V(\boldsymbol{B}^*)^\text{T}\mathit{\boldsymbol{DB}}^*\text{d}V)\hat{\boldsymbol{u}}^*=\\ &t_a\left[∫_{\hat{A}_a}(\boldsymbol{B}^*)^\text{T}\hat{\sigma }_0\text{d}\hat{A}_a+(∫_{\hat{A}_a}(\boldsymbol{B}^*)^\text{T}\mathit{\boldsymbol{DB}}^*\text{d}\hat{A}_a)\hat{\boldsymbol{u}}^*\right] \end{split}$

(22)

其中 $\hat{\boldsymbol{f}}^*=\left\{\hat{f}_{bx}\ \hat{f}_{by}\ \hat{f}_{cx}\ \hat{f}_{cy}\ \hat{f}_{dx}\ \hat{f}_{dy}\right\}^\text{T}，$ $t_a$ 是单元厚度。另外两个分量 $\hat{f}_{ax},\hat{f}_{ay}$ 可由平面膜单元的静力平衡条件得到：

$\begin{split} &∑\hat{F}_{\hat{x}}=0,\hat{f}_{ax}=－(\hat{f}_{bx}+\hat{f}_{cx}+\hat{f}_{dx})\\ &∑\hat{F}_{\hat{y}}=0,\hat{f}_{ay}=－(\hat{f}_{by}+\hat{f}_{cy}+\hat{f}_{dy}) \end{split}$

(23)

上式求得的是单元在变形坐标系下虚拟状态时的节点内力分量，需将其转化为t时刻的内力。首先将变形坐标系下的单元节点内力向量扩展为三维向量形式 $\hat{\boldsymbol{f}}^*_i=\left\{\hat{f}_{ix}\ \hat{f}_{iy}\ \hat{f}_{iz}\right\}^\text{T}，(i = a,b,c,d)，$ 其中 $\hat{f}_{iz}=0$ ，再通过转换矩阵得到整体坐标系下的内力，然后经过正向运动转换（包括刚体转动和刚体平移，后者不影响转换）得到 $t$ 时刻的内力：

$\begin{split} \boldsymbol{f}_i=\left(\boldsymbol{R}_{ip}\boldsymbol{R}_{op}\right)\hat{\boldsymbol{Q}}^T\hat{\boldsymbol{f}}^*_i \left(i=a,b,c,d\right) \end{split}$

(24)

以上所求得的 $\boldsymbol{f}_i$ 是单元节点i由于单元纯变形产生的内力，反向作用于质点上即得到4节点膜单元传给质点的力 $\boldsymbol{f}_\alpha ^{\text{int}}。$

1.3 应力转换计算
跳转到：本文顶部跳转到：本文顶部　1 基本公式推导　    1.1 单元节点位移计算　    1.2 单元节点内力计算　    1.3 应力转换计算　2 数值计算　    2.1 4节点膜单元的位置模式处理　    2.2 单元节点内力计算的数值积分　3 程序编制与算例　4 结论　参考文献

式(11) 中的纯变形是在虚拟位置时的变形坐标系下得到的，而在初始时刻 $t_0$ 的初始应力 $\sigma _0$ 是整体坐标系下的应力分量，因此需要转换到虚拟状态的变形坐标系分量：

$\begin{split} \hat{\mathit{\boldsymbol{\sigma }}}_0=\hat{\boldsymbol{Q}}\mathit{\boldsymbol{\sigma }}_0\hat{\boldsymbol{Q}}^\text{T} \end{split}$

(25)

经上式处理后的应力可以与计算得到的应力增量在虚拟状态相加减，得到的总应力再转换回到整体坐标系，并经过正向运动转换（包括刚体转动和刚体平移，后者不影响转换）回到 $t$ 时刻，有：

$\begin{split} \mathit{\boldsymbol{\sigma ′}}=\mathit{\boldsymbol{\hat{Q}^\text{T}\left(\hat{\sigma }_0+\Delta \hat{\sigma }\right)\hat{Q}}}=\mathit{\boldsymbol{\sigma _0+\hat{Q}^\text{T}\Delta \hat{\sigma }\hat{Q}}} \end{split}$

(26)

$\begin{split} \mathit{\boldsymbol{\sigma }}=\boldsymbol{R}^\text{T}_p\sigma ′\boldsymbol{R}_p=\boldsymbol{R}^\text{T}_p\sigma _0\boldsymbol{R}_p+\left(\mathit{\boldsymbol{\hat{Q}R}}_p\right)^\text{T}\Delta \hat{\mathit{\boldsymbol{\sigma }}}\left(\mathit{\boldsymbol{\hat{Q}R}}_p\right) \end{split}$

(27)

式中， $\boldsymbol{R}_p=\boldsymbol{R}_{ip}\boldsymbol{R}_{op}。$

2 数值计算
跳转到：本文顶部跳转到：本文顶部　1 基本公式推导　    1.1 单元节点位移计算　    1.2 单元节点内力计算　    1.3 应力转换计算　2 数值计算　    2.1 4节点膜单元的位置模式处理　    2.2 单元节点内力计算的数值积分　3 程序编制与算例　4 结论　参考文献

2.1 4节点膜单元的位置模式处理
跳转到：本文顶部跳转到：本文顶部　1 基本公式推导　    1.1 单元节点位移计算　    1.2 单元节点内力计算　    1.3 应力转换计算　2 数值计算　    2.1 4节点膜单元的位置模式处理　    2.2 单元节点内力计算的数值积分　3 程序编制与算例　4 结论　参考文献

1）求解纯变形量和单元内力的节点位置模式Ⅰ

求解从 $t_0$ 到t时刻单元节点纯变形量和内力时， $t_0$ 和 $t$ 时刻单元节点的位置，即位置模式Ⅰ，均需在同一平面上。而4节点膜单元经历运动与变形前后，4个节点均不位于同一平面，即形成空间4边形。通过取用较小的时间步长，一般情况下其变形前后形状与对应平面均可保证不大的差异，此时可以用节点在垂直于单元法线 $(\mathit{\boldsymbol{n,n}}′)$ 且通过原单元形心 $(C,C′)$ 的平面上的投影位置 $(\mathit{\boldsymbol{\bar{x}}}_{i },\mathit{\boldsymbol{\bar{x}}}′_{i })$ 作为位置模式Ⅰ，以实现将空间4边形转换到平面4边形上。

2）中央差分式中的质点位置模式Ⅱ

利用中央差分式求解式(1) 时，各质点的位置, 即位置模式Ⅱ $(\boldsymbol{x}_{a },\boldsymbol{x}′_{a }),$ 需是唯一值。而采用平面投影方式获得的位置模式Ⅰ中质点对应的各单元节点并不重合，本文取位置模式Ⅰ的质点对应的各单元节点的位置平均值作为位置模式Ⅱ (即 $\boldsymbol{x}_a=\frac{1}{k}∑\limits^k_{i=1}\mathit{\boldsymbol{\bar{x}}}_i，\boldsymbol{x}′_{a }=\frac{1}{k}∑\limits^k_{i = 1}\mathit{\boldsymbol{\bar{x}}}′_{i })$ ，其中k为质点a相连各单元的对应节点的总个数。

2.2 单元节点内力计算的数值积分
跳转到：本文顶部跳转到：本文顶部　1 基本公式推导　    1.1 单元节点位移计算　    1.2 单元节点内力计算　    1.3 应力转换计算　2 数值计算　    2.1 4节点膜单元的位置模式处理　    2.2 单元节点内力计算的数值积分　3 程序编制与算例　4 结论　参考文献

由式(22) 求解4边形膜单元的节点内力时，由于单元内部的应力应变是随位置变化的，因而求解该式时将涉及到数值积分运算。引入传统有限元中二维高斯积分方案^[19]进行求解，表 1给出了单积分点和4积分点的积分点坐标及权系数。本文采用4积分点方案计算，可获得较好的计算结果。

$\begin{split} &\ \ \ \ \ \ \ \ \Delta \hat{\boldsymbol{f}}^*=\left[t_a∫_{\hat{A}_a}\left(\boldsymbol{B}^*\left(\hat{x},\hat{y}\right)\right)^\text{T}\mathit{\boldsymbol{DB}}^*\left(\hat{x},\hat{y}\right)\text{d}\hat{A}_a\right]\hat{\boldsymbol{u}}^*=\\ &\left[t_a∫^1_{－1}∫^1_{－1}\left(\boldsymbol{B}^*\left(\xi ,\eta \right)\right)^\text{T}\mathit{\boldsymbol{DB}}^*\left(\xi ,\eta \right)\left|\boldsymbol{J}\left(\xi ,\eta \right)\right|\text{d}\xi \text{d}\eta \right]\hat{\boldsymbol{u}}^* \end{split}$

(28)

表 1

表 1 4节点单元的数值积分 积分点数 坐标ξi 坐标ηi 权系数Hij 1 0 0 4.0 4 $±\frac{1}{\sqrt{3}}$ $±\frac{1}{\sqrt{3}}$ 1.0, 1.0, 1.0, 1.0

表 1 4节点单元的数值积分

其中， $\left|J\left(\xi ,\eta \right)\right|$ 是雅可比行列式。

令 $F\left(\xi ,\eta \right)=\left(\boldsymbol{B}^*\left(\xi ,\eta \right)\right)^\text{T}\mathit{\boldsymbol{DB}}^*\left(\xi ,\eta \right)|\boldsymbol{J}\left(\xi ,\eta \right)|，$ 则有：

$\begin{split} &\ \ \ \ \ \ \ \ \Delta \hat{\boldsymbol{f}}^*=\left[t_a∫^1_{－1}∫^1_{－1}F\left(\xi ,\eta \right)\text{d}\xi \text{d}\eta \right]\hat{\boldsymbol{u}}^*≅\\ &∑\limits^k_{i,j=1}t_aH_{ij}F\left(\xi _i,\eta _j\right) \end{split}$

(29)

其中， $k$ 是每个单元坐标方向上的积分点数。

3 程序编制与算例
跳转到：本文顶部跳转到：本文顶部　1 基本公式推导　    1.1 单元节点位移计算　    1.2 单元节点内力计算　    1.3 应力转换计算　2 数值计算　    2.1 4节点膜单元的位置模式处理　    2.2 单元节点内力计算的数值积分　3 程序编制与算例　4 结论　参考文献

基于前文推导的4节点膜单元基本理论公式，本文采用Matlab编制了4节点膜单元的向量式有限元分析程序。首先通过算例1验证理论公式和所编制程序的正确性和有效性，进而考察气枕充气和布料悬垂等膜结构大变形大转动问题。

算例1：内压下的双抛物线膜结构

图 6所示为受内压的薄膜结构，其初始形状由以下双抛物线曲面函数所确定：

$\begin{split} &z=\frac{16f}{a^2b^2}\left[\left(\frac{a}{2}\right)^2－\left(x－\frac{a}{2}\right)^2\right]\left[\left(\frac{b}{2}\right)^2－\left(y－\frac{b}{2}\right)^2\right],\\ &0≤x≤a;0≤y≤b \end{split}$

图 6

图 6 双抛物线膜结构及其网格划分

取矩形膜片边长a=b=1 m，矢高f=0.1 m，四边固定，膜材质量密度ρ=1.0×10⁶ kg/m³, 弹性模量为E=6.0 GPa，泊松比υ=0.267，膜厚度t_a=1.0×10^-3 m。薄膜结构表面受均布内压作用而向上发生膨胀变形，为尽快获得静力收敛解，内压p₀采用斜坡-平台方式（见图 6）并施加阻尼进行缓慢加载，以消除动力震荡效应。采用20×20网格划分，共400个四边形膜单元和441个质点（见图 6），时间步长取h=2.0×10^-4 s，阻尼参数取α=100。本例在薄膜内部分别施加内压p₀为0.5、1.0和1.5 MPa。

图 7给出了内压p₀=1.5 MPa作用时薄膜结构的膨胀变形图。图 8(a)为p₀=1.5 MPa时膜顶点z向位移计算的收敛过程，表现出很好的收敛性；图 8(b)为3种内压下本文计算结果与已有文献的比较，可见结果基本一致。本文3组结果与文献[20]结果的误差分别为6.0%、1.8%和0.2%，与文献[21]的误差略大。文献[20-21]的计算分别采用模拟方程法（AEM）和无网格迦辽金法（EFG），由于EFG法仅考虑平面弹性应变模型和应力-应变关系而忽略膜的曲率，AEM法与现有解析解相比更加准确^[20]。因此本文方法具有较高的计算精度。计算效率方面，本文方法通过加入阻尼对质点采用动力学过程获得最终静力收敛解，在CPU Q8200 @2.33 GHz和4 GB内存的计算机上，本算例耗时350 s。应该说明，本文方法基于动力学理论，对于静力问题，其计算速度与一般有限元法相比并无优势，而在动力问题中可体现出更高的计算效率。总体而言，本算例验证了向量式有限元4节点膜单元理论公式及所编制程序在膜结构大变形分析中的正确性和有效性。

图 7

图 7 内压1.5 MPa时变形图（细线为变形前形状）

图 8

图 8 内压下顶点位移计算结果

算例2：气枕充气问题

图 9所示为等边六边形气枕模型，由上、下两片六边形薄膜封闭而成，边长L=0.25 m，厚度t_a=1.0×10^-3 m，初始状态为平放且边界无约束，仅在一连续均匀填充气压q作用下作充气膨胀变形；气压q采用斜坡-平台方式缓慢施加，在t₀=0.05 s时达到最大值q₀=15.0 Pa。材料弹性模量E=7.0×10⁴ Pa，泊松比υ=0.3，质量密度ρ=27 kg/m³。采用4边形单元进行网格结构化划分，共2 400个4边形薄膜单元和2 401个质点（见图 9），分析时时间步长取h=5.0×10^-5 s，阻尼参数取α=5.0。

图 9

图 9 六边形气枕及其网格划分

图 10给出了六边形气枕模型几个典型时刻（t= 0.01、0.02、0.03、0.04、0.07、0.12 s）的变形图。由图 10可知，随着气压q的增大并达到最大值，气枕在t=0.03 s前处于膨胀变形逐渐增大状态，在t=0.03 s时膨胀变形已开始扩展到整个气枕并呈现大变形, 在t=0.04 s时达到最大膨胀变形，之后由于振荡效应的存在又开始来回往复变形（如t=0.07 s时）；由于气枕内压力达到饱和状态，在t=0.12 s时气枕周围边界上产生了明显的皱折效应。计算结果符合气枕充气的实际变形情况。本算例体现了本文4节点膜单元程序在膜结构大变形大转动分析中的有效性。

图 10

图 10 若干典型变形图

算例3：布料悬垂问题

以图 11所示方形布料悬垂问题为例进行分析，该布料边长L=0.5 m，厚度t_a=1.0×10^-3m，初始状态为水平放置且中点固定，仅在重力作用下作悬垂运动；材料弹性模量E=7.0×10⁴ Pa，泊松比υ=0.3，质量密度ρ=27 kg/m³。采用较为精细的40×40网格划分，共1 600个四边形薄膜单元和1 681个质点（见图 11），分析时时间步长取h=2.0×10^-4s，阻尼参数取α=5.0。不考虑布料自身的接触行为。

图 11

图 11 布料及其网格划分

图 12给出了布料在重力作用下几个典型时刻（t=0.10、0.20、0.25、0.30、0.35、0.40 s）的悬垂变形图。由图 12可知，布料在重力作用下四角首先开始出现对称性悬垂下摆，并出现4条明显的折痕线；接着布料四角下摆变形逐渐增大，各部分交界处相互靠近趋于明显；最后布料下摆变形继续增大直至角端处于最低点位置的变形状态为止。

图 12

图 12 若干典型变形图

由于本例并未考虑布料自身的接触行为，在t=0.30 s后布料各部分交界已出现非真实的相互穿透现象，需通过加入碰撞机制才能进行消除。尽管如此，本例仍体现了本文方法可有效模拟物体运动而进行结构大变形仿真运动分析。

4 结论
跳转到：本文顶部跳转到：本文顶部　1 基本公式推导　    1.1 单元节点位移计算　    1.2 单元节点内力计算　    1.3 应力转换计算　2 数值计算　    2.1 4节点膜单元的位置模式处理　    2.2 单元节点内力计算的数值积分　3 程序编制与算例　4 结论　参考文献

1) 基于向量式有限元，推导了4节点膜单元的基本公式，详细阐述了通过逆向运动处理膜单元的平面内、外刚体位移从而获得单元节点纯变形位移的过程，以及进一步通过变形坐标系获得单元节点内力的求解方法；同时对4节点膜单元的位置模式和内力计算的数值积分等问题提出了合理可行的处理方法。

2) 编制了向量式有限元4节点膜单元的计算分析程序，并通过算例分析验证了理论推导和所编制程序的正确性和有效性。进而将本文方法应用于气枕充气和布料仿真运动模拟，跟踪获得其大变形大转动全过程。

3) 向量式有限元可有效克服传统有限元中由于刚体运动及大变形引起的刚度矩阵奇异和计算不收敛等问题，在膜结构的大变形大转动分析中具有一定的优势。

参考文献跳转到：本文顶部跳转到：本文顶部　1 基本公式推导　    1.1 单元节点位移计算　    1.2 单元节点内力计算　    1.3 应力转换计算　2 数值计算　    2.1 4节点膜单元的位置模式处理　    2.2 单元节点内力计算的数值积分　3 程序编制与算例　4 结论　参考文献跳转到：本文顶部跳转到：本文顶部　1 基本公式推导　    1.1 单元节点位移计算　    1.2 单元节点内力计算　    1.3 应力转换计算　2 数值计算　    2.1 4节点膜单元的位置模式处理　    2.2 单元节点内力计算的数值积分　3 程序编制与算例　4 结论　参考文献

[1]	Bletzinger K U, Ramm E. A general finite element approach to the form finding of tensile structures by the updated reference strategy[J]. International Joumal of Space Struetures, 1999, 14(2): 131–145. DOI:10.1260/0266351991494759
[2]	Valdés J G, Miquel J, Onate E. Nonlinear finite element analysis of orthotropic and prestressed membrane structures[J]. Finite Elements in Analysis and Design, 2009, 45(6/7): 395–405.
[3]	伞冰冰, 武岳, 沈世钊. 膜结构有限元分析中的平面单元与曲面单元的比较[J]. 工程力学, 2008, 25(2): 168–173. San B B, Wu Y, Shen S Z. Comparison between plane and curved elements for the analysis of membrane structures by finite element method[J]. Engineering Mechanics, 2008, 25(2): 168–173. (in Chinese)
[4]	Li L, Volkov V. Cloth animation with adaptively refined meshes[C]//ACSC, 2005:107-114.
[5]	Pauletti R M O, Pimenta P M. The natural force density method for the shape finding of taut structures[J]. Computer Methods in Applied Mechanics and Engineering, 2008, 197(49/50): 4419–4428.
[6]	Han S E, Lee K S. A study of the stabilizing process of unstable structures by dynamic relaxation method[J]. Computers & Structures, 2003, 81: 1677–1688.
[7]	Maurin B, Motro R. Surface stress density method as a form-finding tool for tensile membranes[J]. Engineering Structures, 1998, 20(8): 712–719. DOI:10.1016/S0141-0296(97)00108-9
[8]	张华, 单建. 张拉膜结构的动力松弛法研究[J]. 应用力学学报, 2002, 19(1): 84–86. Zhang H, Shan J. Dynamic relaxation method study of membrane structures[J]. Chinese Journal of Applied Mechanics, 2002, 19(1): 84–86. (in Chinese)
[9]	刘凌霞, 宋强. 基于简化的指点-弹簧模型织物变形仿真研究[J]. 计算机仿真, 2011, 28(5): 406–409. Liu L G, Song Q. Cloth-object deformation simulation based on a simplified mass-spring model[J]. Computer Simulation, 2011, 28(5): 406–409. (in Chinese)
[10]	Natsupakpong S, Çavuȿoǧlu M C. Determination of elasticity parameters in lumped element (mass-spring) models of deformable objects[J]. Graphical Models, 2010, 72(6): 61–73. DOI:10.1016/j.gmod.2010.10.001
[11]	Bridson R, Fedkiw R, Anderson J. Robust treatment of collisions, contact and friction for cloth animation[J]. ACM Transactions on Graphics, 2002, 21(3): 594–603.
[12]	Teng J G, Chen S F, Hu J L. A finite-volume method for deformation analysis of woven fabrics[J]. International Journal for Numerical Methods in Engineering, 1999, 46(12): 2061–2098. DOI:10.1002/(ISSN)1097-0207
[13]	Ting E C, Shih C, Wang Y K. Fundamentals of a vector form intrinsic finite element. Part Ⅰ:Basic procedure and a plane frame element[J]. Journal of Mechanics, 2004, 20(2): 113–122. DOI:10.1017/S1727719100003336
[14]	Ting E C, Shih C, Wang Y K. Fundamentals of a vector form intrinsic finite element. Part Ⅱ:Plane solid elements[J]. Journal of Mechanics, 2004, 20(2): 123–132. DOI:10.1017/S1727719100003348
[15]	Ting E C, Shih C, Wang Y K. Fundamentals of a vector form intrinsic finite element. Part Ⅲ:Convected material frame and examples[J]. Journal of Mechanics, 2004, 20(2): 133–143. DOI:10.1017/S172771910000335X
[16]	Wu T Y, Wang C Y, Chuang C C, et al. Motion analysis of 3D membrane structures by a vector form intrinsic finite element[J]. Journal of the Chinese Institute of Engineers, 2007, 30(6): 961–976. DOI:10.1080/02533839.2007.9671324
[17]	Wu T Y, Ting E C. Large deflection analysis of 3D membrane structures by a 4-node quadrilateral intrinsic element[J]. Thin-Walled Structures, 2008, 46: 261–275. DOI:10.1016/j.tws.2007.08.043
[18]	江见鲸, 何放龙, 何益斌, 等. 有限元法及其应用[M]. 北京: 机械工业出版社, 2006.
[19]	王勖成. 有限单元法[M]. 北京: 清华大学出版社, 2008.
[20]	Tsiatas G C, Katsikadelis J T. Large deflection analysis of elastic space membranes[J]. International Journal for Numerical Method in Engineering, 2006, 65(2): 264–294. DOI:10.1002/(ISSN)1097-0207
[21]	Noguchi H, Kawashima T, Miyamura T. Element free analysis of shell and spatial structures[J]. International Journal for Numerical Method in Engineering, 2000, 47(6): 1215–1240. DOI:10.1002/(ISSN)1097-0207

X