全矩阵无积分间断有限元方法实现

非结构网格下2D Riesz分数阶方程的Galerkin有限元方法卜玮平【摘要】讨论了2D Riesz分数阶扩散方程的Galerkin有限元方法. 基于非结构网格,采用Lagrange线性分片多项式作为基函数,详细描述了分数阶扩散方程的有限元实现. 与现有方法相比, 该方法有效地降低了计算成本, 提高了刚度矩阵的精度.最后,数值算例验证了所提方法的有效性.【期刊名称】《纯粹数学与应用数学》【年(卷),期】2019(035)002【总页数】13页(P169-181)【关键词】Riesz分数阶扩散方程;有限元方法;非结构网格【作者】卜玮平【作者单位】湘潭大学数学与计算科学学院,湖南湘潭 411105【正文语种】中文【中图分类】O1781 引言最近,分数阶微分方程数值方法的研究越来越受关注.有限差分方法(FDM)是求解分数阶微分方程最常见的数值方法之一.文献[1]是研究FDM求解分数阶微分方程数值方法的先驱.随后,利用有限差分方法人们对各种分数阶微分方程进行了求解[2-5].谱方法(SM)和谱元方法(SEM)是求解分数阶微分方程的重要方法,它们通常具有很高的收敛精度.目前,SM和SEM求解分数阶微分方程的工作有文献[6-8]等.相比于FDM方法和SM(或SEM)方法,有限元方法(FEM)的主要特征是它能够较易处理复杂区域且对解的光滑性要求较低.文献[9]首次尝试用有限元方法求解分数阶微分方程,并给出了有限元逼近的理论框架.此后,文献[10]建立了有限元求解一维和二维分数阶对流弥散方程的理论框架.随后,在有限元求解分数阶微分方程方面涌现了越来越多的工作,其中包括文献[11-15]等.近年来,关于二维分数阶微分方程的有限元方法也有一些研究成果.文献[10]考虑了基于分数阶方向导数的二维分数阶对流-弥散方程的有限元方法.通过使用矩阵转换技术,文献[16]讨论了基于分数阶拉普拉斯算子的二维分数阶扩散方程有限元方法.文献[17]发展了基于分数拉普拉斯算子的分数阶扩散方程自适应有限元方法.为了求解二维Risez/Riemann-Liouville分数阶扩散方程(2DRFDE/2DRLFDE),基于一致三角网格剖分,文献[18-19]考虑了Galerkin有限元方法.然而,对于不规则区域上2DRFDE/2DRLFDE的有限元方法,通常需要采用非结构网格.基于非结构网格,文献[20]利用三角形上定义的Lagrange多项式,建立了2DRLFDE的间断Galerkin方法.文献[21]考虑了非线性2DRFDE的有限元方法,并描述了有限元方法的实现.文献[22]讨论了二维时空分数阶波方程在非规则凸域上的有限元方法.尽管在文献[20-22]中已有非结构网格下有限元求解分数阶微分方程的工作,但是这些工作使用的是相同的有限元实现技巧.值得注意,上述文中提到的有限元实现方法有不足之处.因此,在这篇文章中将提出一些新的技巧,以改进现有的有限元实现方法.本文的主要贡献如下:首先,对于刚度矩阵的计算降低了计算花销,原来的计算花销为O(Ne3),现在的花销为O(Ne2),Ne为总剖分单元数;其次,提高了刚度矩阵的元素的精度,对于三角形单元的高斯积分,不同于现有的计算方法,高斯积分区域为被积函数的非零的区域,这比由以往方法得到的刚度矩阵元素更精确;第三,将Riemann-Liouville导数转化为Caputo导数,简化了内积的计算.本文的结构安排如下:在第2节,给出了分数阶导数的定义、模型问题及有限元全离散格式的推导;在第3节,首先介绍了现有的刚度矩阵计算方法,然后详细描述了有限元的实现方法,并与现有方法进行了比较;在第4节,给出了数值实验来证明方法的有效性;最后,对本文进行了总结.2 准备工作令Ω⊂R2,则x,y方向的左右Riemann-Liouville分数阶导数定义如下:其中γ,n−1<γ≤n,n∈N,a(y),b(y),c(x),d(x)定义如图1.图1 关于Ω上a(y),b(y),c(x),d(x)的定义进一步,定义x,y方向的γ(γ1,3,···)阶 Riesz分数阶导数如下类似地,x,y方向γ阶右Caputo分数阶导数定义如下在文献[23]中已经讨论了Riemann-Liouville分数阶导数和Caputo分数阶导数的等价关系,如果u(a(y),y)=0,u(x,c(x))=0,0<γ<1,则本文考虑如下分数阶扩散方程其中,β<1,P,Q,A为常数.对上述方程,首先考虑其变分形式和有限元全离散格式.对于任意γ≥0,记γ(Ω)为Hγ(Ω)的子集且其元素在Ω 外的零扩张属于Hγ(R2).令V=α(Ω)∩ β(Ω).于是,由文献[18]中引理5可得如下变分问题:寻找u∈V使得其中(·,·):=(·,·)L2(Ω),F(v):=(f,v),且为了得到上述变分问题的全离散格式,先将Ω进行剖分.令{Th}是Ω的一个正则的三角剖分,h为所有三角形单元的最大直径.定义如下有限元空间这里PK(x,y)为K次多项式,接下来,定义(10)式的全离散格式:寻找uh∈XhK使得3 有限元方法的实现为了计算全离散格式(14),需要计算刚度矩阵和荷载向量.由于分数阶导数为非局部算子,与整数阶微分方程相比,其主要区别在于分数阶微分方程刚度矩阵的计算非常复杂.因此,方法实现的重点将放在刚度矩阵的计算上.令 ,则uh可写成,其中ϕi(x,y)是属于剖分节点i的Lagrange线性基函数,N为剖分节点总数目.如果三角形单元e包含顶点i且它的三个顶点逆时针排列依次为i,j,k,则有这里∆e为三角形单元e的面积,(xi,yi),(xj,yj),(xk,yk)是对应顶点i,j,k的坐标.接下来,考虑B(uh,vh)的计算.令这里S={Sij}N×N为刚度矩阵,其元素为3.1 已有的实现方法目前,文献[20-22]已经考虑了Sij的计算.然而,应该注意这些文章关于Sij的计算方法是类似的.这里简单对文献[20-22]中Sij的计算方法进行描述.由于(17)中四个内积的计算具有相似性,因此仅仅以为例进行讨论.考虑分数阶算子的非局部性并利用高斯积分可得这里Ge为三角形单元e上的高斯点,ωl是高斯点(xl,yl)对应的权重.上面(18)式中的方法可以用来计算Sij,然而该方法有三点不足之处.令Ne为剖分三角形单元的总数.首先,对于Sij的计算需要考虑如下积分在每个三角形单元的高斯积分因此,为了得到刚度矩阵的一个元素Sij,上述积分需要进行Ne次,而为了获得整个刚度矩阵,需要计算次,这将使得刚度矩阵的计算花销随Ne的增涨而快速增涨.其次,对于下列积分考虑被积函数在三角形单元上的非零区域.由于分数阶算子具有非局部性,显然存在满足下列情形的三角形单元:被积函数在三角形的某些部分为零,在其他部分非零.因此,被积函数在满足上述条件的三角形单元上为间断函数.如果对这些单元,在整个三角形上运用高斯积分势必达不到数值积分的相应精度.图2给了描述上述情形的例子(图形(a)描述了的非零区域,其中ϕi为节点i的基函数,ϕk为节点k的基函数;图形(b)描述了被积函数的非零区域,它为支集的交集).图2 的非零区域和的非零区域最后,对于被积函数 ,将说明它在某些区域光滑性差的特点.假设ϕi(x,y)和点(xp,yp)由图3给出在(xp,yp)的计算涉及的三角形的边界点p0,p1,p2,它们为积分路径与三角形单元e1,e2,e3边界的交点,它们在x轴方向的坐标分别设为x0,x1,x2).图3 在(xp,yp)的计算结果计算ϕi(x,y)的左Riemann-Liouville分数阶导数在(xp,yp)点的值从公式(21)可以看出, 在支集ϕi(x,y)的某些点处光滑性较差.类似地,也可以证明存在光滑性较差的区域.因此,为了保持(19)式的精度,计算时需要取大量的高斯点. 3.2 实现方法针对已有方法的不足,设计一种新的求解二维Riesz分数阶扩散方程的有限元实现方法.因为,显然因此,为了计算(16)中刚度矩阵S,根据对称性仅需要计算Sx和Sy.下面以为例来描述计算方法.假设ϕi(x,y),ϕj(x,y)的支集定义如图4(Lagrange线性基函数ϕi(x,y),ϕj(x,y)的支集如图4,其中为ϕi(x,y)支集上的三角形单元).图4 Lagrange线性基函数ϕi(x,y),ϕj(x,y)的支集这里是ϕi(x,y)关于变量x的系数.比较 (26)式与 (18)式,显然的计算花销将要下降.为了得到刚度矩阵,此时仅仅需要计算次积分.为了计算积分Ii,已有的方法是在整个三角形单元ei上使用高斯积分.然而,上面已经提到这样做将降低高斯积分的效果,因为分数阶算子具有非局部性, 的非零区域在某些单元可能只占有一部分(即在该单元为间断函数),如图5.图5 的支集的支集包含三角形单元因此,I1,I2能够在上选择高斯点来计算.然而, 仅仅含有支集的一部分.因此,可得其中为在上的支集.由于是一个三角形区域,因此在其上可以选用高斯点来计算I3.由于是四边形区域,因此要计算I4,I5,一种方法是在四边形上选用高斯点,另一种方法是将四边形分解成两个三角形然后分别作高斯积分.考虑图6所描述的情形,即是一个多边形区域.因此对I5的计算,一种方法是将该区域剖分成三角形与四边形然后进行高斯积分,另一种方法可以将积分区域全部剖分成三角形在每个三角形上进行高斯积分.图6 的支集注意到,对于ϕi(x,y), 由 (9)式可得于是利用高斯积分有这里e为的非零区域,ai,∆e定义在(15)式中.假设ϕj(x,y)及其积分路径被定义在图3中,对于高斯点(xp,yp),有显然(30)式比(21)式更容易计算.4 数值实验本节给出一个数值例子来验证方法的有效性.例 4.1在模型方程(10)中,取P=Q=2,A=4.考虑如下两种情形:(a)假设考虑问题区域为[0,1]×[0,1],其精确解为u=100x2(1−x)2y2(1−y)2;(b)假设考虑问题区域为,其精确解为,相应的右端函数分别定义如下当选择不同的α和β时,表1-表4分别列出了情形(a)和情形(b)所得的数值结果. 可以看出,所得误差的收敛率是最优收敛率.表1 基于Lagrange线性多项式与α=0.6,β=0.6计算情形(a)的数值误差与收敛率h‖uh −u (x,y) ‖ 0 收敛率1 4 2.193 4e-2 –1 8 4.730 6e-3 2.213 1 1 16 1.106 0e-3 2.096 7 1 32 2.489 8e-4 2.151 2表2 基于Lagrange线性多项式与α=0.7,β=0.8计算情形(a)的数值误差与收敛率h‖uh −u (x,y) ‖ 0 收敛率1 4 2.259 1 e− 2 –1 8 4.998 3 e− 3 2.176 2 1 16 1.173 6 e− 3 2.090 5 1 32 2.906 8 e− 4 2.013 4表3 基于Lagrange线性多项式与α=0.7,β=0.6计算情形(b)的数值误差与收敛率h‖uh −u (x,y) ‖ 0 收敛率1 4 4.567 0 e− 3 –1 8 1.106 5 e− 3 2.045 2 1 16 2.374 1 e− 4 2.220 6 1 32 5.610 1 e− 5 2.081 3表4 基于Lagrange线性多项式与α=0.6,β=0.8计算情形(b)的数值误差与收敛率h‖uh −u (x,y) ‖ 0 收敛率1 4 4.712 8 e− 3 –1 8 1.133 7 e− 3 2.055 5 1 16 2.467 7 e− 4 2.199 8 1 32 5.767 4 e− 5 2.097 25 总结本文研究了非结构网格下利用Lagrange线性基函数求解2D Riesz分数阶扩散方程的有限元方法实现.首先,描述了现有有限元全离散格式的实现方法,并指出了现有方法的不足之处.随后,针对这些缺点设计了一种新的实现方法,提高了有限元方法的计算效率和刚度矩阵的精度.最后,给出了数值算例,数值结果验证了本文所提方法的有效性.参考文献【相关文献】[1]Lubich C.Discretized fractional calculus[J].SIAM J.Math.Ana.,1986,17(3):704-719.[2]Liu F,Zhuang P,Anh V,et al.Stability and convergence of the difference methods for the space-time fractional advection-diffusion equation[J]pu.,2007,191(1):12-20.[3]Sun Z,Wu X.A fully discrete difference scheme for a diffusion-wave system[J].Applied Numerical Mathematics,2006,56(2):193-209.[4]Ding Z,Xiao A,Li M.Weighted finite difference methods for a class of space fractional partial differential equations with variablecoefficients[J].Appl.Math.,2010,233(8):1905-1914.[5]Tian W Y,Zhou H,Deng W.A class of second order difference approximations for solving space fractional diffusion equations[J].Mathematics of Computation,2015,84(294):1703-1727.[6]Zayernouri M,Karniadakis G E.Fractional Sturm-Liouville eigen-problems:theory and numerical approximation[J].Journal of Computational Physics,2013,252:495-517.[7]Li X,Xu C.Existence and uniqueness of the weak solution of the space-time fractional diffusion equation and a spectral method approximation[J].Phy.,2010,8(5):1016-1021.[8]Zeng F,Liu F,Li C,et al.A Crank–Nicolson ADI spectral method for a two-dimensional Riesz space fractional nonlinear reaction-diffusion equation[J].SIAMJ.Num.Ana.,2014,52(6):2599-2622.[9]Fix G J,Roof J P.Least squares finite-element solution of a fractional order two-point boundary value problem[J].Computers&Mathematics with Applications,2004,48(7-8):1017-1033.[10]Roop J P.Variational solution of the fractional advection dispersion equation[D].South Carolina:Clemson University,2004.[11]Wang H,Yang D.Wellposedness of variable-coefficient conservative fractional elliptic differential equations[J].SIAM Journal on Numerical Analysis,2013,51(2):1088-1107. [12]Jin B,Lazarov R,Pasciak J,et al.Error analysis of a finite element method for the space-fractional parabolic equation[J].SIAM Journal on Numerical Analysis,2014,52(5):2272-2294.[13]Deng W.Finite element method for the space and time fractional Fokker-Planck equation[J].SIAM journal on numerical analysis,2008,47(1):204-226.[14]Li C,Zhao Z,Chen Y Q.Numerical approximation of nonlinear fractional differential equations with subdiffusion and superdiffusion[J].Compu.Math.Appl.,2011,62(3):855-875.[15]Ma J,Liu J,Zhou Z.Convergence analysis of moving finite element methods for space fractional differential equations[J].Journal of Computational and Applied Mathematics,2014,255:661-670.[16]Yang Q,Turner I,Liu F,et al.Novel numerical methods for solving the time-space fractional diffusion equation in two dimensions[J].SIAM .,2011,33(3):1159-1180.[17]Chen L,Nochetto R H,Otrola E,et al.A PDE approach to fractional diffusion:a posteriori error analysis[J].Journal of Computational Physics,2015,293:339-358.[18]Bu W,Tang Y,Yang J.Galerkin finite element method for two-dimensional Riesz space fractional diffusion equations[J].Journal of Computational Physics,2014,276:26-38. [19]Zhao Y,Bu W,Huang J,et al.Finite element method for two-dimensional space-fractional advection-dispersion equations[J].Applied Mathematics andComputation,2015,257:553-565.[20]Qiu L,Deng W,Hesthaven J S.Nodal discontinuous Galerkin methods for fractional diffusion equations on 2D domain with triangular meshes[J].Phy.,2015,298:678-694.[21]Yang Z,Yuan Z,Nie Y,et al.Finite element method for nonlinear Riesz space fractionaldiffusion equations on irregular domains[J].Journal of ComputationalPhysics,2017,330:863-883.[22]Fan W,Liu F,Jiang X,et al.A novel unstructured mesh finite element method for solving the time-space fractional wave equation on a two-dimensional irregular convex domain[J].Fractional Calculus and Applied Analysis,2017,20(2):352-383.[23]Bu W,Tang Y,Wu Y,et al.Finite difference/ finite element method for two-dimensional space and time fractional Bloch-Torrey equations[J].Journal of Computational Physics,2015,293:264-279.。

