二维无网格伽辽金法简单例子程序
应用无网格伽辽金方法分析结构大变形问题
5 影响域
影响域是影响最小二乘精度和计算量的另一个重要参数，在移动最小二乘法的收敛性 证明中，当节点的密度趋向于无穷大的时候，影响半径应该趋向于零。也就是说，影响半 径是和节点密度密切相关的。而且，其选取对结果影响非常明显，尤其是对应力的影响。
结构工程师,Structural Engineers,2003. 66. 增刊/Sup,18~22
因此，式（1）可改写为
u h ( x) = ∑ ni ( x)ui*
i =1
n
(6)
这里 ni ( x) 即为 i 节点的形函数在 x 点的值
ni ( x) = ∑ p j ( x)[ A −1 ( x) B( x)] ji
j =1
m
(6a)
形函数关于坐标的偏导数为
1 −1 ni ,k ( x) = ∑ { p j ,k ( x)[ A −1 ( x) B( x)] ji + p j ( x)[ A,− k ( x ) B ( x ) + A ( x ) B, k ( x )] ji } (7)
1 引
言
无网格伽辽金方法（EFGM）是近年来兴起的一种新型数值计算方法，其基本思路是 利用移动最小二乘法，根据积分点附近一定影响范围（称作影响域）内的节点的位移，用 最小二乘插值得到积分点附近的近似位移场函数。它突破了传统有限元分析中单元网格的 限制，极大地简化了前后处理工作。并且，因为不使用单元网格，在结构发生大变形的情 况下，也不会出现网格畸变的问题，因此，近年来得到广泛的重视和迅速发展。 在无网格伽辽金方法中，移动最小二乘法（MLS）是影响计算结果的一个关键问题， 在实际计算过程中，移动最小二乘法有三个关键性的参数需要事先确定：第一个是基函数
j =1 m
二维弹性力学问题的光滑无网格伽辽金法
二维弹性力学问题的光滑无网格伽辽金法马文涛【摘要】计算效率低的问题长期阻碍着无网格伽辽金法(element-free Galerkin method,EFGM)的深入发展.为了提高EFGM的计算速度,本文提出一种求解二维弹性力学问题的光滑无网格伽辽金法.该方法在问题域内采用滑动最小二乘法(moving least square,MLS)近似、在域边界上采用线性插值建立位移场函数;基于广义梯度光滑算子得到两层嵌套光滑三角形背景网格上的光滑应变,根据广义光滑伽辽金弱形式建立系统离散方程.两层嵌套光滑三角形网格是由三角形背景网格本身以及四个等面积三角形子网格组成.为了提高方法的精度,由Richardson外推法确定两层光滑网格上的最优光滑应变.几个数值算例验证了该方法的精度和计算效率.数值结果表明,随着光滑积分网格数目的增加,光滑无网格伽辽金法的计算精度逐步接近EFGM的,但计算效率要远远高于EFGM的.另外,光滑无网格伽辽金法的边界条件可以像有限元那样直接施加.从计算精度和效率综合考虑,光滑无网格伽辽金法比EFGM具有更好的数值表现,具有十分广阔的发展空间.【期刊名称】《力学学报》【年(卷),期】2018(050)005【总页数】10页(P1115-1124)【关键词】无网格法;光滑应变;两层嵌套光滑三角形网格;计算效率【作者】马文涛【作者单位】宁夏大学数学统计学院,银川750021【正文语种】中文【中图分类】O241;O343引言为了克服有限元法面临的网格生成、重构以及其他与网格有关的困难,Belytschko 等[1]于1994年提出了无网格Galerkin法(element-free Galerkinmethod,EFGM).基于散乱数据的灵活插值、任意阶的光滑形函数、超收敛的数值结果以及在断裂力学问题中的成功表现,使EFGM很快成为了计算力学和工程领域众多研究者关注的热点,同时也掀起了无网格方法研究的热潮.到目前为止,已提出几十种无网格方法,并在偏微分方程数值解、金属冲压成形、高速冲击、裂纹动态扩展、流固耦合和局部化等诸多领域取得令人满意的成果[2-7].尽管无网格法种类繁多,但极少有数值方法能同时兼顾EFGM的灵活性、高精度和稳定性.因此,EFGM 仍然具有十分重要的研究价值和地位.然而,历经二十多年的发展,EFGM还仅仅是实验室高精度结果的生成者,并没有像有限元那样真正成为解决复杂工程问题的参与者.阻挠EFGM深入发展的根本原因是EFGM的计算效率非常低,严重影响了其处理实际问题的能力.其他类型的无网格法,如重构核质点法(reproducing kernel particle method,RKPM)[3],无网格局部Petrove-Galerkin法(meshless local Petrove-Galerkin method,MLPG)[4,8]、径向点插值无网格法(radial point interpolation meshless method,RPIM)[9]和一些结合式的无网格法[10]等,也存在着同样的问题.众所周知,EFGM采用滑动最小二乘法(moving least square method,MLS)建立近似函数及其导数,涉及矩阵求逆,计算量很大.其次,EFGM的形函数不具有Kronecher delta性质,导致本质边界条件施加困难.虽然Lagrange乘子法、耦合法和罚函数等[2,5]都能解决EFGM边界条件施加的问题,但同时也会带来诸如引入附加变量、改变系统矩阵性质等一系列新的问题,导致计算复杂度和计算时间的增加.另外,EFGM生成的是非多项式形式的近似解,要想获得稳定和高精度的弱形式数值解需要在每个背景网格中采用高阶高斯积分.Duan等[11]研究表明对于弹性力学问题,每个三角形背景积分网格中至少需要16个积分点才能使EFGM达到2阶精度.相比有限元法,EFGM的计算量要大得多.如何提高EFGM的计算效率成为近年来无网格法研究的热点和难点问题.为了提高EFGM的计算效率,Bessel等[12]提出了节点积分方案,也就是仅在近似节点上计算积分.该方法是不稳定的,在一些问题中会出现虚假的近奇异模态.Chen等[13-14]利用应变光滑技术,根据线性分片实验推导出Galerkin无网格法对应的积分约束,提出了稳定一致节点积分(stabilized conforming nodal integration,SCNI).SCNI很好地继承了节点积分的优势,同时完全避免了形成整体刚度矩阵时复杂形函数导数的计算过程.另外,SCNI能保证线性精度,而与无网格离散模式无关.为了弱化场函数近似的一致性要求,Liu及其团队[15-16]推广了应变光滑技术,提出了G空间理论和广义光滑Galerkin(GS-Galerkin)弱形式,建立了两类光滑数值方法的(weakened weak,W2)理论基础.这两类方法分别称为光滑有限元法(smoothed finite element method,S-FEM)和光滑点插值法(smoothed point interpolation method,S-PIM).S-FEM[17-20]的基本思想是:在有限元网格上创建的光滑区域内光滑化场函数导数,进而获得比传统有限元“更软”的光滑刚度矩阵,提高解的精度和收敛性.由于不需要计算形函数导数,避免了等参映射过程,因此S-FEM对畸形网格具有极强的适应性.S-PIM与S-FEM的基本原理相同,不同之处在于S-PIM使用多项式或径向基函数建立节点形函数.Liu[21]指出,S-FEM仅仅是S-PIM的特殊形式.Cui等[22]将MLS与GSGalerkin弱形式结合,提出了光滑伽辽金法.但与Galerkin无网格法相比,其计算精度较低.杜超凡等[23]采用S-PIM分析了旋转柔性梁的频率,得到固有频率的下界值.Liu等[24]在传统数值流形方法中引入光滑应变,改善其计算精度.Duan等[25-27]将SCNI中的线性积分约束推广到二阶情况,提出了二阶一致节点积分(quadratially consistent nodalintegration,QCI).QCI满足二阶精度,但为了求解形函数导数,必须在每个光滑区域内额外求解2个3×3的代数方程,计算代价并不小.Wang等[28]将三角形背景网格划分为两水平嵌套光滑区域,利用梯度光滑技术和Richardson外推法,提出了嵌套子域梯度光滑积分(NSGSI)算法.NSGSI算法将SCNI的精度提高到二阶.与高斯积分相比,NSGSI仅需要6个积分点,极大地提高了计算效率.为了得到二阶精度的数值算法,NSGSI必须使边界积分、外力项积分与刚度积分保持一致.显然,NSGSI的一致性要求大大降低了方法的灵活性.本文目的是提出一种能精确、高效地求解弹性力学问题的光滑无网格Galerkin法.基本思想是利用MLS近似位移场函数,在两层嵌套光滑子域上计算最优光滑应变,基于GS-Galerkin弱形式推导系统方程.两层嵌套网格由三角形背景积分网格以及连接三角形网格的三条边中点组成的4个三角形子网格构成.利用Richardson外推法,得到两层网格下对应的最优光滑应变.另外,为了通过线性分片试验条件,GS-Galerkin弱式要求在所有的边界上(包括自然边界和本质边界上)采用线性插值建立近似函数.因此,边界条件可以像有限元方法那样直接施加.1 MLS形函数在EFGM中,MLS方法被用于建立无网格形函数,具体形式为其中N(x)=pT(x)A−1(x)B(x)是MLS形函数;n是支持域内的节点数目;u=[u1u2L un]T为名义节点值向量,pT(x)为基函数向量,通常对于连续问题对于线弹性断裂问题其中(r,θ)为裂尖极坐标系下定义的位置参数.其他矩阵的定义分别为其中wI(x)=w(x−xI)为权函数,本文中使用3次样条权函数.2 系统离散方程2.1 光滑应变广义梯度光滑算子[13]作用于协调应变ε(x)其中为节点xc处的光滑应变,Φ为光滑函数,满足∫ΩcΦ(x;x−xc)dΩ =1. 简单起见,取其中为光滑区域Ωc的面积.将式(8)代入式(7),有其中nx和ny分别为光滑区域边界Γc上的单位外法线沿坐标轴方向的分量;u和v 分别为沿坐标轴方向的位移分量.将方程(1)代入方程(9),则光滑应变为其中称为光滑应变矩阵,具体形式为其中2.2 系统离散方程考虑二维线弹性固体力学问题,问题区域Ω内受体力b作用,边界Γt上受给定外力t 作用,边界Γu上满足u=.用光滑应变代替协调应变ε,得到GS-Galerkin弱式为将方程(1)和(10)代入方程(13)得其中为光滑刚度矩阵,由所有光滑区域组装得到.具体形式为式中,nsc为光滑区域总数.2.3 施加边界条件由以上推导可以看出,光滑应变运算不涉及任何形函数导数计算.因此,问题域内近似函数的不连续性并不会给光滑应变的计算带来任何困难.与标准Galerkin弱形式相比,形函数的一致性要求大大降低了.因此,在光滑无网格Galerkin法中,问题域内的点和边界上的点往往可以采用不同的近似形式.对于任何在问题域内的计算点均采用无网格近似函数;而对于任何在问题域边界上的计算点则采用线性插值函数.在边界上引入线性插值函数的目的是使光滑无网格法能够顺利通过线性分片试验,同时也使边界条件的施加变得和在有限元中一样容易.3 两层嵌套光滑积分网格为了计算光滑应变,光滑区域的选择至关重要.Liu等[21]基于背景积分网格本身、网格边界和网格顶点依次建立了不同类型的S-PIM.很显然,将背景三角形积分网格作为光滑区域是最简单、最直接的,不需要任何的附加工作.根据式(12)和Liu等[21]的研究结论,光滑区域内的光滑应变是个常数,导致了网格型光滑无网格法仅能达到线性精度,也就是说与传统三角形有限元法的精度相当.为了提高计算精度,进一步细分三角网格是非常有必要的.然而,细分过多的光滑子网格又会大大降低光滑无网格法计算效率.因此,寻找到最优的光滑网格细分方案是发展高效和高精度光滑无网格法的关键.本文提出两层嵌套光滑区域解决背景网格细分的问题.首先将问题域Ω离散为三角形背景积分网格.每个三角形网格称为父网格(见图1(a)).然后依次连结父网格的3条边的中点,形成4个等面积的子网格,如图1(b)所示.在每个光滑子网格边界上取一个高斯积分点,则根据方程(12),图1(a)所示的父网格对应的光滑应变矩阵为其中lJ和xJ分别是第J条边的边长和高斯积分点(即边界中点),nxJ和nyJ分别是第J条边单位外法线分别沿x,y方向的分量.组装所有父网格的光滑刚度矩阵可得由方程(17),可以很容易计算第m个子网格对应的光滑应变矩阵为其中Acm=Ac/4.组装所有子网格的光滑刚度矩阵为由图1可知每个光滑子网格的特征长度减小为父网格特征长度的一半,根据经典的Richardson外推法理论[28-29],方程(18)和(20)的最优线性组合可以生成更高精度的解.方程(21)是在同时考虑粗和细两种网格的基础上建立的,因此称为二层嵌套光滑区域.图1 两层嵌套三角形光滑网格(•场节点, 积分点,_高斯点Fig.1 Two-level nesting triangular smoothing cells(•node, grad e point,_Gauss point)在编程计算式(18)和式(20)的过程中,由于每个光滑积分区域的边界上采用1个积分点(线段中点)时,该点处的形函数值正好可以用线段两端点的形函数值的平均值表达.这样的话,每个三角形背景网格中需要3个顶点和3个线段中点,总共6个积分点即可.4 数值算例为了研究本文方法的精度,L2范数下的位移误差和能量误差分别定义为式(22)、式(23)中,ue,un分别代表位移精确解和数值解;εe,εn分别代表应变精确解和数值解.为了书写的方便,使用父网格(1个网格)、4个子网格以及嵌套网格(1个网格及其4个子网格)的3种光滑Galerkin无网格法分别简记为FSMM,SSMM和NSMM.EFGM采用拉格朗日乘子法施加位移边界条件.笔者根据Duan等[11]的研究结论,在EFGM的每个三角形背景积分网格中采用16个高斯积分点.上述4种方法,在悬臂梁、无限大开孔平板和双连拱隧道算例中均使用线性基函数(见式(2)),在含中心裂纹的无限大板算例中则使用内部扩展基函数(见式(3)).4种方法都在同一台电脑(处理器:Intel(R)Core(TM)*******************)上采用matlab语言编程实现.4.1 悬臂梁如图2所示,悬臂梁左侧固定,右侧受剪切作用.在平面应力条件下,解析解[30]为图2 悬臂梁Fig.2 The cantilever beam其中 I=D3/12.梁的材料参数取为E=3×107MPa,v=0.3,P=1000 N.在梁上分别布置17×5,33×9和65×17个节点.图3比较了3种节点分布下5种方法求解悬臂梁问题时的位移误差.可以看出,FSMM与传统三角形FEM的精度几乎一致;SSMM的精度要高于FSMM的;而NSMM与EFGM的计算精度几乎一致,且远远高于其他3种方法.图3 悬臂梁问题的位移误差比较Fig.3 Displacement error comparison for the cantilever beam problem图4则给出了5种方法的能量误差比较.从图4可以看出,无网格法的精度都要高于FEM的,NSMM比FSMM和TSMM的精度高;虽然EFGM的精度高于3种光滑无网格Galerkin法的,但其收敛率却最低.这些结论很好地证明了随着光滑区域的增加,积分点的个数也随之增加,光滑无网格Galerkin法的计算精度会逐步升高.图4 悬臂梁问题的能量误差比较Fig.4 Energy error comparison for the cantilever beam problem图5给出了5种方法的CPU时间比较.可以看出,在节点分布相同时,FSMM的计算耗时比FEM略长;NSMM与SSMM的计算耗时相当,比FSMM的要长,但要远远短于EFGM的.从精度和效率两个方面综合考虑,NSMM是5种方法中表现最好的. 图5 悬臂梁问题的CPU时间比较Fig.5 CPU time comparison for the cantilever beam problem4.2 无限大开孔平板设一无限大平板,中心开有半径为a的圆孔,在远离孔心的位置沿水平方向受σ0=1 Pa的轴向拉伸作用.该问题的解析解[30]为其中,r,θ为以孔心为原点的极坐标考虑对称性,仅取平板右上角四分之一进行数值计算,见图6.设板长、宽均为b=5 cm,圆孔半径a=1 cm,弹性模量E=30 MPa,泊松比v=0.3.在底部和左侧边界上分别给定位移边界条件uy(x,0)=0和ux(0,y)=0;在板右端(x=b)和上部(y=b)边界按精确解施加自然边界条件.在平板模型上分别布置9×9,17×17和33×33个节点.图7∼图9分别给出了3种节点分布下不同方法对应的位移误差、能量误差和花费的CPU时间.从图7∼图8可以看出,本文提出的NSMM的精度与EFGM的基本相同的,要高于其他几种方法.图9再次证明了NSMM的效率要远远高于EFGM.图6 无限大开孔平板模型Fig.6 Model of the infinite plate with a circle hole图7 无限大开孔平板问题的位移误差比较Fig.7 Displacement error comparison for the problem of the infinite plate with a circle hole图8 无限大开孔平板问题的能量误差比较Fig.8 Energy error comparison for the problem of the infinite plate with a circle hole图9 开孔平板问题的CPU时间比较Fig.9 CPU time comparison for the problem of the infinite plate with a circle hole4.3 含中心裂纹的无限大板考虑一无限大板,中心包含一长度为2a的直裂纹,在远端受单向拉应力σ作用,见图10.计算过程中,计算区域ABCD取为10 mm×10 mm,a=100 mm;E=300 MPa,v=0.3,σ =1 MPa;在计算区域内均匀分布20×20个节点.分别采用EFGM,FSMM,SSMM和NSSM求解该问题.图11为4种方法计算得到的裂纹前端应力与精确解[31]的比较.可以看出,FSMM,SSMM和NSSM的应力曲线呈现依次逼近解析解的特征,说明随着光滑区域数量的增加,光滑Galerkin无网格法的精度也随之提高;NSMM和EFGM的计算结果与解析解吻合得非常好.图10 含中心裂纹的无限大板的几何结构和载荷Fig.10 Geometry and loads of infinite plate with a center crack图11 裂纹前端(r>0,θ=0)应力比较Fig.11 Stresses comparison ahead of the crack tip(r>0,θ=0)for the near-tip crack4.4 双连拱隧道设有一双连拱形隧道,基本结构如图12所示.分别采用本文方法(NSMM)和EFGM,按照平面应变假设分析隧道围岩和衬砌在自重作用下的变形和应力分布.计算过程中,围岩及混凝土力学参数见表1;三角形背景网格划分见图13.将单元角点作为节点,共5145个背景单元,2729个节点.根据对称性,模型左侧边界上各点处x方向固定,y 方向自由,而在下边界上各点x方向自由,y方向固定.图12 双连拱隧道结构Fig.12 The structure sketch of twin-arched tunnel表1 材料参数性能Table 1 Parameters of material propertyMaterial typeE/GPa v Bulk density/(kN ·m−3)III type surrounding rock 2 0.25 22 surrounding rock of reinforced area 2.6 0.2 23 lining(C25) 28.5 0.2 25 shotcrete(C25) 28.5 0.2 23 backfilling concrete(C10) 18.5 0.2 22图13 三角形背景网格和边界条件Fig.13 Triangular background mesh and boundary conditions图14 第一主应力图Fig.14 The first principal stress图15 第二主应力Fig.15 The second principal stress图14和图15分别给出了由NSMM的计算结果绘制的第一、第二主应力云图.由图14可以看出,拱圈顶部附近、仰拱拱底附近、边墙和中墙基础底部都出现了较大的拉应力.从图15可以看出,几个形状突变较为明显的区域,由于应力集中导致了较大压应力的出现.图14和图15的结论与文献[32]的结论是一致的.NSMM与EFGM的计算结果十分接近,但EFGM的CPU用时为197.48 s,而本文方法仅为35.44 s,进一步说明本文方法具有极高的计算效率.5 结论EFGM是一种十分重要的无网格数值方法,但其计算效率低的问题长期困扰着研究工作者.本文基于广义梯度光滑技术,提出了光滑无网格Galerkin法,试图解决这一难题.本文方法的基本思想是利用MLS近似位移场函数,在两层嵌套光滑子域上计算最优光滑应变,基于GS-Galerkin弱形式推导了系统方程.两层嵌套网格由三角形背景积分网格以及连接三角形网格的3条边中点组成的4个三角形子网格构成.利用Richardson外推法,得到两层网格下对应的最优光滑应变.几个数值算例验证了本文方法的精度和计算效率.从这些数值结果可以看出,本文方法具有以下优势: (1)计算效率高.本文方法充分利用了光滑梯度技术,将区域积分转化为区域边界积分,完全避免了形成整体刚度矩阵时复杂形函数导数计算.因此,极大地提高了的计算效率.数值结果表明,在节点分布相同的情况下,NSMM比EFGM花费的计算时间要少的多.(2)计算精度高.FSMM类似SCNI,仅有线性精度:SSMM将FSMM的光滑区域细分为4个子区域,相当于引入了更多的积分取样点,因此SSMM比FSMM的精度更高:Richardson外推法又从理论上保证NSMM比FSMM和SSMM二者的精度更高.4个数值算例证明了NSMM具有几乎和EFGM相同的计算精度.(3)易于施加边界条件.本文方法在实现过程中,域内计算点采用MLS近似,而在域边界上的计算点则采用线性插值.一方面,本文方法可以自然满足线性分片试验;另一方面,边界条件的施加过程就和FEM一样了.这样,无网格边界条件施加困难的问题也就迎刃而解了.综上所述,本文提出的光滑Galerkin无网格法,既保持了无网格法的精度高的优势,又极大地减少了计算用时,同时也解决了施加边界条件的难题,具有十分广阔的发展空间.将本文方法进一步推广到结构动力学分析、大变形分析以及裂纹动态扩展等问题将是未来的研究方向.参考文献【相关文献】1 Belytschko T,Lu YY,Gu L.Element-free Galerkin methods.International Journal for Numerical Methods in Engineering,1994,37:229-2562 Belytschko T,Kronggauz Y,Organ D,et al.Meshless methods:An overview and recent putational Methods in Applied Mechanics and Engineering,1996,139:3-473 Liu WK,Jun S,Zhang YF.Reproducing kernel particle methods.International Journal for Numerical Methods in Fluids,1995,20:1081-11064 Atluri SN,Shen SP.The Meshless Local Petrov-Galerkin(MLPG)Method.Tech Science,20025 Nguyen VP,Rabczuk S,Bordas S,et al.Meshless Methods:A review and computer implementation aspects.Mathematics and Computers in Simulation,2008,79:763-8136 Rabczuk T,Bordas S,Zi G.A three-dimensional meshfree method for continuous multiplecrack initiation,nucleation and propagation in statics andputational Mechanics,2007,40:473-4957 马文涛,许艳,马海龙.修正的内部基扩充无网格法求解多裂纹应力强度因子.工程力学,2015,32(10):18-24(Ma Wentao,Xu Yan,Ma Hailong.Solving stress intensity factors of multiple cracks by using a modified intrinsic basis enriched meshless method.Engineering Mechanics,2015,32(10):18-24(in Chinese))8 Atluri SN,Zhu T.A new meshless local Petrov-Galerkin(MLPG)approach in computational putational Mechanics,1998,22:117-1279 Wang JG,Liu GR.A point interpolation meshless method based on radial basis function.International Journal for Numerical Methods in Engineering,2002,54:1623-1648 10 杨建军,郑健龙.无网格局部强弱法求解不规则域问题.力学学报,2017,49(3):659-666(Yang Jianjun,Zheng Jianlong.Meshless local strong-weak(MLSW)method for irregular domain problems.Journal of Theoretical and Applied Mechanics,2017,49(3):659-666(in Chinese)) 11 Duan QL,Li XK,Zhang HW,et al.Second-order accurate derivatives and integration schemes for meshfree methods.International Journal for Numerical Methods in Engineering,2012,92:399-42412 Beissl S,Belytschko T,Nodal integration of the element-free Galerkin puter Method in Applied Mechanics and Engineering,1996,139(19):49-6413 Chen JS,Wu CT,Yoon S,et al.A stabilized conforming nodal integration for Galerkin meshfree methods.International Journal for Numerical Methods inEngineering,2001,50:435-46614 Chen JS,Yoon J,Wu CT.Nonlinear version of stabilized conforming nodal integration for Galerkin meshfree methods.International Journal for Numerical Methods in Engineering,2002,53:2587-261515 Liu GR.A generalized gradient smoothing technique and smoothed bilinear form for Galerkin formulation of a wide class of computational methods.International Journal of Computational Methods,2008,5(2):199-23616 Liu GR,Zhang GY.A normed G space and weaked weak(W)formulation of a cell-based smoothed point interpolation method.International Journal of Computational Methods,2009,6(1):147-17917 Liu GR,Nguyen TT,Smoothed Finite Element Method.Boca Raton:CRC Press,201018 Liu GR,Zhang GY.Smoothed Point Interpolation Methods:G Space Theory and Weakened Weak Forms.New Jersey:World Scientific,201319 Bordas PAS,Rabczuk T,Hung NX,et al.Strain smoothing in FEM and puters and Structures,2010,88:1419-144320 Chen L,Rabczuk T,Bordas PAS,et al.Extended finite element method with edge-based strain smoothing(ESm-XFEM)for linear elastic crack putational Methods in Applied Mechanics and Engineering,2012,209-212:250-26521 Liu GR,Smoothed finite element methods(S-FEM):An overview and recent developments.Arch Computat.Methods Eng,2018,25:397-43522 Cui XY,Li GY.Smoothed Galerkin methods using cell-wise strain smoothing technique.Engineering Analysis with Boundary Elements,2012,36:825-83523 杜超凡,章定国.光滑节点插值法:计算固有频率下界值的新方法.力学学报,2015,47(5):839-847(Du Chaofan,Zhang Dingguo.Node-based smoothed point interpolation method:A new method for computing lower bound of natural frequency.Chinese Journal of Theoretical and Applied Mechanics,2015,47(5):839-84724 Liu F,Yu CY,Yang YT.An edge-based smoothed numerical manifold method and its application to static,free and forced vibration analyses.Engineering Analysis with Boundary Elements,2018,86:19-3025 Duan QL,Wang BB,Gao X,et al.Quadratically consistent nodal integration for second order meshfree Galerkin putational Mechanics,2014,54:353-36826 Wang BB,Duan QL,Shao YL,et al.An efficient nodal integration with quadratic exactness for three-dimensional meshfree Galerkin methods.Engineering Analysis with Boundary Elements,2016,70:99-11327 邵玉龙,段庆林,李锡夔等.功能梯度材料的二阶一致无网格.工程力学,2017,34(3):15-21(Shao Yulong,Duan Qinglin,Li Xi kui,et al.Quadratically consistent meshfree method for functionally graded materials.Engineering Mechanics,201734(3):15-21(in Chinese))28 Wang DD,Wu JC.An efficient nesting sub-domain smoothing integration algorithm with quadratically exactness for Galerkin meshfree puter Methods in Applied Mechanics and Engineering,2016,298:485-51929 Brezinski C,Zaglia MR.Extrapolation Methods:Theory and Practice.North-Holland,199130 Timoshenko SP,Goodier JN.Theory of Elasticity,3rd edition.New York:McGraw-Hill,198731 Williams ML.On the stress distributions at the base of a stationary crack.Journal of Applied Mechanics,1957,24:109-11332 徐荣桥.结构分析的有限元法与MATLAB程序设计.北京:人民交通出版社,2006(Xu Rongqiao,Finite Element Method and Matlab Programing Design for Structural Analysis.Beijing:China Communications Press,2006(in Chinese))。
无网格伽辽金法_EFGM_模拟不连续面
u ′= {u ′ 1 v′ 1 u′ 2 v′ 2 u′ 3 v′ 3 u′ 4 v′ 4} u = {u 1 v 1 u 2 v 2 u 3 v 3 u 4 v 4 } u′ = Lu
2 在 EFGM 中引入节理单元
211 节理单元应变能
L 为坐标转换矩阵 r L =
0
r
0 0
r
0 0 0
r
考虑一长度为 l 的节理单元[ 10 ] , 节理切向与 x 轴夹角为 Η , 局部坐标系 x ′ y′ 与整体坐标系 x y 如图 1 所示。 节理单元的 4 个节点是 1, 2, 3, 4, 节点位移 在局部坐标系与整体坐标中分别是 u ′ , v′ 和 u , v。
100429665 2000 08 ( 03) 20364205 J ou rnal of E ng ineering Geolog y 工程地质学报
无网格伽辽金法 (EFGM ) 模拟不连续面
庞作会① 葛修润② 王水林②
Ξ
( ①北方交通大学土木建筑工程学院 100044) ; ( ②中国科学院武汉岩土力学研究所 武汉 430071)
l
为方便讨论, 这里写出 EFGM 的控制方程[ 8 ] ζ = G 0 Κ
T
K G
u
3
f q
2x
l
〕
( 1)
( 13)
式中各符号意义见 [ 8 ]。
3 EFGM 控制方程中的未知量是 u , 而上述应变
+ (u ′ u′ 1+ 3 2 )〔
〕 ]
2x
l
∃v ′ =
1 [ (v ′ v′ 14 1 )〔 2 2x
庞作会等: 无网格伽辽金法 ( EFGM ) 模拟不连续面
紧积分算子特征值无网格Galerkin算法
而， 有
口（ ）＝ Ａ （ ） ｓ ｑ， Ｓ Ｓ Ｂ（ ）
（ ．） ２ ５
其 中，
［ （） Ａ ｓ］ ∑ ｃ（） ｓ （） 业＝ ｕｓ ｓ， （ ） ｊ
征 向 量 在 该 算 法 下 的 收敛 阶 ． 关 键 词 ： 征 值 问 题 ； 积 分 算 子 ； Ｓ 无 网格 特 紧 ＭＬ ； 中图分类号 ： ７ ． Ｏ１ ５ ３ 文献标识码 ： Ａ
１ 引 言
紧 积分算 子特征 值 问题在 数 学 、 理 以及 工程 等领 域有 着很 广泛 的应 用 ， 物 因此 其数 值计 算方法 一 直
艿 ， （ Ｚ）＝ｓｐ ｄ ｔｙ Ｚ）： ，ＹＩ ｝ ｕ ｛ｉ （ ， ｓ ＹＥ Ｉ『 ｌ ＝１ ，
那么 ， 称
（ ， ）＝ ｍａ ｛ ， ） ， ） Ｚ ｘ （ Ｚ ， （ Ｙ ｝ 为空 间 与 的 间隙 ． 定 理 ２ １ 当 充 分大 时 ， 在 与 无 关 的常数 Ｃ 使得 ． 存 ，
数． 如果 Ｐ ＝一
．
１ ） ｚ 关 算 的 投 算 ，秩 ，么Ｆ 的 点 （ 一 ｄ是 于 子 谱 影 子且 为 那 ，内 谱 由
－， ． ． 我们 记 ＝ 为这 个 特 征值 的代数 均 来自个特 征值 组成 ， 别记 为 分
值 ．
对 于 的非 零子 空 间 和 ， 令
［ ＳＪ Ｂ（ ） ＝∞ｌ ｓ （ｌ） ０ ≤ 志，≤ 志 志 Ｘ） ＾∈Ａ（ ） ｆ ） ｓｋ ，≤ （ １ ≤ （ ， Ｓ，
ｑ ＝ ｖＩ ＾ ＾．
将 （ ． ） 入 （ ．） ２５ 代 ２ ４ 得
无单元伽辽金法程序设计和算例
子程序 PHIXY BASEP BASPXY AXYINV FMABXY FORMAB
FORMAB FMABXY FMNODE FMWXY BASEP FORMW BASEP
FMNODE FMWXY BASEP
图 3 子程序 PHIXY 的结构示意图
5.1.5 其它的几个主要子程序及其功能 FOMXYZ：形成节点坐标(对于不规则形状物体则读入坐标值)。 PUTIN：输入基本数据。 EQSTR：计算等效应力σ 。 FMRCXY：划分背景网格。 GAUSDT：形成高斯积分的积分点坐标和积分权系数。 DISPST：计算节点位移和应力。 FORMDEP：计算形成弹塑性矩阵[D]ep 。 KPUT：从高斯积分点形成对整体刚度矩阵的贡献。 NODALSTRESS：计算弹塑性分析时的节点应力。 FORMK：计算完全弹性线性分析时的刚度矩阵。 FORMF：计算完全弹性线性分析时的载荷列阵。 STIF2D：计算弹塑性分析时的刚度矩阵[K ({σ }i−1 )] 和前一增量加载步末与结构中的应
来计算；
(7)求解相应的基本方程(4.24)求得位移增量，进而计算应变增量、应力增量，并把位移、 应变和应力增量叠加到加载前的水平上去；
ቤተ መጻሕፍቲ ባይዱ
(8)如还未加载到全部载荷(即加载次数< n ),则回复到步骤(5)，继续加载；否则计算停止； (9)输出有关信息。
1.2 程序的主要变量及数组说明 WIDX： x 方向的长度。 WIDY： y 方向的长度。 E： 弹性模量。 ET：弹塑性切线模量。 HL：材料塑性模量，又称为硬化参数。 V：泊松比 SSS：屈服极限。 NITR: 循环次数。 NGCELL：背景网格内高斯积分点的个数。 IPLANE：问题类型变量。当 IPLANE = 0 时为平面应变问题；当 IPLANE = 1时为平面 应力问题。
无网格Galerkin法在电磁场计算中的应用研究
无网格Galerkin法在电磁场计算中的应用研究∗曹素，龚曙光，刘翔，刘新湘潭大学机械工程学院，湖南湘潭（411105）E-mail：gongsg@摘要：有限元法是偏微分方程数值计算的强大工具，但它以网格单元为基础，存在着某些不足。
无网格法作为一种新兴的数值方法，解除了节点的网格束缚，能够消除由于网格存在所带来的缺陷。
本文以电磁场数值计算的泊松方程边值问题为研究对象，建立了无网格Galerkin法求解的离散方程，编写了MatLab程序，完成了2个电磁场问题的数值计算，所得结果与有限元法计算结果进行了比较，显示无网格Galerkin在电磁场计算中具有更好的数值精度和稳定性。
关键词：无网格Galerkin法，电磁场，有限元 MATLAB中图分类号：TM151. 引言经过近半个世纪的发展和完善，有限元(Finite Element—FE)算法已成为一个成熟而强大的计算工具，特别是在电磁场的数值模拟方面也取得了许多可用的成果[1-2]。
然而有限元法在电磁场的数值计算中并不是不存在缺点，由于有限元法对单元网格必须要满足一定的的形状要求，这使得一些特别应用的地方，有限元法就遇到了困难，如反求形状优化、移动导线以及裂纹等计算中出现的几何大变形，这时需要在计算中不断产生单元重构才能使计算顺利进行，另一方面，在微小气隙及场中有极薄铁板的时候，限于计算机的容量无法形成合理的有限单元，因此这两个方面已限制了有限元法在电磁场数值计算的中应用[3-5]。
目前一种新兴的不需要网格的数值方法即无网格方法(Meshless Method)已经出现，它为解决上述问题提供了新的希望[6]。
无网格法起源于20世纪70年代，但直到近几年移动最小二法(Moving Least-Squares—MLS)近似的引入，才使该方法在工程界得到广泛关注。
在无网格法方法中，由于形函数的构造方式不同，出现了近20多种无网格方法，如光滑质点流体动力学法、再生核质点法、无网格Galerkin(Element-free Galerkin—EFG)法、单位分解法等，其中EFG法与其它方法相比具有数值稳定、后处理方便、精度高、收敛快等特点，并认为是工程应用中发展良好和最具有应用价值的无网格方法之一[7-9]。
计算电磁学3-有限元法、里兹法、伽辽金法、矩量法
电磁波方程
Yee格式及蛙跳机制
电磁波方程的离散
激励源
Mur吸收边界条件
解的数值稳定性
Yee格式及蛙跳机制
n d 2 l E dl = 0 dt A H dS 1 = 0 H n1 dS H n dS A A t d H d l = E dA J dA 0 l A dt A
t H x 0
E
n 1 z i , j , k 1/2
Hx z
n 1 2 i , j 1/2, k 1/2
Hz
n 1 2 i 1/2, j 1/2, k
Hz x
n 1 2 i 1/2, j 1/2, k
n 1 2 J Source _y
f x x
xi
1 2 f x x f x x O x i i 2x
离散
计算机处理
1.积分 f xi x
温度对薄壁筒结构影响的伽辽金无网格法分析
￡ ￡ ， ￡ ， ， ， ］ ＝［ ￡ ，： ＝Ｌ（ “ ）
域 内任 一 点 ｚ处 的应 力 为 ：
引起 的节点位移 ， 就可以求得热应力 ｝ 。 ２ Ｊ 计算应力时应包括初应变项 ：
ｄ ＝Ｄ（ 一￡） ￡ ０ （） ２
ｄ：［ ，ｙ ｄ ， ， ， ］ ａ ，： ｒ ｒ ：Ｄｅ Ｌ（ ） ：Ｄ ３ “ ２ 则 系统势能 的泛 函表达式 为 ：
］一 ｄ Ｏ
（） ８
其中 ， ａ为材 料 的热膨 胀 系数 ， ／； ０为结 构 的初始 温 度 ｌ ｔ；
对 于线 弹性 问题 ， 变分原理为平衡位移使 系统 总势能取驻 其 值， ： 即
１２ 从 变分原 理推 导基 本方程 ．
设弹性体域为 ｎ， 力边界 为 ｓ ， 位移 边界 为 ， 则平衡 方程
）－Ｍ ）－ ）・ 训 （ ｄ ５）
｛ ＝［ ｛ Ｋ］ Ｔ｝ Ｑ｝ Ｃ］ Ｔ｝ ｛
（） １
其 中 ，Ｋ］ ［ 为传导矩 阵 ， 包含导热 系数 、 流 系数及辐射 率和 对
形状 系数 ； Ｃ］ ［ 为比热矩 阵 ， 考虑 系统 内能的增加 ； Ｔ｝ ｛ 为节 点温 度向量 ；Ｔ｝ ｛ 为温度对时间 的导数 ； Ｑ｝ ｛ 为节点 热流率 向量 ， 包含
热生 成 。
∑
ｒ （） ０ ｚ
： ２
Ｎ（ ＝［ （２Ｎ２ …Ｎ ） ｘ） Ｎ１．） （ ３ ） （ ］
（ｅ ５）
０ ＿ １
物体 由于温度变化而引起的应力成 为热应力或 温度 应力… ，
当弹性体 的温度场 由瞬态传热分析求得 时 ， 可以进一步求 出弹 就 性体各部分 的热应 力。物体 由于热膨 胀产 生变形 而引 起 的应 变
求解偏微分方程三种数值方法
数值模拟偏微分方程的三种方法介绍(有限差分方法、有限元方法、有限体积方法)I.三者简介有限差分方法(Finite Difference Methods)是数值模拟偏微分方程最早采用的方法，至今仍被广泛使用O该方法包括区域剖分和差商代替导数两个步骤。
首先将求解区域划分为差分网格，用有限个网格节点代替连续的求解区域。
其次，利用Taylor级数展开等方法将偏微分方程中的导数项在网格节点上用函数值的差商代替进行离散，从而建立以网格节点上的值为未知量的代数方程组。
该方法是一种直接将微分问题变为代数问题的近似数值解法，数学概念直观，表达简单，是发展较早且十分成熟的数值方法。
差商代替导数后的格式称为有限差分格式，从格式的精度来考虑，有一阶格式、二阶格式和高阶格式。
从差分的空间离散形式来考虑，有中心格式和迎风格式。
对于瞬态方程，考虑时间方向的离散，有显格式、隐格式、交替显隐格式等。
目前常见的差分格式，主要是以上几种格式的组合，不同的组合构成不同的差分格式。
差分方法主要适用于结构网格，网格的大小一般根据问题模型和Courant稳定条件来决定。
有限元方法(Finite Element Methods)的基础是虚位移原理和分片多项式插值。
该方法的构造过程包括以下三个步骤。
首先，利用虚位移原理得到偏微分方程的弱形式，将计算区域划分为有限个互不重叠的单元（三角形、四边形、四面体、六面体等），在每个单元上选择合适的节点作为求解函数的插值点，将偏微分方程中的变量改写成由各变量或其导数的节点值与所选用的分片插值基函数组成的线性表达式，得到微分方程的离散形式。
利用插值函数的局部支集性质及数值积分可以得到未知量的代数方程组。
有限元方法有较完善的理论基础，具有求解区域灵活（复杂区域）、单元类型灵活（适于结构网格和非结构网格）、程序代码通用（数值模拟软件多数基于有限元方法）等特点。
有限元方法最早应用于结构力学，随着计算机的发展已经渗透到计算物理、流体力学与电磁学等各个数值模拟领域。
无网格伽辽金方法在线弹性断裂力学中的应用研究
山东大学硕士学位论文无网格伽辽金方法在线弹性断裂力学中的应用研究摘要（在处理裂纹扩展这类动态不连续性问题时，传统的计算方法如有限元法、＼有限差分法等常需要网格重构。
这样不仅增加了计算工作量，而且会使计算精度严重受损。
无网格方法中，由于采用基于点的近似，网格可以彻底或部分地、消除，因此可以完全抛开网格重构，从而保证了计算精度广本文在系统分析了前人所做工作的基础上，对无网格伽辽金方法（ＥＦＧＭ）做了部分改进，并用算例对其正确性和有效性进行了验证。
本文中主要的研究成果和结论有：／基于移动最小二乘近似的ＥＦＧＭ是目前应用最广泛的无网格方法，由于移动＼最小二乘形函数一般不具有常规有限元形函数所具有的插值特性，即ＥＦＧＭ的近＼，似函数不通过节点变量，本质边界条件的处理成为ＥＦＧＭ实旖中的一个难点可本文在处理本质边界条件时，采用了再生核质点方法中的完全变换法，实现了本质边界条件在节点处的精确施加。
权函数的使用是ＥＦＧＭ和其它无网格方法的精华所在，本文采用了一种基于ｔ一分布的新型权函数，在一定程度上提高了ＥＦＧＭ的计算精度。
ｒ／影响半径大小的取值对最终的场函数近似解或其导数有较大影响，传统方＼｜法是在整个求解域内使用统一的影响半径。
、ｊ本文针对裂纹扩展中的实际情况，√对动态影响半径法作了进一步的补充和改进。
即在均匀分布节点区域，采用与基函数相对应的规定节点数来确定影响半径的大小；而在局部加密节点邻域，根据节点加密情况，相应地增加确定影响半径所需的节点数。
本文分别计算了单一型和复合型裂纹的应力强度因子，计算结果表明使用部分扩展基函数不仅能获得较高的计算精度，而且积分围线对它的计算结果影第１页山东大学硕士学位论文响较小，计算稳定性好。
用ＥＦＧＭ模拟了拉剪复合型裂纹的扩展行为，由于避免了有限元方法中网格重构的繁琐，大大简化了裂纹扩展的模拟工作。
，计算结果证明本文模拟的裂纹扩展轨迹与前人的研究结果符合得较好。
通过本课题的研究工作，进一步发展和完善了ＥＦＧＭ，为其在断裂力学问题以及其它结构计算问题中的应用奠定了良好的理论基础；此外，也为进一步研究复杂的断裂问题，如弹塑性材料的裂纹扩展问题、三维裂纹扩展问题、动态裂纹扩展问题以及界面断裂力学问题等做了一些有益的基础准备。
典型无网格法
无网格法典型无网格法V伽辽金型无网格法V配点型无网格法V基于局部弱形式和边界积分方程的无网格法V最小二乘无网格法V物质点法V EFG¾MLS¾Galerkin V RKPM¾RK¾Galerkin V PIM¾PIG l ki¾Galerkin等效积分弱形式（虚位移原理）MLS 近似：()()u =x N x dV 计算量大精度高稳定性V 精度高，稳定性好V 需要背景网格进行积分V 系数矩阵对称V 不易施加本质边界条件处理V背景网格积分,,d I i J j N N ΩΩ∫()()()()=ΔΩ+ΔΓP N x f x N x t x 11I I I I I II I ==∑∑零能模态V 单位分解积分11.函数ψk (x )只定义在子域Ωk 上；2.子域Ωk 相互重叠，且它们完全覆盖了域Ω；3l=3.函数ψk (x )满足单位分解条件1()d ()()d k k k f f ψΩΩ∩Ω=ΩΩ∑∫∫x x x配点型无网格法V FPM¾MLS¾CollocationV SPH¾KA¾Collocationp es ess c ouds V Hp meshless clouds ¾PUC ll ti¾Collocation基于局部弱形式和边界积分方程的无网格法V MLPG¾MLS¾LPGV LBIE¾MLS¾LBIENV BNM¾MLS¾BIEgn pn 1iI IJ iJJ p m v ==∑1IJ p Ip Jpp m m N N ==∑质量阵求逆iIp iJv?质量阵求逆！gpn n 11I IJ p IpJ p m m m N ====∑∑对角质量阵iI I iIp m v =V已知t k 时刻的物理量，求t k +1时刻的物理量11.更新网格结点数据kk Ip Ipm m N=∑ppn k k kiIip Ip p m v N=1p ppp =∑int,ext,kk k iI iIiIf ff=−2.在背景网格结点上积分动量方程并施加固定边界条件1k k k pp f t+=+ΔiIiI iI10,0k k iIiIp f +==在固定边界上6.更新密度，应力k k k 1/(1)p p iipρρε+=+Δ1k k kkij ij ij ij+=ijpijp ijp ijp r σσσ+Δ+Δkk k kk r σσΔ=ΔΩ−ΔΩ其中ijpijpijpijpijpk ijpσΔ根据弹塑性本构关系更新7.进行下一个时间步循环。
无网格——精选推荐
第一章绪论计算流体力学的发展现状计算流体力学(Computational Fluid Dynamics)是现代流体力学中的一个重要学科分支。
作为一门多学科交叉融合而形成的新兴学科，它是流体力学、计算数学和计算机科学相结合的产物。
随着计算机性能的飞速提高以及数值计算方法的不断发展，计算流体力学技术正在逐渐走向成熟。
计算流体力学经历了数值求解拉普拉斯方程、小扰动速势方程、全速势方程、Euler方程和Navier-Stokes方程等发展阶段。
20世纪80年代以前，由于受到计算机技术的限制，计算流体力学的数值模拟主要以求解拉普拉斯方程、小扰动速势方程、全速势方程为主，其中有代表性的是基于拉普拉斯方程的面源法以及有限差分法求解小扰动速势方程和全速势方程。
在随后的二十多年中，在计算机技术发展的推动和广大计算流体力学工作者的努力下，计算流体力学在求解Euler方程和Navier-Stokes方程以及数值模拟复杂流场方面都取得了重大突破。
在此期间，计算流体力学数值模拟的方法以有限差分法、有限体积法、有限元法为主。
随着诸如TVD格式、ENO格式、NND格式等高阶精度、高分辨率差分格式的提出，计算流体力学对激波、漩涡等复杂问题的模拟能力也有了很大的提高。
目前，计算流体力学工作者正致力于研究和发展更高精度(二阶以上)的计算格式和方法，以适应更精细、更复杂的流动研究和设计的需要。
计算流体力学研究的一个重要分支是计算网格的生成技术，它是计算流体力学走向工程实用阶段所必须面临的关键技术之一。
一般来讲，适合工程使用的网格生成技术应该具备以下特点:(1) 网格生成过程直观明了、简单易行、效率高、自动化程度好。
(2) 通用性、普适性好，对复杂外形、复杂流动的适应能力强。
(3) 网格几何灵活性好，尺度变化易于控制，网格自适应加密简便易行。
目前，已经成熟并走向工程实用中的计算网格有结构网格、非结构网格以及结构非结构的混合网格。
在结构网格方面，出现了代数生成网格法、解微分方程生成网格法、保角变换法等多种网格生成方法，网格类型也由单一的C型网格、0型网格、H型网格发展到嵌套网格和多块对接网格等。
10_无网格方法
a( x) A1 ( x) B( x)u (8)
把式(8)代入式(1)得近似函数为
u ( x, x) pT ( x) A1 ( x) B( x)u N T ( x, x)u (9)
U T RA T a d V T RB T a d 0
式中 RA a A a 和 RB a B a 为余量。 虽然 U 和V
T T T Tຫໍສະໝຸດ
为任意权函数，但在实际应用时，不可能也不需要取无穷多个权函数。 与试函数表达式类似，可以把权函数也写成已知基函数的组合，即
pT ( x) [1, x, y, x 2 , xy, y 2 ]
m6
m 10
基函数的个数 m 、 基函数中包含的完备多项式的最高阶数 n 和问题的维数 nDim 之间的关系为
m
(n 1)(n 2) (n nDim ) nDim !
(3)
在移动最小二乘近似的（MLS）中，坐标 ai ( x ) 的选取应该使近似位移
p ( x) [1, x, y, r cos , r sin , r sin sin , r cos sin ] 2 2 2 2
T
式中：r 为某点距裂纹尖端的距离， 为该点与裂纹尖端连线与裂纹线的夹角。
若把式(1)作为有限单元容许位移函数，则 p1 ( x) 1 可以保证单元容许位
求解方法。只要试函数是利用离散点来建立的，则由紧支试函
数加权余量法导出的各种近似方法都称为无网格方法。
紧支近似函数 紧支近似函数是定义在局部区域(支撑区域)中的函数，它只 在支撑域中有定义，而在支撑域外为零。在二维问题中，支撑 域一般为圆形或矩形，与求解域相比，支撑域是一个很小的区 域，并且是可以互相重叠。有限元网格表示的区域是不能彼此 重叠的。
有限元-伽辽金法
单元结点温度列阵
e
12
3
Ni
(x,
y)
ai
bi x 2A
ci
y
i、j、m
e [Ni
Nj
N e
Nm
] ij
N
e e
Ni
Nj
Nm
m
9.3二维稳态热传导有限元方程
二、三结点三角形单元
Ni
(x,
y)
ai
bi x 2A
ci
y
i、j、m
Nr (x, y) ( ar br x cr y ) br
3
NeTQd NeTqds+ NeTh ds
e
2
3
9.3二维稳态热传导有限元方程
一、有限元方程
1
2
x 2 e
y22NxeT
Ne
Qx=0
NeT Qd
Ne y
T
Nye ,xdnx
n hNeT
3 ,y
Ne
y
ds
q
NeTqds+,xnNxeThd,ysny h
x (a ,b)
x=a
x=b
b
b
b
x,1 x,2 x (x)Du(x)dx+1(x)B1 udx+2 (x)B2 udx 0
a
a
a
9.1 伽辽金方法
一、加权余量法
b
b
b
x,1 x,2 x (x)Du(x)dx+1(x)B1 udx+2 (x)B2 udx 0
a
a
a
b
b
b
x,1 x,2 x (x)Du%(x)dx+1(x)B1 u%dx+2 (x)B2 u%dx 0
断裂力学裂纹扩展
断裂力学裂纹扩展做裂纹扩展仿真确实比较难，目前一般都是以弹性断裂力学为基础，二维裂纹扩展容易一些，三维裂纹比较复杂，如果仅是要获得扩展寿命，裂纹长度，可以自己编程做，我是这样做的。
如果要想获得不同裂纹前沿的应力应变场和K，模拟结构裂纹随载荷的动态真实变化，可能要借助软件：（1) Beasy，边界元软件，将三维问题解化为二维问题，比较方便。
（2) Fatigue软件，也还可以，但对复杂结构很难胜任。
（3) FE-fatigue 也不错（4) FRANC3D。
至于计算，常用的方法有：（1）Prescribed Method特点：裂纹只能沿单元边界扩展。
（2）Analytical Geometry Method特点：将几何和载荷、约束分解为简单的解析形式。
（3）Known Solution Method特点：查表求已知解。
两个重要软件：NASGRO and AFGROW（4）Meshfree method美国西北大学做的最好。
优点是不需重新划分网格。
（5）Adaptive BEM/FEM自适应网格边界元/有限元，用的较广。
（6）Lattice method格子方法（7）Atomic method一般使用分子动力学方法。
（8）Constitutive method在本构方程里引入破坏准则，无需预先引入裂纹。
如本人上篇帖子。
（9）Cohesive element使用cohesive element。
断裂学科研究的新趋向第十届国际断裂大会（ICF10）的情况介绍四年一届的国际断裂大会(Int. Conference of Frature, ICF-10)于2001年12月3日～12月6日在美国夏威夷召开。
与会的有来自44个国家的代表约610人。
中国参加会议的代表并有论文在论文集上发表的计34人（含中国香港10人），其中部分代表因故未能到会。
此次会议的举办是成功的，现将会议的简要情况与参加会议的体会及有关建议分别作简单汇报于下。
水下爆炸船体结构响应间断伽辽金法数值模拟
水下爆炸船体结构响应间断伽辽金法数值模拟于福临;郭君;姚熊亮;任少飞【摘要】In order to solve the underwater explosion flow field with large discontinuities, Level Set method was applied to track the interface position of the multi⁃medium flow, Ghost Fluid method was used to calculate the physical parameter of both sides of the interface, time and space were discretized by Runge⁃Kutta Discontinuous Galerkin Method, Euler equations of the flow field were solved. One⁃dimensional andtwo⁃dimensional assessments were conducted by RKDG approach. The results reflect the phenomena of underwater explosion shock wave generation, propagation, reflection and explosion products expansion. Finally, the shock responses and damage characteristics of hull plates under shock load were simulated with the nonlinear FEM softeware ABAQUS. The RKDG method can be applied to simulate the hull plates response with high accuracy. The response of hull plates is inversely proportional to the blast center distance.%为求解水下爆炸强间断流场，采用Level Set方法定位多相流界面位置，应用虚拟流体方法处理邻近界面两侧物理量，并用RKDG方法进行空间和时间的离散，求解流场的Euler方程，并进行一维、二维评价，计算结果能够较好地反映水下爆炸冲击波产生、传播、反射和爆炸产物的膨胀等现象。
变形网格上的非多项式Galerkin投影
变形网格上的非多项式伽辽金投影图一：我们的方法可以减少飞鸟周围的流体模拟，比相应的全模拟快2000倍以上。
减少在该建筑场景的辐射计算，比相应的全辐射快113倍以上。
摘要：本文将伽辽金投影扩展到图形中常见的大量非多项式位函数。
我们通过在两个截然不同的问题上的应用证明方法的广泛适用性，即流体模拟和辐射渲染都采用变形网格。
标准的伽辽金投影不能有效地近似这些现象。
我们的方法与此不同，能使这些复杂的非多项式系统紧凑表示和逼近，其中包括商数及多项式的根。
我们依靠表示各函数的模型缩减作为张量积，矩阵求逆和矩阵根的组成部分。
一旦某函数在该形式中表示出来，它就可以很容易的模型缩减，并且它的降阶形式能够随时进行评估，存储器成本只依赖于降维空间的维数。
CR种类：I.3.7[计算机图形]:三维图像和现实动画，辐射算法；I.6.8[模拟和建模]模拟的种类——动画；关键字：缩减模型，流体模拟，固流耦合，辐射。
1、简介伽辽金投影在图形加速上令人吃惊。
然而，尽管其有广泛的应用——从全局照明到流体——该方法的关键限制是底层现象必须是多项式，这种约束限制了其在计算机图形上的适用性。
本文提出伽辽金投影在组成任何初等代数运算函数的有效延伸——在算术加有理根的四则运算——从而在图形中贯穿这种模式缩减方法的适用性。
为了证明我们的方法的广泛适用性，我们将其用于两个显著不同的问题上：辐射渲染和流体模拟。
尽管这两种现象都可以通过固定网格法用多项式格式表示，但是我们认为实现几何畸变需要非多项式计算，以表现动力和外观的改变。
理论上可以在这些函数中应用标准的伽辽金投影，但是这不会提高任何的运行速度。
我们的技术与此不同，能够高效地模拟这些复杂的非多项式系统。
与标准伽辽金投影相似，我们的方法不仅保留关键的最优保证，而且保证在缩减空间方面有着紧凑和易分析的模型。
我们的方法有广泛的应用范围。
环境交互方面，如在视频游戏和建筑设计应用上逐渐地融入物理模拟。
我们的技术可以围绕包括人物角色和动物在内的操纵对象加入流动效果。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
!*******************************************************! 无网格伽辽金法程序主要计算模块简单例子! Kernel modules for element-free Galerkin Method! By Lu Xinzheng of Tsinghua University!*******************************************************module TypeDefuse Lxz_Toolsimplicit none;!integer,parameter :: MLS_m=6; !定义基为二次基integer,parameter :: MLS_m=3; !定义基为1次基integer,parameter :: NInt=4; !定义高斯积分点数目integer,parameter :: NDOF=2; !定义自由度数目integer,parameter :: MLS_n=6;type :: typ_Kcol !定义变带宽总体刚度矩阵real(rkind),pointer :: row(:)end type typ_Kcoltype :: typ_GValue !定义总体控制变量integer(ikind) :: NNode; !总节点数integer(ikind) :: NCell; !总积分区域数integer(ikind) :: NIntPoint; !总高斯积分点数integer(ikind) :: NPress; !总压力数integer(ikind) :: NSupport; !总位移约束数real(rkind) :: dm; !节点影响半径real(rkind) :: c; !权函数相对权重参数real(rkind) :: Plenty !罚函数大小end type typ_GValue;type :: typ_Node; !定义节点类型real(rkind) :: X,Y; !节点坐标real(rkind) :: dX, dY; !节点位移end type typ_Node;type :: typ_IntPoint; !定义高斯积分点类型!以下为高斯点几何数据real(rkind) :: CenterX, CenterY; ! 中心点坐标real(rkind) :: Len, Wide; ! 高斯积分点的积分面积integer(ikind) :: NodeNo(MLS_n); ! 包含节点的编号!以下为移动最小二乘所需数据real(rkind) :: a,b !x,y的最大距离real(rkind) :: dm; !节点影响半径real(rkind) :: c; !权函数相对权重参数real(rkind) :: Pasi(MLS_n); !形函数矩阵, 维数是 Pasi(MLS_n)real(rkind) :: Pasidx(MLS_n); !形函数对x求导, 维数是 Pasidx(MLS_n) real(rkind) :: Pasidy(MLS_n); !形函数对y求导, 维数是 Pasidy(MLS_n) !以下为高斯点力学数据real(rkind) :: Strain(3); !高斯点的应变 Strain=B*uireal(rkind) :: D(3,3); !高斯点材料性质矩阵real(rkind) :: Stress(3); !高斯点的应力 Stress=D*Strainreal(rkind) :: E,mu;end type typ_IntPoint;type :: typ_Load !定义节点荷载integer(ikind) :: NodeNo !高斯点编号real(rkind) :: Value(2) !荷载大小end type typ_Loadtype :: typ_Support !定义位移约束integer(ikind) :: NodeNo !约束节点号real(rkind) :: Value(2) !位移大小end type typ_Support! containsend module TypeDefmodule MLSuse Lxz_Tools;use TypeDef;implicit none;contains!SetIntPoint(IntPoint,Node) !设定积分区域内各个参数的数值subroutine SetIntPoint(IntPoint,Node) !设定积分区域内各个参数的数值type(typ_IntPoint) :: IntPoint;type(typ_Node) :: Node(:);integer(ikind) :: I,J; !循环变量real(rkind) :: d,x1,y1,dx,dy; ! d节点和积分点之间的距离, x1,y1节点和积分点之间x,y方向的距离real(rkind) :: temp1; ! 中间变量real(rkind) :: tempAdx(MLS_m,MLS_m),tempAdy(MLS_m,MLS_m) !对A求导时的中间变量 real(rkind) :: Ptemp1(1,MLS_m),Ptemp2(MLS_m,1); !为了向量相乘设定的二维数组中间变量real(rkind) :: P(MLS_m,MLS_n); !基函数数值, 维数应该是P(MLS_m,MLS_n)real(rkind) :: Px(MLS_m);real(rkind) :: Pdx(MLS_m); !基函数数值对x求导, 维数应该是Pdx(MLS_m)=(0,1,0)real(rkind) :: Pdy(MLS_m); !基函数数值对y求导, 维数应该是Pdy(MLS_m)=(0,0,1)real(rkind) :: A(MLS_m,MLS_m); !矩阵Areal(rkind) :: Adx(MLS_m,MLS_m); !矩阵A对x求导real(rkind) :: Ady(MLS_m,MLS_m); !矩阵A对y求导real(rkind) :: InvA(MLS_m, MLS_m); !A的逆矩阵real(rkind) :: B(MLS_m,MLS_n); !矩阵B, 维数应该是B(MLS_m,MLS_n)real(rkind) :: Bdx(MLS_m,MLS_n); !矩阵B对x求导, 维数应该是Bdx(MLS_m,MLS_n)real(rkind) :: Bdy(MLS_m,MLS_n); !矩阵B对y求导, 维数应该是Bdy(MLS_m,MLS_n)real(rkind) :: Weight(MLS_n); !权函数, 维数是Weight(MLS_n)real(rkind) :: Weightdx(MLS_n); !权函数对x求导, 维数是Weightdx(MLS_n)real(rkind) :: Weightdy(MLS_n); !权函数对y求导, 维数是Weightdy(MLS_n)P=0.0d0;B=0.0d0Bdx=0.0d0;Bdy=0.0d0;Weight=0.0d0;Weightdx=0.0d0;Weightdy=0.0d0;do I=1, MLS_nx1=Node(IntPoint%NodeNo(I))%X;y1=Node(IntPoint%NodeNo(I))%Y;! 求节点和高斯积分点之间的距离d=sqrt( (x1-IntPoint%CenterX)**2 + (y1-IntPoint%CenterY)**2);!d=d/IntPoint%dm;! 防止有点超出积分区if(d<IntPoint%dm) then! 设定积分点P的值, 二次基, 参见李卧东论文公式(2.2b)Px=(/1.0d0,IntPoint%CenterX,IntPoint%CenterY/);!设定积分点关于x,y求导时的值, 二次基,Pdx=(/0.0d0,1.0d0,0.0d0/);Pdy=(/0.0d0,0.0d0,1.0d0/);! 设定节点二次基的值P(:,I)=(/1.0d0,x1,y1/);! 计算各个节点的权重, 参见李卧东论文公式(2.12)Weight(I)=(exp(-(d/IntPoint%c)**IntPoint%k)-exp(-(IntPoint%dm/IntPoint% c)**IntPoint%k))&/(1.0-exp(-(IntPoint%dm/IntPoint%c)**IntPoint%k));! 计算各个节点权重的导数temp1=IntPoint%CenterX**2-2.0d0*IntPoint%CenterX*x1+x1**2+IntPoint%Cent erY**2-2.0d0*y1*IntPoint%CenterY+y1**2;Weightdx(I)=-0.5d0*(temp1**0.5d0/IntPoint%c)**IntPoint%k*IntPoint%k*& (2.0d0*IntPoint%CenterX-2*x1)*exp(-(temp1**0.5d0/IntPoint%c)**IntPo int%k);Weightdx(I)=((Weightdx(I)/temp1))&/(1.0d0-exp(-(IntPoint%dm/IntPoint%c)**IntPoint%k));Weightdy(I)=-0.5d0*(temp1**0.5d0/IntPoint%c)**IntPoint%k*IntPoint%k*& (2.0d0*IntPoint%CenterY-2*y1)*exp(-(temp1**0.5d0/IntPoint%c)**IntPo int%k);Weightdy(I)=((Weightdy(I)/temp1))&/(1.0d0-exp(-(IntPoint%dm/IntPoint%c)**IntPoint%k));elseWeight(I)=0; Weightdx(I)=0; Weightdy(I)=0end ifend doA=0;Adx=0; Ady=0;do I=1, MLS_n! 得到P'do J=1,MLS_m;Ptemp1(1,J)=P(J,I);end doPtemp2(:,1)=P(:,I);A=A+Weight(I)*matmul(Ptemp2,Ptemp1);B(:,I)=Weight(I)*P(:,I);Adx=Adx+Weightdx(I)*matmul(Ptemp2,Ptemp1);Ady=Ady+Weightdy(I)*matmul(Ptemp2,Ptemp1);Bdx(:,I)=Weightdx(I)*P(:,I);Bdy(:,I)=Weightdy(I)*P(:,I);end do;InvA=matinv(A);TempAdx=-matmul(InvA,matmul(Adx,InvA));TempAdy=-matmul(InvA,matmul(Ady,InvA));do I=1, MLS_n;IntPoint%Pasi(I)=0;IntPoint%Pasidx(I)=0; IntPoint%Pasidy(I)=0;IntPoint%Pasi(I)=dot_product(Px,matmul(InvA,B(:,I)));IntPoint%Pasidx(I)=dot_product(Pdx,matmul(InvA,B(:,I)))+&dot_product(Px,matmul(TempAdx,B(:,I))+matmul(InvA,Bdx(:,I)));IntPoint%Pasidy(I)=dot_product(Pdy,matmul(InvA,B(:,I)))+&dot_product(Px,matmul(TempAdy,B(:,I))+matmul(InvA,Bdy(:,I)));end do;returnend subroutine SetIntPointend module MLSmodule M_Int !高斯积分区力学数据计算use Lxz_Tools;use TypeDef;use MLS;implicit none;contains! GetD(IntPoint) ! 计算材料性质矩阵D, 注意: 可以根据材料性质加以改变! GetM_K(IntPoint) ! 计算积分点的刚度矩阵subroutine GetD(IntPoint) ! 计算材料性质矩阵D, 注意: 可以根据材料性质加以改变 type(typ_IntPoint) :: IntPoint;IntPoint%D=0.0d0;!采用弹性平面应力IntPoint%D(1,:)=(/1.0d0,IntPoint%mu,0.0d0/);IntPoint%D(2,:)=(/IntPoint%mu,1.0d0,0.0d0/);IntPoint%D(3,:)=(/0.0d0,0.0d0,(0.5d0*(1.0d0-IntPoint%mu))/);IntPoint%D=IntPoint%E/(1.0d0-IntPoint%mu**2)*IntPoint%Dreturn;end subroutine GetD;function GetM_B(IntPoint) result(M_B)type(typ_IntPoint) :: IntPoint;real(rkind) :: M_B(3,2*MLS_n);Integer(ikind) :: IM_B=0.0d0;do I=1,MLS_nM_B(1,2*I-1)=IntPoint%Pasidx(I);M_B(2,2*I)=IntPoint%Pasidy(I);M_B(3,2*I)=IntPoint%Pasidx(I);M_B(3,2*I-1)=IntPoint%Pasidy(I);end doreturnend function GetM_Bfunction GetM_K(IntPoint) result(M_K)! 计算积分点的刚度矩阵type(typ_IntPoint) :: IntPointreal(rkind) :: M_B(3,2*MLS_n),M_K(2*MLS_n,2*MLS_n); !,M_B1(:,:)M_B=0.0d0; M_K=0.0d0;M_B=GetM_B(IntPoint);M_K=matmul(transpose(M_B),matmul(IntPoint%D,M_B));M_K=M_K*IntPoint%Len*IntPoint%Wide;returnend function GetM_Ksubroutine GetStress(IntPoint, Node)type(typ_IntPoint) :: IntPoint(:)type(typ_Node) :: Node(:)real(rkind) :: DispVecter(2*MLS_n)real(rkind) :: Xinteger(ikind) :: I,J;do I=1, size(IntPoint)do J=1,MLS_nX=Node(IntPoint(I)%NodeNo(J))%dX;DispVecter(J*2-1)=X;X=Node(IntPoint(I)%NodeNo(J))%dY;DispVecter(J*2)=X;end doIntPoint(I)%Strain=matmul(GetM_B(IntPoint(I)),DispVecter);IntPoint(I)%Stress=matmul(IntPoint(I)%D,matmul(GetM_B(IntPoint(I)),DispVect er));end doreturnend subroutine GetStressend module M_Int。