不连续伽辽金时域面积分方程方法的研究

二维无网格伽辽金法简单例子程序

!******************************************************* ! 无网格伽辽金法程序主要计算模块简单例子 ! Kernel modules for element-free Galerkin Method ! By Lu Xinzheng of Tsinghua University !******************************************************* module TypeDef use Lxz_Tools implicit 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_Kcol type :: 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*ui real(rkind) :: D(3,3); !高斯点材料性质矩阵 real(rkind) :: Stress(3); !高斯点的应力 Stress=D*Strain real(rkind) :: E,mu;

重构的无单元伽辽金形函数

第!"卷第#期$%%&年$月西安建筑科技大学学报’自然科学版( )*+,-./0/,1234567829:;782’<.=>6.?@7,;/7;A B,=,3/( C3?2!"<32# D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D E;F2$%%&重构的无单元伽辽金形函数 孟闻远#G$G李秀芹$G卓家寿# ’#2河海大学土木工程学院G江苏南京$#%%H"I$2华北水利水电学院土木工程系G河南郑州J K%%%"(摘要L分析了目前已有的无单元法在发展过程中存在的问题L如本质边界条件的处理M太多的求逆矩阵运算 而产生过度的计算量等问题G这些问题难以克服且互相牵涉2针对这些问题G从本质上研究了产生这些问题的内在 原因G然后在继承已有无单元法优秀内核的基础上G借助传统数学分析的工具G开发了新的形函数G从根本上解除了 这些障碍I对非凸区影响域问题G提出了新的准则2算例证明L新的形函数及准则非常有效2 关键词L无单元形函数I移动插值技术I影响域 中图分类号L N!J文献标识码L5文章编号L#%%&O P H!%’$%%&(%#O%%%K O%J Q 无单元法’R;S8?;S S R;=83B(T#U在断裂力学分析等方面显示了强劲的生命力2尤其是无单元伽辽金法’A E V(等方法G因具有坚实的数学物理基础而受到人们的青睐G但这些依W移动最小二乘’RX@(T$U Y 技术为基础的无单元法到目前为止一直被两大问题所困扰G并成为无单元法发展中的障碍L一是本质边界条件的处理问题G是由W RX@法Y形成的逼近函数无W过点插值性Y引起的I另一个是太多的计算量问题2这是因为RX@法确定逼近函数系数时要对每个计算点进行繁重不堪的求逆运算2 #加权移动最小二乘法技术内涵与继承 移动最小二乘’RX@(法用具有局部紧支撑性的权函数来形成局部逼近G再通过局部逼近在区域内’Z(的任意性回归到整体区域’Z(2W移动性Y M W局部紧支性Y M W高阶连续性Y是它的突出优点G但遗憾的是该形函数没有过点拟合性质G这对本质边界条件处理很困难2人们试图通过各种办法来解决这一问题G如无单元与有限元结合形成混合法处理本质边界条件G使得无单元法的实现变得复杂化I X.[6./[;乘子法增加了大量的未知量G也破坏了刚度阵的带状性I修正的X.[6./[;乘子法T$U G其物理量的大小难以确定I配点法显然有失精度2针对这个问题G文献T!U提出了用奇异权函数的概念G由此形成的移动最小二乘技术G具有过点插值的性质2这种处理技术实现了过点插值G但对每点的求逆运算异常繁重2为减少计算工作量G文献T!U试图对基函数进行正交化处理G但对每个\处的计算作正交化处理同样需要大量的运算2也就是说便于本质边界条件处理的过点插值技术G因依赖RX@技术仍因工作量繁重而无法完成2因此即要有W插值性Y又要避免W求逆运算Y G而且还要保持无单元法优点的新形函数迫切需要重构出来2RX@W插值Y技术虽然不易实现G但却给出了重要的暗示L局部紧支集上的奇异权函数具有W移动插值Y功能2若继承这种技术且避免求逆运算G将会减少工作量且实现W插值Y性质G这自然使我们想到了@8;].6B插值L即取^’#(_#一项作为基函数的移动最小二乘法2但仅用@8;].6B函数无法构造完备的多项式成分G仅有‘%阶连续性2为此G我们借助泰勒多项式基把多项式高阶基函数构造出来G或者说把@8;].6B函数看作%阶基函数的插值成分G而#阶以上的基函数G由泰勒展开中的多项式基来完成2由此弥补 Q收稿日期L$%%!O##O$J修改稿日期L$%%K O%P O%K 基金项目L国家自然科学基金资助项目’K%%P H%%K( 作者简介L孟闻远’#H&K O(G男G河南郾城人G教授G博士研究生G从事教学M科研及管理工作2研究方向L数值计算新方法研究M结构静M动态分析及稳定性分析等2 万方数据

三种方法概述