riks算法原理
RIKS（Reduced Integration with Known Stiffness）算法是一种用于非线性分析的有限元方法。

它在处理非线性材料和几何非线性时，相比传统的全积分有限元方法具有更高的计算效率。

RIKS算法的基本原理是利用已知刚度矩阵进行迭代求解非线性问题。

以下是RIKS算法的主要步骤：
1、定义初始状态：给定初始位移（或应变）和载荷条件，初始化位移（或应变）和载荷向量。

2、计算刚度矩阵：根据当前位移（或应变）状态计算刚度矩阵。

这可以通过线性化处理或使用增量刚度矩阵方法来实现。

3、求解增量位移（或应变）：利用已知刚度矩阵和外部载荷向量，求解增量位移（或应变）。

4、更新位移（或应变）：将增量位移（或应变）叠加到当前位移（或应变）上，得到更新后的位移（或应变）。

5、判断收敛：检查收敛标准，如果满足要求则停止迭代，否则返回第3步继续迭代。

RIKS算法的核心思想是通过在每个加载步骤中只考虑刚度矩阵的部分积分，从而减少计算量。

在每个迭代步骤中，RIKS算法使用一个称为“刚度比例因子”的参数来控制刚度矩阵的缩放。

这个刚度比例因子可以根据问题的特性和收敛行为进行调整。

