§4 拟协调单元
一、
协调问题与拟协调单元
等参数单元有一个缺点,就是节点与插值次数无论如何增加,各单元之间都只能保证插值函数本身的连续性,这对于积分式仅含待定函数及其一阶导数的变分问题(例如空间弹性力学问题,其变形能中应变是位移的一阶偏导数)来说是适用的,但是对于包含高阶导数的变分问题(例如杆、板的弯曲问题,其变形能积分就含有待定的挠度函数的二阶偏导数),对插值函数往往要求在整个区域有一阶的连续性。而这类所谓协调问题,等参数单元是无法解决的。这种光滑性的要求当然要首先在节点上反映出来,也就是说,单元的节点参数值应当包含待定函数的有关导数值。这一类单元可统称为拟协调单元。
[型函数] 假定在某个单元上有p个节点,其局部编号为i=1,2,…,p。先考虑待定函数u=u(x,y,z)及其一阶偏导数在节点i的参数值
每一节点参数值个数r=4。于是每个单元共有4p个节点参数值。如果插值函数为三元n次多项式,则其项数或各项系数一般有个,显然次数n必须满足
这时对应于各节点参数值,可定义型函数为如下的4p个n次多项式(i=1,2,…,p) :
(i)在节点i
(ii)在其余p-1个节点j≠i,上述十六个函数即及其一阶偏导数都等于零。
(iii)
如果能把这些型函数构造出来,那末u的插值多项式就可以表示为
如果待定函数是(u,v,w),其他条件同前,则其插值函数()可表示为
或用矢量表示
[待定系数法] 待定系数法是从4p个节点参数值直接解出插值多项式的各项系数。一般地说,完全的n次多项式系数{}的个数N=会大于4p,需要对多项式补加N-4p个条件才能唯一地确定这些系数。最简单的方式是象§3中那样,限制多项式的形式(例如限制多项式为三k次多项式,即对x、y或z都是k次的,系数的个数就减为)或除去某些高次项(即令相应的系数为零)使得系数个数为4p,假定经过限制后的插值多项式改写为
从4p个节点参数值可得4p个方程
(i=1,2,…,p)
式中等分别表示由多项式各项(项数为4p)及其偏导数在节点i的数值所构成的列矢量。把这单元的全部节点参数值排成列矢量,并以U表示上式右端p个4×4p的系数矩阵依序排列所构成的4p×4p的系数矩阵,则4p个方程可简写为
从此可解出各项的系数
(15)
可以看出,在同样的限制下,型函数的各项系数实际上就是中一个分量取1,其余分量取零的解。
另一种方式是对插值函数一般表达式
按光滑性要求或物理条件附加一定的约束,并假定这些约束可表示为N-4p个关于的线性方程
(16)
式中Q为(N-4p)×N的系数矩阵。对于这完全的n次插值多项式,同样可从4p个节点参数值得出如下4p个方程
(17)
式中表示4p×N的系数矩阵,只要选取约束条件适当*,一般解联立方程组(16)(17)可得
,
把N×N矩阵G分成前4p列与后(N-4p)列二子矩阵:G=[]。则
(15)或(17)表示u的节点参数值与插值多项式的系数之间的对应关系,对待定函数v,w也有类似的关系式。如果写成矢量形式,则插值多项式为
(18)
[广义节点参数] 如果直接用插值多项式(18)代替(1)式作单元分析,则变分方程(2)应改为
式中各项{R}相当于型函数,而系数{a}起着节点参数值的作用,称为单元的广义节点参数.在上式和号后的单元系数矩阵与等价外力也是广义的,分别记作
于是有
式中A为(15)中的矩阵.对于插值多项式也有类似结果.
[节点参数值的变换] 节点参数值中含有偏导数,它们在直角坐标系和局部坐标系中的数值是不同的.例如,直角坐标系中的节点参数值
在局部坐标系中为
(20)
式中四阶矩阵右下标i表示在节点i的值.反之
记局部坐标系中的型函数为,直角坐标系中的型函数为,从表达式
可知它们具有如下的线性关系
(21)
实际上,凡是包含微分运算的部分,例如在的积分式中包含微分运算的矩阵B,边界的方向导数等等,在坐标变换下都要作相应的变换.
二、
一维单元的高次插值
[三次插值] 一元三次多项式有四项,其系数一般可由四个节点参数值来确定。现在取线段二端点为节点1,2,待定函数及其导数的节点值为节点参数值共四个,即
(i=1,2)
取距离坐标()的为局部坐标,这样线元二端点的坐标分别为=0, =1.定义满足如下条件的三次多项式为型函数:
(i) 在节点i上
(ii) 在节点j(≠i)上
(iii)
利用距离坐标的对称性从这些条件可分别定出
于是在局部坐标系中,插值多项式可表示为
如果要在直角坐标系中表示,可用坐标变换
(22)
代入上式并展开就得如下的埃尔米特三次插值多项式
[五次插值] 一元五次多项式有六项,其系数一般可由六个节点参数值来确定.现在取待定函数及其一、二阶导数的节点值为节点参数值,即
(i=1,2)
同样取距离坐标()的为局部坐标,并定义满足如下条件的五次多项式为型函数:
(i) 在节点i上
(ii) 在节点j(≠i)上及其一、二阶导数都等于零.
(iii)
从这些条件并利用距离坐标的对称性可分别定出
于是在局部坐标系中,插值多项式可表示为
如果要在直角坐标系中表示 ,同样可用坐标变换(22)代入并展开就得出如下埃尔米特五次多项式
式中在前出现因子就是由于导数在不同的坐标系中具有如下关系:
[型函数与待定系数法] 以上述五次插值多项式为例.写出五次插值函数的一般形式
式中六个待定系数可由六个边界条件
(i=1,2)
来确定,它可归结为一组关于的线性方程,其常数项就是这些节点参数值.
型函数不过是取单位矢量的特殊插值函数.例如在局部坐标中,取
则六个边界条件可表示为
由此解得
于是这个特殊的插值函数可写成
即上一小段的型函数,在局部坐标系中的其余型函数也可同样求出.
对一维的情况,用待定系数法求型函数是比较容易的.但对二维的情况,特别是对于不完全的高次插值,常常因为限制条件的补充使得型函数的定义与构成不那末简单.为了避免型函数的直接构成,可采用广义节点参数的办法.
三、
三边形单元的高次插值
局部坐标系取面积坐标,并令()=(),型函数的定义与构成只在局部坐标系中进行。
[二次插值] 二元()二次多项式或三元()二次齐次式共六项.因此需要六个节点参数值才能做到完全的二次插值.除了在三角形的顶点给函数值外,可以在三个边中点给定其法向导数值,这里表示在边=0上的外法向(图19.14).
这是最简单的拟协调板元,在六个节点各取一参数值,其型函数在局部坐标系中可定义如下:
设有六个函数与满足条件:
(i) 在顶点i上=1, =0,(i=1,2,3)
在边=0的中点上
(ii) 在其余有关节点上, ,的外法向导数都等于零.
(iii)
现在先来求在各边外法向导数的表达式.由于取()=(),而
式中表示边=0的线段长,即二顶点j,k的距离
从§2三边形单元的逆变换矩阵可得在各边的外法向导数值
式中A为三角形面积.例如,在=0边上
同理
其中表示到底边的高度.结果可统一写成
有了(23)就不难导出这几个型函数.例如,对由条件(i),(ii)可写成
式中是待定系数.令,再由于在三边的中点都等于零,可得
从而可联立解出
利用面积坐标的循环性可得.对于同样很容易得出其表达式,结果型函数及插值函数可统一写成
(i=1,2,3)
显然,这种插值函数的导数沿单元的公共边界的连续性是无法保证的,但由于它能通过所谓分片检验的收敛准则*,因而可以作为非协调板元的位移(挠度)模式.
注意,由于这些节点参数值的定义与坐标系无关,因此上述的型函数在局部的面积坐标和整体的直角坐标系中是一致的,即 (i=1,2,3).
[三次插值] 三次的二元()多项式或三元()的齐次式共十项,因此需要十个节点参数值,才能做到完全的三次插值,对每个顶点i(i=1,2,3),取节点参数值为,共九个,其余一个或取在形心O的函数值,或改为一个限制条件.对后一情况,可采用同九节点等参数单元一样的限制(§3),即要求对于三元二次齐次多项式是完全的,对于()的齐次式
这条件可表示为系数的线性方程
(24)
对十节点参数值的情况,可定义型函数(i=1,2,3)如下:
(i) 在形心O上 (i=1,2,3)
(ii) 在节点i上
(iii) 在其余节点j(≠i)上,
及其一阶偏导数都等于零.
(iv)
利用待定系数法可得局部坐标的型函数
u的插值多项式可写成
转到直角坐标系,由于节点参数值不变,再按(21)式(二维情况)可得其型函数
(i=1,2,3)
对九节点参数值加限制的情况,把形心节点与去掉,型函数定义可照搬.至于其构成则可仿照九节点等参数单元那样办法,对上述每个型函数补加一项,再要求满足限制条件(24),从而定出a .结果,可得出九个型函数
它们与直角坐标的型函数之间的关系以及插值多项式的表达式同前一情况一样,只要去掉就行了.
注意,这两种单元都不是协调的,虽然它们适用于平面弹性问题,但要作为板元还需要满足分片检验收敛的准则.因此,对九节点的三边形单元还要加上所谓平行三角剖分(即区域内各单元的三边都平行于三个固定的方向)的限制.才能作为板元.当然这种限制还可以放松到只要在Ω的各子区域作平行三角剖分就行了.
[五次插值] 完全的二元五次多项式或三元五次齐次式共有个系数,因此需要对单元给出21个插值条件.首先在每个顶点给定直到二阶导数的节点参数值:
(i=1,2,3)
共18个,其余三个条件可分别在三边上给定,使得插值函数的导数在单元之间保持连续.在边界上插值多项式的外法向导数一般是四次多项式,应由五个独立的条件来唯一确定,而节点参数值对它仅提供与在两端的数值,即只有四个独立的条件.为了在单元之间保持连续性,可对各边的加以限制.一般限制的方式有二种:
(i)
要求在边界为三次多项式,即其四次项的系数等于零.
(ii)
在边界中点的值为一
节点参数值(图19.16).
这样就在三边上得出三个方程,连同原18个方程就可以唯一确定插值多项式.下面介绍在面积坐标系中的广义节点参数法.
设()的五次齐次项排列如下:
相应的系数为
而插值多项式为
利用公式(23)对求边的外法向导数,可得
依方式(i),以代入上式并要项的系数为零,则得
式中
考虑其余二边=0,=0的限制,同理可得另二方程,经过简化,三个方程可写成如下形式
式中
对于方式(ii),则可令.代入上述的表达式可得它在中点的值;在其他二边界中点的值可同样得到.结果三个中点节点的参数值与系数之间的关系可写成
按(16),(17)的记号不难列出相应的3×21矩阵Q,而与18节点参数值相应的18×21矩阵则可写成
式中为6×6矩阵:
由于其逆矩阵
的后三行对任一方式都容易得到:先解出再代入后三式就可解出.于是得到完全五次多项式系数与等的线性关系:
左端就是广义节点参数,但对方式(ii)右端的{b}还不是节点参数值,它们要改为在三边中点的值.从(25)可知对G后三列还得分别乘上因子
注意,这时两种坐标系的节点参数值的变换已不是(20),而要改为
中点的外法向导数与坐标系无关.
四、
四边形单元
[双三次插值(协调板元)] 二元双三次插值多项式共有个系数,其一般形式可写成
只要在每个顶点给定节点参数值就可以唯一确定这些系数.
考虑局部坐标(ξ,η)中的单位矩阵,对一维单元的三次插值的型函数表达式作适当的组合,不难看出其型函数依序可写成
容易证明,它们满足如下的条件:
(i) 在顶点i上等于1.
(ii) 其余函数及其一阶与混合二阶导数在各节点上都等于零.
(iii)
于是在局部坐标系中,插值多项式可写成
注意,在两坐标系中,节点参数值之间的变换应取
而型函数相应地有如下关系
(i=1,2,3,4)
这种单元是协调单元.
[不完全的双三次插值] 节点参数值取,共有十二个插值条件,因此要唯一确定插值多项式,必须对它加以限制.一般是从上述的完全双三次多项式中去掉后面高次项,即等四项,而仍保持对称性.注意,这种插值函数虽然少了四个高次项,但仍然包含完全的三次多项式,只是由于混合导数不取作节点值,在两单元之间沿边界法向导数就无法保证其连续性。尽管如此,它还满足分片检验的收敛准则,因此可以作为板元的挠度模式.
利用待定系数法,可以求出与局部坐标()相应的型函数为
在局部坐标系中,插值多项式可写成
* 对一或二的情况应改为或.如果规定为双次多项式,则取等等.
* 注意,附加的约束应当结合单元的剖分来选定:有的约束对某种单元并不适用.例如,对三边形板元作三次插值,如果在三个顶点取函数及其一阶偏导数为节点参数值,则,于是,即需加一个条件.为了收敛性,要包含完全的二次多项式,而由于对称性,一般可令与二系数相等,但这种约束对于平行坐标轴的等腰直角三角形单元,可以论证相应的矩阵是奇异的.
* 如果在任何单元片(由若干相邻单元所组成)上,用某种插值函数进行有限元解法能得出常“应变”状态的解,则称这种插值函数能通过分片检验,而收敛性得到保证。对薄板弯曲问题,如果在任一单元片的边界上给定对应于完全二次多项式的边界条件,而用某种插值函数求解的结果,得出挠度,则称它能通过分片检验。