水工结构的三维阶谱有限元分析.doc

水工结构的三维阶谱有限元分析程 昭 陈胜宏(武汉水利电力大学水电学院)摘 要 根据p型有限单元法的阶谱特点，详细论述了三维阶谱单元法的基本的分析过程和具体的实现路径，包括基函数的构造、边界约束条件的处理、刚度矩阵和荷载列阵的形成、提高数值积分效率的途径等。引入了三维阶谱单元法的虚结点和广义结点的概念，并把三维阶谱单元法应用于水工结构计算。关键词 有限元法，阶谱单元，三维分析，p型，水工结构。本文于1999年5月24日收到，系国家自然科学基金资助项目(59979021). 有限元法是研究工程结构问题最为广泛的数值方法。从逼近真实解的途径分类，有限元法可以分为：(1)h型，即传统的有限元法，通过减小单元尺寸h来提高有限元解的精度；(2)p型，通过增加基函数多项式的阶数p来提高有限元解的精度；(3)hp型，综合了(1)(2)两种方法。文献1给出了p型有限元法的理论分析，指出p型有限元的收敛速度比h型快，而对于奇点问题，至少是h型的两倍。文献2阐述了p型有限元法中的阶谱概念及其优点。基函数具有阶谱特性的p型有限元法，称为阶谱有限单元法或阶谱单元法。当前,p型有限元法一般都采用阶谱单元，所以在不作特别说明的情况下，p型有限元法和阶谱单元法是同一个概念。文献3对一维阶谱单元法作了较为详细的研究。目前对阶谱有限单元法的研究主要是针对二维问题，且偏重于理论分析，对工程计算的具体过程如何实现论述较少。比如，约束如何体现，如何提高数值积分效率等等。本文详细地给出了三维阶谱单元法的基本的分析过程和实现路径，包括基函数的构造、边界约束条件的处理、刚度矩阵和荷载列阵的形成、提高数值积分效率的途径等。引入了三维阶谱单元法中的虚结点和广义结点的概念，并应用于水工结构工程计算，取得了满意的效果。1 三维阶谱单元的基函数1.1 基函数的特点 以Hpi表示单元尺寸不变时pi阶阶谱单元逼近空间，阶谱的要领即低阶单元逼近空间是高阶单元逼近空间的一个子集：(1)因此，低阶单元刚度矩阵是高阶单元刚度矩阵的子块，当为提高精度而升阶时，可继续利用已经计算出来的低阶单元刚度矩阵3,4。 三维阶谱单元的基函数包括点基函数、棱基函数、面基函数和体基函数。点基函数满足在本点值为1，在其余各点值为0；棱基函数满足在本棱的端点处值为0，在其余各棱上的值为0；面基函数满足在本面的棱边上的值为0，在其余各面上的值为0；体基函数满足在各面上的值为0，体基函数有时也称为内部基函数。 本文以六面体单元为例进行分析，四面体单元和楔形单元5分析过程类似。1.2 基函数的形式点基函数：(2)其中，0,0,0为i点在母单元坐标系中的坐标，0=(-1)i,0=(-1)i/2+0.5,0=(-1)i/4+0.75。这和八结点等参单元完全相同。 棱基函数(p2):(3)其中，,Lp是p阶单位Legendre多项式。 面基函数(p4)：(4)其中，i,j2,i+j=p. 体基函数(p6):(5)其中，i,j,k2,i+j+k=p. 以上相当于Serendity族基函数，共有基函数的个数：(6) 如果不限制面基函数的i+j=p和体基函数的i+j+k=p，而分别改为i,jp和i,j,kp即相当于Lagrange族基函数。Serendity族基函数比Lagrange族基函数个数少，能减轻计算量，且性能很好，本文以前者为例作分析。 其实，p可以是一切在一维阶谱单元中可行的阶谱基函数。如果让上面的p(t)=(1-t2)J(2,2)p(t),其中Jp(2,2)是p阶Jacobian单位(2,2)多项式，就可以构成基于Jacobian多项式的Serendity族和Lagrange族阶谱基函数。2 广义结点和边界约束条件 三维阶谱单元法的拓扑信息包括：(a)各结点的坐标；(b)各个单元的点编号信息；(c)各个棱的点编号信息；(d)各个面的棱编号信息；(e)各个单元的面编号信息。其中(c)-(e)可由(a)和(b)通过程序生成。由于阶谱单元法即使是粗糙的网格也能达到很高的精度，不像传统的有限元法需要精细地剖分网格，能大大减小前处理的工作量。 传统有限元法把面和棱的约束转化为面和棱上点的约束。由于阶谱单元法中面和棱有各自的基函数，因此对单元贡献了自由度，可以看成是虚结点的自由度，这些虚结点的自由度和实结点的自由度同等对待，虚结点和实结点可以统称为广义结点。这样，点基函数、棱基函数、面基函数和体基函数就统一为广义结点的基函数。有了广义结点的概念，许多问题就很容易理解。例如，固定面的约束体现为该面的四个角点为固定约束，该面的棱所代表的虚结点为固定约束，该面所代表的虚结点也为固定约束；其它约束情况可以此类推。因为基函数的阶数p小于4时没有面基函数，p小于2时没有棱基函数，故当p小于4时不考虑面的约束，当p小于2时不考虑棱的约束。3 刚度矩阵与荷载移植3.1 单元刚度矩阵 p=1时，点18分别作为第18个结点，相应的基函数为N8；p=2时，增加棱112作为第920个结点(虚结点，下面省略)，相应的基函数为E1-12；p=3时，增加棱112作为第2132个结点，相应的基函数为E3112；p=4时，增加棱112作为第3344个结点，相应的基函数为E4112，再增加面16作为第4450个结点，相应的基函数为4F(2,2)16；p=5时，增加棱112作为第5162个结点，相应的基函数为E5112，再增加面16作为第6374个结点，相应的基函数为5F(2,3)16和5F(3,2)16；p=6时，增加棱112作为第7586个结点，相应的基函数为E6112，再增加面16作为第87104个结点，相应的基函数为6F(2,4)16、6F(3,3)16和6F(4,2)16，再增加一个体结点作为第105个结点，相应的基函数为6B(2，2，2)；如此类推。基函数可以统一记为i，表示第i个广义结点的基函数，其中1ife(p).(7)坐标的插值仍同等参单元法：(8)于是(9)这样，单元刚度矩阵就是阶谱形式。3.2 整体自由度的整和 设结构的单元数、实结点个数、棱的个数和面的个数分别为N、Nn、Ne、Nf，则总的结点的个数为：(10)将结构看作点数、棱数、面数和体数分别为Nn，Ne，Nf，N的一个大的单元，将Nn，Ne，Nf，N分别对应3.1节中单元的点数8，棱数12，面数6和体数1，以3.1节相同的方法编排整体结点，并去掉约束所对应的行和列，能使总体刚度矩阵也成为阶谱形式。3.3 荷载移植 荷载移植与传统的有限元法类似：(11)需要注意的是，由于虚结点的存在，移植后的结点荷载和与原荷载不存在平衡关系，只是实结点的移植荷载和与原荷载平衡。这并不妨碍问题的解决。经3.2节方法对整体自由度进行整和后，整体荷载列阵F也是阶谱形式。4 数值积分 当阶数升高时，基函数将是高阶多项式，为保证精度，需要更多的高斯积分点，数值积分的运算量将随着阶数升高成倍增加，对三维问题尤为突出。利用阶谱单元的特点，减少运算量的措施有： (1)充分利用B的稀疏性。对于弹性问题，还可以利用D的稀疏性，将66的矩阵BiTDBj的各个元素先算出来，再对各个元素作数值积分。 (2)单元刚度矩阵的不同元素采用不同阶高斯积分来计算。实际上，如果对高阶阶谱单元，各个元素全部都用高阶高斯积分，计算时间简直难以接受。 (3)另外从基函数的构造可以看到，即使一个方向上的阶数很高，另两个方向上阶数可能很低，因此不同方向上也可以采用不同阶的高斯积分：(12)其中，r,s，设P(f,)表示多项式函数f在方向上的阶数，(13) 表1对不同积分方案的效率进行了比较，从表1可以看出，采用措施(3)能在(2)的基础上进一步大大减少运算时间，当阶数较高时，甚至能减少90%以上的运算时间。 (4)采用Howard E.Hinnant提出的矢量积分法6，也能大大减少计算时间，这里不再赘述。表1 单刚形成时间比较(单位：s)p=1p=2p=3p=4p=5方法(2)0.4259.73575.495644.4955140.450方法(3)0.4254.93024.101101.845467.128注：使用 233微机，一般编程技巧。5 应用举例 为便于分析整理成果，取段面如图1的混凝土重力坝，坝高100m，坝顶宽10m，上游坡面垂直，下游坡面斜率10.75，坝体容重2400kg/m3，弹性模量25GPa，泊松比1/6，坝基容重2600kg/m3，弹性模量75GPa，泊松比0.25.坝基上游取1.5倍的坝高，下游取2倍的坝高，坝基的深度取2倍的坝高。断面上的单元网格如图1，再取10m的坝段宽构成直四棱柱形式的三维六面体网格，一共114个单元，284个实结点。以坝踵为坐标原点，x轴指向下游，z轴垂直向上，y轴由右手法则确定。两侧坝面无y向位移，坝基上下游面无x向位移，坝基底面无z向位移。假设坝基地应力场由坝基重力场产生，坝体一次浇筑完成，水库一次库满，下游无水。图1 单元网格的侧视(实际是空间网格) 控制总能量范数误差小于10%，实际计算当p=3时即达到要求。p=1,2,3的主要数据列表比较如表2，从表1、图2和图3可以看出，整体自由度随着阶数的升高迅速增大，能量范数收敛较快，p=2和p=3时位移和能量范数已相当接近；坝踵处的应力等值线随着阶数的升高变得更密；不同的阶数情况下，应力的计算结果差别最大的地方是在坝踵处，体现了坝踵处应力集中的特点。表2 主要数据比较单元广义结点整体广义结点整体自由度总刚密度最大位移/mm能量范数能量范数误差p=182845105.786%5.992.117E+1327.04%p=22093618754.049%6.742.135E+1316.69%P=332158832403.795%6.802.141E+139.93%(a)p=1时第一主应力等值线(b)p=2时第一主应力等值线(c)p=3时第一主应力等值线图2 第一主应力(MPa)(a)坝基面上的水平应力x(b)坝基面上的垂直应力z(c)坝基面上的剪应力xz图3 坝基面上的应力(MPa)6 结 语 阶谱有限元法与传统有限元法的分析过程有一些不同的地方，主要表现在基函数的形式、边界约束条件处理、积分方案等方面。引入虚结点和广义结点的概念可以不必区分三维阶谱有限元中点棱面体等不同的拓扑量。阶谱有限元法可以在较粗糙的网格下，通过升高阶数达到较高的精度，不必重新划分网格，并且高阶可以利用低阶的计算成果，大大减轻了前处理量和计算量，自适应过程相对容易实现。阶谱有限元法适宜于处理有应力集中等奇点问题。参考文献1 I.Babuska,B.A.Szabo and I.N.Katzs.The p-version of the finite element method.SIAM J.Numer.Anal.,1981,18(3).2 O.C.Zienkiewicz,J.P.De,R.Gago and D.W.Kelly.The hierarchical concept in finite element analysis.Computers & structures,1983,(16).3 诸德超。升阶谱有限元法。北京：国防工业出版社，1993.4 O.C.Zienkiewicz,FRS and R.L.Taylor.The finite element method (Fourth edition)。Volume 1:Basic formulation and linear problems.McGraw-Hill Book Company(UK).5 P.Carnevall,R.B.Morris,Y.Tsuji and G.