RIKS算法也可用于处理失稳问题，例如杆件的屈曲分析。

通过在加载过程中引入惯性力项，可以模拟出杆件局部稳定性失效的情况。

总之，RIKS算法通过使用已知刚度矩阵和部分积分的方法，在非线性问题求解中提供了一种高效的选择。

它可以减少计算量并加快收敛速度，对于大型和复杂的非线性问题具有重要的应用价值。

南京理⼯⼤学研究⽣有限元⽅法理论及应⽤考试个⼈答案⽬录1等参单元及其应⽤ (1)1.1概述 (1)1.1.1等参单元的概念、原理 (1)1.1.2等参单元对有限元法⼯程应⽤的意义 (1)1.2等参单元的数值积分⽅法 (1)1.2.1等参单元刚度矩阵的数值积分⽅法 (1)1.2.2确定积分阶的原理 (2)1.2.3全积分单元与减缩积分单元讨论和评价 (3)1.3线性等参单元 (3)1.3.1全积分、减缩积分线性等参单元有关问题的分析讨论 (3)1.4等参单元的应⽤ (5)2分析与计算 (6)2.1四节点平⾯等参单元的收敛协调性 (6)2.2⼋节点平⾯等参单元 (8)2.33节点平⾯三⾓形单元 (9)2.420节点六⾯体等参单元 (10)2.520节点六⾯体等参单元 (11)3上机实验 (15)3.1实验⼀ (15)3.1.1实验题⽬ (15)3.1.2实验⽬的 (15)3.1.3建模概述 (15)3.1.4计算结果分析与结论 (16)3.1.5实验体会与总结 (32)3.2实验⼆ (33)3.2.1实验题⽬ (33)3.2.2实验⽬的 (33)3.2.3建模概述 (33)3.2.4计算结果分析与讨论 (34)3.2.5实验体会与总结 (36)3.3实验三 (36)3.3.1实验题⽬ (36)3.3.2实验⽬的 (36)3.3.3建模概述 (37)3.3.4计算结果分析与结论 (37)3.3.5实验体会与总结 (44)1 等参单元及其应⽤1.1 概述1.1.1 等参单元的概念、原理普通单元受到两个⽅⾯的限制：（1）单元的精度。