有限差分 有限差分方法(Finite Difference Method)是计算机数值模拟最早采用的方法，至今仍被广泛运用。该方法将求解域划分为差分网格，用有限个网格节点代替连续的求解域。有限差分法以Taylor级数展开等方法，把控制方程中的导数用网格节点上的函数值的差商代替进行离散，从而建立以网格节点上的值为未知数的代数方程组。该方法是一种直接将微分问题变为代数问题的近似数值解法，数学概念直观，表达简单，是发展较早且比较成熟的数值方法。对于有限差分格式，从格式的精度来划分，有一阶格式、二阶格式和高阶格式。从差分的空间形式来考虑，可分为中心格式和逆风格式。考虑时间因子的影响，差分格式还可以分为显格式、隐格式、显隐交替格式等。目前常见的差分格式，主要是上述几种形式的组合，不同的组合构成不同的差分格式。差分方法主要适用于有结构网格，网格的步长一般根据实际地形的情况和柯朗稳定条件来决定。构造差分的方法有多种形式，目前主要采用的是泰勒级数展开方法。其基本的差分表达式主要有三种形式：一阶向前差分、一阶向后差分、一阶中心差分和二阶中心差分等，其中前两种格式为一阶计算精度，后两种格式为二阶计算精度。通过对时间和空间这几种不同差分格式的组合，可以组合成不同的差分计算格式。 有限元法

有限元方法(Finite Element Method)的基础是变分原理和加权余量法，其基本求解思想是把计算域划分为有限个互不重叠的单元，在每个单元内，选择一些合适的节点作为求解函数的插值点，将微分方程中的变量改写成由各变量或其导数的节点值与所选用的插值函数组成的线性表达式，借助于变分原理或加权余量法，将微分方程离散求解。采用不同的权函数和插值函数形式，便构成不同的有限元方法。有限元方法最早应用于结构力学，后来随着计算机的发展慢慢用于流体力学的数值模拟。 在有限元方法中，把计算域离散剖分为有限个互不重叠且相互连接的单元，在每个单元内选择基函数，用单元基函数的线形组合来逼近单元中的真解，整个计算域上总体的基函数可以看为由每个单元基函数组成的，则整个计算域内的解可以看作是由所有单元上的近似解构成。 常见的有限元计算方法是由变分法和加权余量法发展而来的里兹法和伽辽金法、最小二乘法等。根据所采用的权函数和插值函数的不同，有限元方法也分为多种计算格式。从权函数的选择来说，有配置法、矩量法、最小二乘法和伽辽金法，从计算单元网格的形状来划分，有三角形网格、四边形网格和多边形网格，从插值函数的精度来划分，又分为线性插值函数和高次插值函数等。不同的组合同样构成不同的有限元计算格式。对于权函数，伽辽金(Galekin)法是将权函数取为逼近函数中的基函数；最小二乘法是令权函数等于余量本身，而内积的极小

伽辽金法求简支平板的屈曲载荷

This program calculates the buckling load for a simply supported plate using the Galerkin method. The region considered for the plate in this case is, x=0, x=1, and y=0, y=1 (* Definition of trial function which satisfies the given boundary conditions. The least order defined for this case.*) order=8; poly[n_]: Sum [a[(i+j) (i+j+1) / 2+i+1] x^i y^ j, {i, 0,n},{j, 0,n - i}] (*Generating the trial function using the above function definition*) poly[order] (*Applying the boundary conditions *) eq1 =CoefficientList[poly[order] / . x?>0 , y] eq2 =CoefficientList[poly[order] / . x?>1 , y] eq3 =CoefficientList[poly[order] / . y?>0 , x] eq4 =CoefficientList[poly[order] / . y?>1 , x] 63 eq5 =CoefficientList[D[poly[order] , {x, 2}] +v( D[poly[order],{y, 2}] ) / . x?>0 , y] eq6 =CoefficientList[D[poly[order] ,{x, 2}] +v( D[poly[order],{y, 2}] ) / . x?>1 , y] eq7 =CoefficientList[D[poly[order] ,{y, 2}] +v( D[poly[order],{x, 2}] ) / . y?>0 , x] eq8 =CoefficientList[D[poly[order] ,{y, 2}] +_v( D[poly[order],{x, 2}] ) / . y?>1 , x] (*Flatten command used for grouping all the eight equations together *) Eq = Flatten[ { eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8} ] (* Tabulating all the eight equations using the Table command *) Equn = Table[ Eq[ [ i ]] == 0, {i ,1, Length[Eq] } ] (* Solving all the eight equations and generating the coefficients *) sol = Solve[ Equn, Table[a[i], {i, 1, 60} ]][[1]] [ [ 1] ]是指第一个解，变量是a[i]，用Table列出每一个变量。Equn里每一行都有a[i]的系数，i最大为45，所以Table里i从1到至少45。 (* Substituting the values of coefficients back into the equations *) poly1 = poly[order] / . sol 方程的解的形式为sol，将计算结果替换原方程未知数 {a[1]→0,a[2]→0,a[3]→0,……} 上式中是将已经求出来的a[i]的解替换原来变量。解的形式为a[1]→0等。 (* Simplifying the equation by taking out common terms *)

