3[1].4+基于薄板理论的协调板单元

等计算步骤是标准化了的，但各个公式比较冗长，这里不一一列出。 可以检验这种单元是完全满足协调性要求的。在每个边界上 w 是 5 次变化，两端结点
分别有 w，∂w ∂s ，∂2w ∂s2 ，即共有 6 个结点参数，可以唯一地确定边界上 5 次变化的 w， 所以边界上位移是协调的。另外边界上 ∂w ∂n 是 4 次变化的，两端结点分别有 ∂w ∂n 和 ∂2w ∂n ∂s ，以及边界的中结点有 ∂w ∂n ，即共有 5 个结点参数，可以唯一地确定边界上 4
151
⎧⎛ ⎪⎪⎜⎝
∂w0 ∂n
⎞ ⎟⎠4
⎫ ⎪ ⎪
⎪⎪⎛ ⎨⎪⎜⎝
∂w0 ∂n
⎞ ⎟⎠5
⎪⎪ ⎬ ⎪
=
[Z
]{δ
}e
⎪⎪⎩⎪⎛⎜⎝
∂w0 ∂n
⎞ ⎟⎠6
⎪ ⎪ ⎪⎭
式中下标 4，5，6 分别表示边 23，31，12 的中点。
按原来位移函数计算得到的各个边界两端结点 ∂w0 ∂n 的平均值可表示成
参数，但是各个边界中结点的 ∂w ∂n 值不再作为独立的参数，而是作为限制各个边界上 ∂w ∂n 为三次变化的附加条件。这样一来仍是 21 个条件用以确定位移函数中的 21 个待定
系数。计算实践表明这种单元可以取得较好的计算结果。 多结点参数的四边形协调板单元使用较少，因为只有矩形单元可以较方便地采用
并在边界 23 的中点 4 取单位值。
类似地可以有校正函数φ31 和φ12 。这样一来，则可按一定比例将它们叠加到原来非协调
的位移函数中，则有
w = w0 + γ φ1 23 + γ φ2 31 + γ φ3 12
（4.1）
其中 w0 是非协调的位移函数，即 w = [ N ]{δ }e 式。γ1 、γ 2 、γ 3 是待定常数，可以通过
下图是几种三角形协调板单元。
150
⎛ ⎜ ⎝
w,
∂w ∂y
,
−
∂w ∂x
⎞ ⎟ ⎠
⎛ ∂w ⎞ ⎜⎝ ∂
,
−
∂2w ∂n∂s
⎞ ⎟ ⎠
⎛ ⎜ ⎝
w,
∂w ∂y
,
−
∂w ∂x
,
∂2w ∂x2
,
∂2w ∂y2
,
∂2w ∂x∂y
⎞ ⎟ ⎠
一、 3 结点参数的协调元
前述的 3 结点三角形单元所以是非协调的，是由于在每个边界上法向导数 ∂w ∂n 是二
第四节 基于薄板理论的协调板单元
从上节的讨论中已知，在不少实际问题的分析中非协调单元获得较好的结果，但是收敛 性是以通过分片试验为条件的，使用范围受到限制。此外，即使收敛也并非单调的，不能对 解的上界或下界做出估计。因此在板壳有限元分析的研究中，特别在其早期，协调板单元的 研究受到相当的重视。
在经典薄板理论范围内，使板单元满足协调性要求的方法有二，一是增加结点参数，即 在结点参数中还包含 w 的二次导数项；另一是在保持每个结点有三个参数的前提下采取其 它一些措施，如附加校正函数法、再分割法等。现选择其中有代表性的一、二种加以介绍。
调整它们的大小，使得 ∂w ∂n 在各个边界中点的数值等于各个边界两端结点 ∂w ∂n 值的平
均值，也即使 w 在各个边界的法向导数 ∂w ∂n 成线性变化，因此两端的法向导数值就能唯
一地确定它了，从而使相邻单元交界面上的协调性得到实现。
按原来非协调位移函数计算得到的各个边界中点的 ∂w0 ∂n 值可表示如下：
（4.7）
上式所表示的位移函数是完全满足协调性要求的，而且对原来位移函数 w0 的完全性是
没有干扰的，因为在常应变（即常曲率和常扭率）的情况下，校正函数项恒为 0。 现在的问题是能否找到上述校正函数，回答是肯定的。例如
ε 23
=
(L1
+
L1L22 L23 L2 )(L1
+
L3 )
（4.8）
或
ε 23
以上所述两种建立协调元的方法，都是以位移 w 作为唯一场函数来设计板单元的。总 的来说位移函数比较复杂，且各自还存在一些固有的缺点，近年来的研究工作已提出很多不 限于以 w 为位移场函数的板单元。这些将在以下几节中进行讨论。
155
153
图 4.2
二、多结点参数的协调元
多结点参数是指结点参数中除 w 和它的一阶导数 ∂w ∂x ， ∂w ∂y 而外，还包含 w 的二
阶导数，甚至更高阶导数。这里仅以 21 个自由度和 18 个自由度的三角形单元为例，说明这 种单元的一些特点。
图 4.3 所示的三角形板单元，与总体坐标平行的x，y坐标系的原点放在单元的中心。位 移函数采用x，y的完全 5 次多项式，其中包含 21 个待定系数，可以用 21 个结点参数的条件
=
L1L22L23 (1+ L1) (L1 + L2 )(L1 + L3 )
（4.9）
152
可以检验，它们是满足校正函数要求的，即在全部边界上函数值等于 0；在边界 12 和
13 上一次导数等于 0；在边界 23 上 ∂ε23 ∂n 是二次变化的。现在只要令
φ23
=
ε 23 (∂ε23 ∂n)4
（4.10）
为使叠加校正函数后的位移函数在各边界中心的法向导数值等于两端结点的平均值，应 有下列方程
[Y ]{δ }e = [Z ]{δ }e +{γ }
（4.5）
从上式得到
{γ } = ([Y ] −[Z ]){δ }e
（4.6）
这样就得到经校正后的位移函数表达式
w = [N e ]{δ e} + [φ23 φ31 φ12 ]([Y ] − [Z ]){δ e}
就达到我们的目的。其中 (∂ε23 ∂n)4 是 ε23 在边界 23 中点 4 的法向导数值。类似地，还可
以得到φ31 和φ12 。
关于四边形协调板单元的构造，最简单方法是利用三角形单元的组合，如图 4.1 所示。
图 4.1
关于这类协调单元的性能，可以指出，其优点是在保持三个结点参数的条件下单元协调 性的要求可以完全满足，因而保证了有限元解的收敛性，即在单元尺寸不断减少时，解能单 调收敛于精确解。但在实际计算中单元尺寸总是有限的，因此计算结果常常使结构表现得过 于刚硬。现以简支方板承受中心集中载荷为例，采用不同的三角形单元（其中有些单元在讨 论中未涉及到）进行计算，图 4.2 上给出中心挠度的误差和网格单元数的关系。其中误差是 计算结果与解析解相差，N 是 1/2 边长的单元数。从计算结果可见，对于同是 9 个自由度的 三角形单元，非协调元（2）（3）较之协调元（5）、（6）有更好的精度的收敛性。
决定。现在的情况是每个角结点包含有 6 个参数，即wi，(∂w ∂x)i ，(∂w ∂y)i ，(∂2w ∂x2 )i ， (∂2w ∂y2 )i ，(∂2w ∂x∂y)i ，(i = 1, 2, 3) ，另外每个边界的中结点有一个参数，即 (∂w ∂n)k ， (k = 4, 5, 6) ，由此共有 21 个条件正好用以确定 21 个待定系数。单元刚度矩阵、载荷向量
154
次变化的 ∂w ∂n ，所以边界上 ∂w ∂n 也是协调的。
图 4.3 需要指出，边界中结点的出现一般是不希望的，这是因为通常都是采用直接法求解有限 元的线性代数方程组，边界中结点的出现将较多地增加方程组的带宽，在计算上是不经济的， 因此又出现 18 个自由度的板单元。它的位移函数仍是 5 次多项式，角结点仍各有 6 个结点
（4.2）
⎧⎪⎪⎛⎜⎝
∂w0 ∂n
⎞a ⎟⎠4
⎫ ⎪ ⎪
⎪⎪⎛ ⎨⎪⎜⎝
∂w0 ∂n
⎞a ⎟⎠5
⎪⎪ ⎬ ⎪
=
[Y
]{δ }e
⎪⎪⎪⎩⎛⎜⎝
∂w0 ∂n
⎞a ⎟⎠6
⎪ ⎪ ⎪⎭
（4.3）
（4.1）式中校正函数项在各个边界中点法向导数的数值是
{γ
}
=
⎧γ ⎪⎨γ
1 2
⎫ ⎪ ⎬
⎩⎪γ 3 ⎭⎪
（4.4）
次变化，只由两端结点的 ∂w ∂n 值不足以唯一地决定它。现在设想能找到分别和各个边界
相联系的校正函数，例如和边 23 相联系的校正函数φ23 ，它具有如下性质：
（1）在全部边界上φ23 = 0 ；
（2）在边界 12 和 13 上法向导数 ∂φ23 ∂n = 0 ；
（3）在边界 23 上 ∂φ23 ∂n ≠ 0 ，按二次变化；
Hermite 多项式作为位移函数，但是这种单元不适应一般几何形状的结构。而一般四边形板 单元除取三角形组合的方法而外，直接建立全单元的位移函数也比较复杂。
利用多结点参数的协调元，在某些情况下也可取得好的结果，但是由于总的自由度较多， 表达格式复杂，故计算费用较大。使其应用受到限制的另一重要原因是由于结点参数中包含 高阶导数，如果相邻单元材料性质不同或厚度不同，则保持结点位移高阶导数的连续性，就 不可能保持结点上力矩的连续性，因而不可能得到很好的计算结果。