单元的节点数越多，单元精度越⾼；（2）单元⼏何上的限制。

普通矩形和六⾯体单元都不能模拟任意形状⼏何体，所有⼏种普通单元都是直线边界，处理曲边界⼏何体误差较⼤。

为了解决上述⽭盾，⽅法就是突破矩形单元和六⾯体单元⼏何⽅⾯的限制，使其成为任意四边形和任意六⾯体单元，这类单元位移模式和形函数的构造和单元列式的导出不能沿⽤构造简单单元的⽅法，必须引⼊等参变换，采⽤相同的插值函数对单元的节点坐标和节点位移在单元上进⾏插值。

等参单元及其应用摘要本文主要讲述等参单元的原理及其对有限元法工程应用的意义。

等参单元的数值积分方法，等参单元刚度矩阵的数值积分方法及确定积分阶的原理。

全积分、减缩积分单元讨论和评价。

线性等参单元和非协调元，全积分、减缩积分线性等参单元和非协调元有关问题的分析讨论。

关键词等参单元; 数值积分; 应用1.引言用有限元法划分单元时，单元的节点数越多，单元精度越高。

因此在这一点上，矩形单元优于简单三角形单元，六面体单元优于四面体单元。

但单独使用矩形或长方体单元都不能模拟任意形状几何体，且网格中单元大小无法过渡。

