全矩阵无积分间断有限元方法实现
非结构网格下2D Riesz分数阶方程的Galerkin有限元方法
非结构网格下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算法原理
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.进行优化和设计：利用有限元模拟的结果，可以进行系统的优化和设计改进。
这可以通过改变几何形状、材料属性或边界条件来实现。
梁杆结构几何非线性有限元的数值实现方法
NUMERICAL IMPLEMENTATION OF GEOMETRICALLY NONLINEAR FINITE ELEMENT METHOD FOR BEAM STRUCTURES
CHEN Zheng-qing
(College of Civil Engineering, Hunan University, Changsha 410082, China)
= tσ ij + ∆∗T ij = ∆∗ Eij
(1) (2)
而它在 t+Δt 时刻柯西应变就等于其增量：
t + ∆t t Eij
式中， ∆ Eij 为：
∗
∆∗ Eij = ∆∗ε ij + ∆∗ηij 1 ∆∗ε ij = (∆ui ,j + ∆u j ,i ) 2 1 ∆∗ηij = ∆uk ,i ∆uk ,j 2
———————————————
收稿日期：2013-05-01；修改日期：2014-03-06 基金项目：国家自然科学基金项目(91215302) 作者简介： 陈政清(1947―)， 男， 湖南湘潭人， 教授， 博士， 湖南大学风工程研究中心主任， 主要从事结构振动与控制研究(E-mail: zqchen@).
(3) (4) (5)
44
工
程
力
学
E G [ t kαβ ]{∆qα } = {t+ ∆t Pβ − tψ β } + t kαβ
仍然假定变形体的应变增量是小应变，应 力应变增量关系可以记为：
(14) (15) (16)
′ ∆∗ε kl ∆∗T ij = Cijki
功增量方程如下： ′ = A3 ′ − A4 ′ A1′ + A2 式中：
有限元分析考试
有限元复习资料工程问题的研究对象：离散体：材料力学中的连续梁、建筑结构框架和桁架结构。
我们把这类问题，称为离散系统。
离散问题是可解的连续体：通常可以建立它们应遵循的基本方程，即微分方程和相应的边界条件。
通常只能得到少数简单问题的精确解答。
对于许多实际的工程问题，还无法给出精确的解答。
有限元方法的两大应用:科学计算或数值模拟数字设计或虚拟仿真(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编程实现有限元切线刚度矩阵计算一、引言有限元法是一种用于求解复杂工程问题的数值分析方法，它将连续介质划分为许多小的单元，通过对每个单元进行离散化处理，再用数学方法对这些单元进行组装，最终得到整个结构的解。
在有限元方法中，刚度矩阵是求解结构问题的关键步骤之一，而有限元切线刚度矩阵的计算则是其中的重要内容之一。
本文将介绍如何使用MATLAB编程实现有限元切线刚度矩阵的计算。
二、有限元切线刚度矩阵的基本概念1. 切线刚度矩阵在有限元方法中，切线刚度矩阵是描述结构对外部载荷作用下的应变-应力关系的重要矩阵。
它描述了结构在外部载荷下的变形行为，是求解结构变形和应力的重要工具。
2. 切线刚度矩阵的计算切线刚度矩阵的计算是通过对单元的局部坐标系进行刚度矩阵的求解，并进行坐标变换得到全局坐标系下的切线刚度矩阵。
在实际计算中，需要考虑单元的几何形状、材料性质等因素，以及在单元上施加的外部载荷。
三、MATLAB编程实现有限元切线刚度矩阵的基本步骤1. 单元刚度矩阵的计算我们需要编写MATLAB函数来实现对单元刚度矩阵的计算。
这个函数需要考虑单元的几何形状、材料性质等因素，以及在单元上施加的外部载荷。
通常情况下，我们可以利用数值积分的方法来进行刚度矩阵的计算。
2. 坐标变换矩阵的计算在得到单元刚度矩阵之后，我们需要计算坐标变换矩阵，将单元刚度矩阵从局部坐标系变换到全局坐标系。
这也需要编写一个MATLAB函数来实现坐标变换矩阵的计算。
3. 矩阵组装我们需要将所有单元的切线刚度矩阵组装成整体刚度矩阵。
这通常需要考虑到单元之间的连接关系，以及边界条件等因素。
在MATLAB中，我们可以利用矩阵的组合和相加等运算来实现整体刚度矩阵的计算。
四、编程实例这里我们以一个简单的弹簧-弹簧系统为例，介绍如何使用MATLAB编程实现有限元切线刚度矩阵的计算。
我们需要定义系统的几何形状、材料性质等参数，然后编写MATLAB函数来进行刚度矩阵的计算，坐标变换矩阵的计算，以及矩阵的组装，最终得到整体刚度矩阵，并求解系统的变形和应力。