无单元法与有限元法及工程中的对比分析

无单元法与有限元法及工程中的对比分析马连军南昌大学建筑工程学院，南昌（330029）E-mail：[e-mail]摘要：本文分别介绍了三维无单元伽辽金法和非线性有限单元法及其相关理论。无单元伽辽金法是在移动最二乘法的基础上，建立全计算域的高阶连续可导的插值函数，只利用节点信息来建立离散模型，基于变分原理来推导无单元伽辽金法的平衡方程。非线性弹性有限单元法是在划分单元基础上，构造连续可导的插值函数，根据变分原理，由最小位能原理方程推导平衡方程。论文阐述了移动最小二乘法的基本原理和权函数、节点影响半径的选取问题，同时推导全量形式的非线性弹性本构关系，并且编制了三维无单元伽辽金法和非线性弹性全量有限单元法的程序，对悬臂梁算例进行分析，通过计算值与理论值比较表明，新兴的无单元伽辽金法在某种领域上，计算精度高于有限单元法，已经成为有限单元法的有力补充。关键词：移动最小二乘法，三维无单元伽辽金法，非线性本构关系，非线性弹性有限单元法，程序设计中图分类号：TU375文献标识码：A1.前言有限单元法[1]是工程数值分析领域的昀广泛、重要、成熟的数值方法,伴随计算机科学的飞速发展，越来越发挥它的强大功能，具有对复杂几何构形的适应性、各种物理问题的可应用性、建立严格理论基础上的可靠性及适合计算机实现的高效性。它建立于20世纪,为解决重大的科学工程问题作出卓越的贡献.随着科学的发展,有限单元法遇到了许多的困难问题.例如动态裂纹扩展、材料破坏及失效、应变局部化、高速撞击，以及大变形问题，利用有限单元法分析时出现网格畸变等问题。无单元法[2]是一种新兴的数值方法,它只需要节点信息。与传统的有限单元法相比，无单元法的近似函数不依赖于网格，在分析大变形问题时具有优势，同时它克服了有限元需要划分大量的单元及单元重构的缺点，减少了工作量，具有精度高、收敛快的优点。比较早的无单元法是光滑粒子动力学法[3][4]，又称SPH法。为在分析中得到更佳的结果，将移动昀小二乘法引入Galerkin法中[5]，后来经改进与保留提出无单元Galerkin法[6],即无单元伽辽金法（EFGM）。它是通过几个不相关的数值来模拟函数。无单元伽辽金法是与有限单元法相类似的无单元数值方法，它采用移动昀小二乘法拟合形函数，由拟合函数的自由度远远少于插值函数的基点数，具有较高的精度，又能保证应力场的连续性，同时该方法的前后处理比较方便。2.三维无单元伽辽金法(EFGM)2.1移动昀小二乘法的基本原理移动昀小二乘法是用加权昀小二乘法来近似模拟场函数。设场函数为()ux,其中(,,)Txxyz=为场点，即变量。给定场函数()ux在节点上的值为：()iiuxu=1,2,3i=,……,n(1)用移动昀小二乘法，场函数可近似的表示为：()hux=1()()()()mTjjjpxaxpxax==∑(2)其中()jpx—m维基函数，()jax—m维系数矩阵，对于三维情况()jpx可取：线性基()(1,,,)Tpxxyz=，二次基222()(1,,,,,,,,,)Tpxxyzxxyyyzzzx=(3)由移动昀小二乘法将（2）式改写成下列形式：1()()()nhiiiuxNxuNxu===∑(4)111()()[()()]()()()mTijjijNxpxAxBxpxAxBx−−===∑(5)1()()()()mTiiijAxxxpxpxω==−∑(6)1122()[()(),()(),()()]nnBxxxpxxxpxxxpxωωω=−−−L(7)其中()iNx—i节点的形函数在x点的值，()ixxω−—i节点的权函数在(,,)Txxyz=点的值当()ixxω−为常数时，()hux就是昀小二乘拟合。由式(5)可进一步得形函数关于变量的偏导数。上述方法计算量较大,精度不稳定,所以采用改进的正交基移动昀小二乘法,即以一组正交基函数来代替移动昀小二乘法中的基函数()jpx.在估值点x通过正交,得出一组正交基函数[7]:11121(,)()()(,)()()(,)()()(,)kkkkjjjnikijiikjnijiiqxxpxxqxxxpxqxxxxqxxαωαω−===⎧=−⎪⎪⎪⎨⎪=⎪⎪⎩∑∑∑(8)用(8)式正交基代替式(2)式中基函数()jpx,正交基昀小二乘法得以实现.2.2权函数选择权函数()ixxω−是直接关系到场函数的光滑性,计算量的大小,计算结果收敛的快慢,进而影响无单元伽辽金法的计算精度及复杂性.权函数的选择应遵循三条原则(1)权函数保持非负,(2)某节点(插值点)的权函数在该点取昀大值,向外成单调递减性,在影响半径外取零,(3)权函数在整个域内应连续可导,确保移动昀小二乘法的近似函数光滑。采用三次样条曲线权函数:23322144,()32441()()44,(1)3320,(1)irrrrxxrrrrrωω⎧−+≤⎪⎪⎪−==−+−≤⎨⎪⎪⎪⎩pf（9）imidrd=，iidxx=−,maxmiiddc=(10)其中r—x与ix距离，mid为影响域，maxd—权函数所支撑的域的昀大半径，ic为节点与其它影响域内昀临近的节点的距离。2.3三维无单元伽辽金法的原理本文以三维线性弹性为基础，来推导无单元伽辽金法的基本方程。弹性力学平衡方程为：平衡方程：,0ijjifσ+=，在Ω上（13）；力的边界条件：0ijjintσ−=，在tΓ上（14）位移边界条件：iiuu=，在uΓ上（15）式中if—域Ω中给定的体力，it—给定的面力，uΓ、tΓ—给定的面力与位移的边界，jn—边界tΓ的外法线方向余弦，iu—给定的位移，ijklD—本构张量，ijσ、ijε—应力与应变张量。对于线性弹性问题，变分原理[1]是平衡位移使系统总势能为零。设在域内有n个节点，位移列阵为：[]Tiiiiuuνω=12TTTTnuuuu⎡⎤=⎣⎦Li=1,2,…,n（16）由移动昀小二乘法拟合（4）式，近似位移函数为：()[(),(),()]()TiiidxuxxxNxuνω==（17）其中12()[(),(),()]nNxNxNxNx=Li=1，2，…,n对系统总位能施行变分原理得平衡方程为：Kuf=，12(,,)TTTnffff=L（18）111212123212nnnnnnkkkkkkKkkk=LLMMOML其中ijk为3×3阶矩阵：1()()()()()()unnTTTuuijijijimjmmIkBxDBxdkNxRNxdkNxRNx=Ω=Ω+Γ+∑∫∫∫∫∫i，j=1，2，…，n(0)()()()()()()tTTTiiiiIfBxxdNxfxdNxtxdσΩΩ=−Ω+Ω+Γ∫∫∫∫∫∫∫∫+()()()()()()utunnTuTuTttimimijjmjIkuxNxRrdkuxNxRrNxTxθαθα==Γ++∑∑∫∫i，j=1，2，…，nm=1，2，…，un式中()iNx—形函数，D—为弹性矩阵，()TiBx—为几何矩阵，k—为罚参数[8]，(0)σ－初应力，()tjTx—集中力，uθα()x－面位移约束，uθα()umx—点位移约束，R、r—相对应的转换矩阵。3.非线性弹性有限单元法3.1位移插值函数有限单元分析法是在选定单元形式后,构造位移插值函数.设位移场函数为()ux,其中(,,)xuvω=，则位移函数可表示为:()uuxvω⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦000000000000000000iiiiiiiiiuNNNvNNNNNNω⎡⎤⎡⎤⎢⎥⎢⎥==⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦×Tiiijjjmmmuvuvuvωωω⎡⎤⎣⎦TijmijmINININaaa⎡⎤⎡⎤=⎣⎦⎣⎦eNa=（19）[]Tiiiiauvω=(,,)ijm，式中iN—插值函数矩阵，ea—单元节点位移列阵。3.2应变矩阵单元位移函数确定后，将位移函数代入几何方程εLu=中，得单元应变矩阵：TeexyzxyyzzxijmLuLNaLNNNaεεεεγγγ⎡⎤⎡⎤====⎣⎦⎣⎦eeijmBBBaBa⎡⎤==⎣⎦（20）000000000000000TiiiiixyzNBLNNyxzNzyx⎡⎤∂∂∂⎢⎥∂∂∂⎢⎥⎡⎤⎢⎥∂∂∂⎢⎥==⎢⎥⎢⎥∂∂∂⎢⎥⎢⎥⎣⎦⎢⎥∂∂∂⎢⎥∂∂∂⎣⎦000000000TiiiiiiiiiNNNxyzNNNyyzNNNzyx⎡⎤∂∂∂⎢⎥∂∂∂⎢⎥⎢⎥∂∂∂=⎢⎥∂∂∂⎢⎥⎢⎥∂∂∂⎢⎥∂∂∂⎣⎦其中,,ijmBBBB⎡⎤=⎣⎦，B—应变矩阵，L—微分算子。3.3建立平衡方程根据变分原理，由昀小位能原理方程得出有限元平衡方程：Kap=（21）其上式中TeeKGkG=∑,iiijimeTjijjjmmimjmmkkkkBDBtAkkkkkk⎡⎤⎢⎥==⎢⎥⎢⎥⎣⎦()TeeepGpF=+∑,eaGa=,eeTkBDBtdΩ=Ω∫eeTfpNftdΩ=Ω∫,eeTsspNTtdsσ=∫,00eeTpBtdσσΩ=−Ω∫00eeTpBDtdεεΩ=Ω∫,eeefsPpp=+,00eeeFppσε=+K、p—结构整体刚度矩阵和结构节点荷载列阵，ek、ep—单元刚度矩阵和单元有效节点荷载列阵，a—结构节点位移列阵，G节点自由度转换矩阵。3.4非线性弹性本构矩阵非线性关系有很多种,主要有两种形式:弹性本构关系和弹塑性本构关系[9]。其中非线性弹性本构关系分为全量形式[10]和增量形式[11].在全量形式中采用割线模量,关系式较简单,适用按比例一次加载.在增量形式中,采用切线模量和耦连模量为基础，适用范围广,但计算量大.本文采用非线性弹性全量形式，用弹性模量E和泊松比ν表示本构(一般为混凝土)矩阵:1000100010001200000[]2(1)(12)1200000212000002SSSSSSSSSSSSSSSEDνννννννννννννν−⎡⎤⎢⎥−⎢⎥⎢⎥−⎢⎥−⎢⎥=⎢⎥−+⎢⎥−⎢⎥⎢⎥⎢⎥−⎢⎥⎣⎦（22）式中0(1)2SEEβ=+，22fJJβ=，D—本构矩阵，SE、Sν—即时割线模量和割线泊松，β—材料的非线性指标，22fJJ、—应力状态不变量，2fJ可由材料的破坏准则求出，Sν可由试验公式江见鲸[10]建议公式求得.4.程序的实现4.1三维无单元伽辽金法程序根据三维无单元伽辽金法的理论分析编制计算程序.对于分析过程中积分的实现,目前有三种常用的方法:(1)节点积分.(2)利用有限元网格作背景积分网格.(3)背景积分网格,本文采用此方法.程序流程如下(1)定义结构尺寸及材料参数;(2)定义弹性矩阵D;(3)建立网格节点坐标;(4)定义网格节点的影响域;(5)建立积分网格;(6)建立网格高斯点、权重及雅可比矩阵;(7)循环高斯点:(7.a)定义节点的相临高斯点,(7.b)定义所有高斯点的权函数、形函数及形函数导数,(7.c)定义几何矩阵,(7.d)组集到总刚矩阵K;(8)定义力和位移边界条件;(9)建立边界高斯点;(10)实现积分式;(11)用罚函数法处理边界条件;(12)解方程求节点位移;(13)后处理等.4.2有限单元法程序有限单元程序流程如下(1)定义结构的形状及材料参数;(2)生成单元及单元节点编号;(3)定义非线性弹性全量式本构矩阵D;(4)施加力的边界条;(5)组集刚度;(6)单元循环组集刚度;(7)将单元刚度组集进总刚;(8)施加位移边界条件;(9)求解位移;(10)后处理等.5.实例分析悬臂梁算例:如图1所示，悬臂梁长l=10米,高和宽均为1米,弹性模量E=2.1×10kpa,泊松比µ=0.167.在端固定,在自由端施加向下荷载p=1000N,不考虑自重.分别用以上两种数值分析方法计算,梁的长度方向上划分11个计算点,高和宽分别划分5个计算点.图1悬臂梁算