所有上述单元都是直线边界，处理曲边界几何体误差较大。

解决上述矛盾的途径是突破矩形单元和长方体单元几何上的限制，使其成为平面任意四边形和空间任意六面体单元，如果再增加边中间节点，还可以成为曲边四边形和曲面六面体高精度单元。

任意四边形和任意六面体单元的位移模式和形函数的构造不能沿用前面构造简单单元时采用的总体坐标多项式位移函数插值的方法，必须通过所谓的等参变换建立单元局部坐标，采用相同的插值函数对单元节点的总体坐标和节点位移在单元上进行插值。

这类单元称为等参单元。

等参单元的提出对于有限元法在工程实践中的应用具有重要意义。

2.等参单元的数值积分方法2.1 高斯数值积分的基本概念一维高斯数值积分公式：i ni i H x f dx x f I )()(111∑⎰=-== 其中：积分点-i x ，积分点数目，积分阶-n ，权重系数-i H结论：n 阶高斯积分公式对 2n-1 次多项式被积函数可求得精确积分！ 同理，对二维高斯积分：),(),(111111i i j n i nj i F H H d d F I ηξηξηξ∑∑⎰⎰==--==积分公式对ξ，η方向最高方次为 2n-1 的多项式可求得精确值。

2.2 减缩积分的原理实际应用中选取的积分阶往往可以低于被积函数所有项次精确积分所需要的阶数，这种积分方案称为减缩积分。

有限元方法的求解步骤
1.构建几何模型：首先，需要根据实际问题构建一个几何模型。

这可以通过使用计算机辅助设计（CAD）软件进行建模，或者手动绘制模型。

2.离散化：在几何模型的基础上，需要将其离散化为有限个小元素。

最常用的元素是三角形和四边形，也可以使用更复杂的元素类型。

3.选择数学模型和假设：根据问题的物理特性，需要选择适当的数学模型和假设。

这可能涉及选择适当的方程、边界条件和材料性质等。

4.导出有限元方程：根据选择的数学模型和假设，使用变分原理或其他数学方法，可以导出与离散化模型相对应的有限元方程。

这个方程通常是一个代数方程组。

5.建立刚度矩阵和负载向量：有限元方程可以转化为刚度矩阵和负载向量的形式。

刚度矩阵描述了系统中元素和节点之间的关系，而负载向量描述了外部作用力。

6.施加边界条件：为了解决方程组并确定未知位移，需要施加边界条件。

边界条件可以是位移约束、力约束或其他类型的约束。

7.求解方程：将刚度矩阵和负载向量与边界条件组合起来，可以形成一个线性代数方程组。

可以使用各种数值方法求解线性方程组，例如直接求解、迭代法、预处理方法等。

8.后处理：在求解方程后，可以根据需要进行后处理。

后处理包括计算和输出感兴趣的结果，如应力、位移、应变等。

9.验证和调整：完成有限元求解后，需要验证结果的准确性，并根据需要对模型参数进行调整。

验证可以通过与理论解、实验结果或其他数值方法进行比较来完成。

10.进行优化和设计：利用有限元模拟的结果，可以进行系统的优化和设计改进。

这可以通过改变几何形状、材料属性或边界条件来实现。

有限元复习资料工程问题的研究对象：离散体：材料力学中的连续梁、建筑结构框架和桁架结构。

我们把这类问题，称为离散系统。

离散问题是可解的连续体：通常可以建立它们应遵循的基本方程，即微分方程和相应的边界条件。

通常只能得到少数简单问题的精确解答。

对于许多实际的工程问题，还无法给出精确的解答。

有限元方法的两大应用:科学计算或数值模拟数字设计或虚拟仿真(CAD CAM CAE)有限元法的实质：将复杂的连续体划分为有限个简单的单元体，化无限自由度问题为有限自由度问题。

将连续场函数的（偏）微分方程的求解问题转化为有限个参数的代数方程组的求解问题。

有限元法的基本思路：将连续系统分割成有限个分区或单元，对每个单元提出一个近似解，再将所有单元按标准方法组合成一个与原有系统近似的系统有限元分析的基本步骤：1）建立研究对象的近似模型。

2）将研究对象分割成有限数量的单元(结构离散化)3）用标准方法对每个单元提出一个近似解（单元分析）（1）选择位移模式（2）建立单元刚度矩阵（3）计算等效节点力4）将所有单元按标准方法组合成一个与原有系统近似的系统（单元集成）5）用数值方法求解这个近似系统。

（选择合适计算方法计算节点位移、应力、应变）6）计算结果处理与结果验证（如何显示、分析数据并找到有用的结论）弹性力学的基本假设：（1）连续性（2）均匀性（3）各向同性（4）完全弹性符合（1）-（4）假定的称为理想弹性体。

（5）小变形假定满足以上五个基本假设的弹性力学称为线弹性力学。

弹性力学中的基本量：1）外力（体力面力集中力）弹性体受外力以后，其内部将产生应力。

三种外力应力（正应力剪应力）应力成对出现，坐标面上的应力以正面正向，负面负向为正。

六个应力分量剪应力互等定律：作用在两个互相垂直的面上并且垂直于该两面交线的剪应力是互等的。

(大小相等，正负号也相同)。

弹性体在受外力以后，还将发生变形。

物体的变形状态，一般有两种方式来描述：1、给出各点的位移；（运动）三个位移分量2、给出各微元体的变形（应变）六个应变分量在弹性体内部，三类基本方程：根据微分体上力的平衡条件（法向应力和剪应力即内力与其体力平衡），，建立平衡微分方程；根据微分线段上应变－位移的几何条件，建立几何方程；根据应力－应变间的物理条件，建立物理方程当沿X轴方向的两个对面受有均匀布的正应力时，弹性体在X方向的伸长还伴随有侧收缩，即在y和Z方向的单位缩短。

题目：使用MATLAB编程实现有限元切线刚度矩阵计算一、引言有限元法是一种用于求解复杂工程问题的数值分析方法，它将连续介质划分为许多小的单元，通过对每个单元进行离散化处理，再用数学方法对这些单元进行组装，最终得到整个结构的解。

在有限元方法中，刚度矩阵是求解结构问题的关键步骤之一，而有限元切线刚度矩阵的计算则是其中的重要内容之一。

本文将介绍如何使用MATLAB编程实现有限元切线刚度矩阵的计算。

二、有限元切线刚度矩阵的基本概念1. 切线刚度矩阵在有限元方法中，切线刚度矩阵是描述结构对外部载荷作用下的应变-应力关系的重要矩阵。

它描述了结构在外部载荷下的变形行为，是求解结构变形和应力的重要工具。

2. 切线刚度矩阵的计算切线刚度矩阵的计算是通过对单元的局部坐标系进行刚度矩阵的求解，并进行坐标变换得到全局坐标系下的切线刚度矩阵。

在实际计算中，需要考虑单元的几何形状、材料性质等因素，以及在单元上施加的外部载荷。

三、MATLAB编程实现有限元切线刚度矩阵的基本步骤1. 单元刚度矩阵的计算我们需要编写MATLAB函数来实现对单元刚度矩阵的计算。

这个函数需要考虑单元的几何形状、材料性质等因素，以及在单元上施加的外部载荷。

通常情况下，我们可以利用数值积分的方法来进行刚度矩阵的计算。

2. 坐标变换矩阵的计算在得到单元刚度矩阵之后，我们需要计算坐标变换矩阵，将单元刚度矩阵从局部坐标系变换到全局坐标系。

这也需要编写一个MATLAB函数来实现坐标变换矩阵的计算。

3. 矩阵组装我们需要将所有单元的切线刚度矩阵组装成整体刚度矩阵。

这通常需要考虑到单元之间的连接关系，以及边界条件等因素。

在MATLAB中，我们可以利用矩阵的组合和相加等运算来实现整体刚度矩阵的计算。

四、编程实例这里我们以一个简单的弹簧-弹簧系统为例，介绍如何使用MATLAB编程实现有限元切线刚度矩阵的计算。

我们需要定义系统的几何形状、材料性质等参数，然后编写MATLAB函数来进行刚度矩阵的计算，坐标变换矩阵的计算，以及矩阵的组装，最终得到整体刚度矩阵，并求解系统的变形和应力。

接触问题求解方法(2010-10-02 00:25:00)转载▼标签：杂谈不管在接触边界之间是否有间隙存在，接触作用的出现对结构受载荷之后的接触状态和应力分布都有直接的影响，有些结构正是由于接触作用，使不连续的部分共同工作，从而提高了整个结构的承载力和刚度；而正是由于接触的存在使得有些结构出现局部高应力，很容易使材料屈服或发生裂缝，如果再受到循环荷载的影响，还可能产生疲劳失效。

因此对于接触问题的研究具有重要的工程实际意义。

赫兹接触理论及后来其它学者发展的弹性接触理论称为经典接触力学，它们都是通过封闭的解析解来解决接触问题的，但其应用范围有限，因而接触力学的进展主要与消除这些限制有关。

接触问题属于数学上的混合边值问题，Boussinesq积分方程是其主导方程。

按所用数学方法的不同，可将接触问题的理论解法大致分为经典解法和非经典解法。

以经典的数学工具如积分变换法和复变函数法求解接触问题的方法称为经典解法，以有限元法、边界元法等求解接触问题的方法称为非经典法。

以传统有限元法等数值方法为基础的非经典解法主要有罚函数法、拉格朗日乘子法、增广拉格朗日乘子法、摄动拉格朗日乘子法及数学规划法。

Trefftz有限元模型利用辅助网线位移场或面力场，在一种杂交意义上将单元域内位移场关联起来。

单元域内位移场精确满足控制微分方程，它可表达为微分方程的特解、适当截断的Trefftz完备解系与待定参数乘积的和的形式。

利用定义在每个单元边界上的独立的网线位移场，单元间的连续性就在一种近似意义上得到满足。

在单元一级上消去内部待定参数后即可得到标准的力-位移关系式（即单元刚度方程）。

变分泛函是Trefftz有限元法的核心，它在单元公式推导中起着至关重要的作用。

由于Trefftz有限元法继承了传统有限元法和传统边界元法的优点：(1)公式中只含有对单元的边界积分，这样就可以生成任意多边形单元甚至曲边单元。

因此，Trefftz有限元法可被认为是一种特殊形式的边界型求解方法，单元的边界类似于一种特殊形式的边界单元，其刚度矩阵对称、计算简单，而不像传统边界元那样要进行复杂的计算（如计算边界奇异积分的复杂积分规则，处理非齐次方程的特殊积分，处理间断问题的双节点技术等等）。

一 、20分）（×）1. 节点的位置依赖于形态，而并不依赖于载荷的位置（√）2. 对于高压电线的铁塔那样的框架结构的模型化处理使用梁单元（×）3. 不能把梁单元、壳单元和实体单元混合在一起作成模型（√）4. 四边形的平面单元尽可能作成接近正方形形状的单元（×）5. 平面应变单元也好，平面应力单元也好，如果以单位厚来作模型化处理的话会得到一样的答案（×）6. 用有限元法不可以对运动的物体的结构进行静力分析（√）7. 一般应力变化大的地方单元尺寸要划的小才好（×）8. 所谓全约束只要将位移自由度约束住，而不必约束转动自由度（√）9. 同一载荷作用下的结构，所给材料的弹性模量越大则变形值越小（√）10一维变带宽存储通常比二维等带宽存储更节省存储量。

二、填空（20分）1．平面应力问题与薄板弯曲问题的弹性体几何形状都是 薄板 ，但前者受力特点是： 平行于板面且沿厚度均布载荷作用 ，变形发生在板面内；后者受力特点是： 垂直于板面 的力的作用，板将变成有弯有扭的曲面。

2．平面应力问题与平面应变问题都具有三个独立的应力分量： σx ，σy ，τxy ，三个独立的应变分量：εx ，εy ，γxy ，但对应的弹性体几何形状前者为 薄板 ，后者为 长柱体 。

3．位移模式需反映 刚体位移 ，反映 常变形 ，满足 单元边界上位移连续 。

4．单元刚度矩阵的特点有：对称性 ， 奇异性 ，还可按节点分块。

5．轴对称问题单元形状为：三角形或四边形截面的空间环形单元 ，由于轴对称的特性，任意一点变形只发生在子午面上，因此可以作为 二 维问题处理。

6．等参数单元指的是：描述位移和描述坐标采用相同的形函数形式。

等参数单元优点是：可以采用高阶次位移模式，能够模拟复杂几何边界，方便单元刚度矩阵和等效节点载荷的积分运算。

7．有限单元法首先求出的解是 节点位移 ，单元应力可由它求得，其计算公式为{}{}[][]eD B σδ=。

目录1等参单元及其应用 (1)1.1概述 (1)1.1.1等参单元的概念、原理 (1)1.1.2等参单元对有限元法工程应用的意义 (1)1.2等参单元的数值积分方法 (1)1.2.1等参单元刚度矩阵的数值积分方法 (1)1.2.2确定积分阶的原理 (2)1.2.3全积分单元与减缩积分单元讨论和评价 (3)1.3线性等参单元 (3)1.3.1全积分、减缩积分线性等参单元有关问题的分析讨论 (3)1.4等参单元的应用 (5)2分析与计算 (6)2.1四节点平面等参单元的收敛协调性 (6)2.2八节点平面等参单元 (8)2.33节点平面三角形单元 (9)2.420节点六面体等参单元 (10)2.520节点六面体等参单元 (11)3上机实验 (15)3.1实验一 (15)3.1.1实验题目 (15)3.1.2实验目的 (15)3.1.3建模概述 (15)3.1.4计算结果分析与结论 (16)3.1.5实验体会与总结 (32)3.2实验二 (33)3.2.1实验题目 (33)3.2.2实验目的 (33)3.2.3建模概述 (33)3.2.4计算结果分析与讨论 (34)3.2.5实验体会与总结 (36)3.3实验三 (36)3.3.1实验题目 (36)3.3.2实验目的 (36)3.3.3建模概述 (37)3.3.4计算结果分析与结论 (37)3.3.5实验体会与总结 (44)1 等参单元及其应用1.1 概述1.1.1 等参单元的概念、原理普通单元受到两个方面的限制：（1）单元的精度。

单元的节点数越多，单元精度越高；（2）单元几何上的限制。

普通矩形和六面体单元都不能模拟任意形状几何体，所有几种普通单元都是直线边界，处理曲边界几何体误差较大。

为了解决上述矛盾，方法就是突破矩形单元和六面体单元几何方面的限制，使其成为任意四边形和任意六面体单元，这类单元位移模式和形函数的构造和单元列式的导出不能沿用构造简单单元的方法，必须引入等参变换，采用相同的插值函数对单元的节点坐标和节点位移在单元上进行插值。

有限元方法基础教程第三版答案第二单元1. 诉述有限元法的定义答:有限元法是近似求解一般连续场问题的数值方法2. 有限元法的基本思想是什么答:首先,将表示结构的连续离散为若干个子域,单元之间通过其边界上的节点连接成组合体。

其次,用每个单元内所假设的近似函数分片地表示求解域内待求的未知厂变量。

3. 有限元法的分类和基本步骤有哪些答:分类:位移法、力法、混合法;步骤:结构的离散化,单元分析,单元集成,引入约束条件,求解线性方程组,得出节点位移。

4. 有限元法有哪些优缺点答:优点:有限元法可以模拟各种几何形状复杂的结构,得出其近似解;通过计算机程序,可以广泛地应用于各种场合;可以从其他CAD软件中导入建好的模型;数学处理比较方便,对复杂形状的结构也能适用;有限元法和优化设计方法相结合,以便发挥各自的优点。

缺点:有限元计算,尤其是复杂问题的分析计算,所耗费的计算时间、内存和磁盘空间等计算资源是相当惊人的。

对无限求解域问题没有较好的处理办法。

尽管现有的有限元软件多数使用了网络自适应技术,但在具体应用时,采用什么类型的单元、多大的网络密度等都要完全依赖适用者的经验。

5. 梁单元和平面钢架结构单元的自由度由什么确定答:由每个节点位移分量的总和确定6. 简述单元刚度矩阵的性质和矩阵元素的物理意义答:单元刚度矩阵是描述单元节点力和节点位移之间关系的矩阵单元刚度矩阵中元素aml的物理意义为单元第L个节点位移分量等于1,其他节点位移分量等于0时,对应的第m个节点力分量。

7. 有限元法基本方程中的每一项的意义是什么P14 答:Q——整个结构的节点载荷列阵(外载荷、约束力);整个结构的节点位移列阵;结构的整体刚度矩阵,又称总刚度矩阵。

8. 位移边界条件和载荷边界条件的意义是什么答:由于刚度矩阵的线性相关性不能得到解,引入边界条件,使整体刚度矩阵求的唯一解。

9. 简述整体刚度矩阵的性质和特点P14 答:对称性;奇异性;稀疏性;对角线上的元素恒为正。