应用伽辽金法求解一维热传导问题

应用伽辽金法求解一维热传导问题 解：令~(1)1n i a x x i i Φ≈Φ=?-∑= i a 为待定系数，形函数为(1)i i N x x =-，近似解~ Φ满足边界条件，无余量，近似解~ Φ在域中连续。 ?21201[()]0n j i i x i d w N a Q d x =+=∑? 由伽辽金法，j j w N =，对上式进行分步积分，有 11 001()0j n w j x j j x i x x d d w Qd N a d d d =-=∑?? 当i=1时， 211(1)w N x x x x ==-=- 上式变为1122100()[()](21)0x x x d x x d a x x x d d ---?-=?? 求得114a =- 则一项解为~21()4x x =--Φ 当i=2时， 211222(1)w N x x w N x x ==-==- 原式变为 1 1222120011322221200()[(1)(1)](21)0()[(1)(1)](32)0x x x x x x d x x d a x x a x x x d d d x x d a x x a x x x x d d ?---+-?-=????---+-?-=?? ???? 即变为

12121113612125 6 15192a a a a ?+=-????+=-?? 121332516a a ?=-?????=?? 则二项解为~2323213552313()()3216163232 x x x x x x x =--+-=-+Φ 而本例的精确解为 2131,0282111,1882x x x x x ?Φ=-+≤≤????Φ=-+<≤?? 如图所示，给出伽辽金法的近似一项解，二项解与精确解结果的曲线。

有限元法的基础理论

一、里兹法与迦辽金法（摘自电磁场有限元方法 金建铭） 1. 里兹法 里兹法是一种变分方法，其中边值问题用变分表达式（也称泛函）表示，泛函的极小值对应于给定边界条件下的控制微分方程。通过求泛函相对于其变量的极小值可得到近似解。 2. 伽辽金法 伽辽金法属于残数加权方法类型，它通过对微分方程的残数求加权的方法得到方程的解。 若u 是方程的近似解，将u 代入方程可得到非零的残数： r Lu f =- u 的最佳近似应能使残数r 在Ω内所有点上有最小值。残数加权方法要求： 0i i R rd ωΩ =Ω=? 这里i R 表示残数的加权积分，i ω是所选的加权函数。 在伽辽金法中，加权函数与近似解展开中所用的函数相同。通常，这样可得到最精确的 解。 二、有限元方法 里兹法和伽辽金法中，在整个解域内找出能表示或至少近似表示问题真实解的试探函数是非常重要的。然而对于许多问题，这个步骤是十分困难的，对二维和三维问题尤其如此。为此，我们可将整个区域划分成小子域，并应用定义在每个子域上的试探函数。因为子域是小区域，因而在每一子域内函数的变化不大，所以定义在子域上的试探函数通常比较简单。这正是有限元法的基本思想。应用里兹法的过程通常称为里兹有限元法或变分有限元法，而应用伽辽金方法的过程通常称为伽辽金有限元方法。 有限元法与经典里兹法和伽辽金法的不同之处是在试探函数的公式上。在经典里兹法和伽辽金法中，试探函数由定义在全域上的一组基函数组成。这种组合必须能够（至少近似）表示真实解，也必须满足适当的边界条件。在有限元法中，试探函数是由定义在组成全域的子域上的一组基函数构成。因为子域很小，所以定义在子域上的基函数能够十分简单。 三、关于形函数（摘自有限元法在电磁计算中的应用 张榴晨） 对于一个待求的微分方程，用一组线性独立的尝试函数i ψ和待定系数i C 来表示方程的近似解，并用加权余数法（迦辽金法）来求解这些待定系数。求解待定系数的代数方程组为： 1 []1,2,,n i j i j i d C q d j n ψψψΩ Ω =??Ω=Ω =∑? ? 这里j ψ为所选择的加权函数，应用迦辽金法时，所选取的加权函数即为尝试函数。 有限元中应用的尝试函数代表了单元上近似解的一种插值关系，它决定了近似解在单元上的形状。因此尝试函数在有限元法中又称为形函数。对于一维有限元来说，形函数为一个直线段；对一维高阶有限元来说，形函数为一个曲线段；对二维一阶有限元来说，形函数为一个平面；对二维高阶有限元来说，形函数为一个曲面；三维有限元来说，形函数为多维平面或曲面。选择形函数时可以使一个任意元上的函数只与该元所对应的节点势函数值有关，而与其它各点的值无关。 1. 一维有限元