四节点等参单元xfem程序设计
四节点四边形单元悬臂梁的matlab有限元编程-概述说明以及解释
四节点四边形单元悬臂梁的matlab有限元编程-概述说明以及解释1.引言1.1 概述悬臂梁是一种常见的结构形式，在工程领域中被广泛应用。
四节点四边形单元是有限元分析中常用的元素类型，能够准确地模拟悬臂梁的受力情况。
Matlab是一种强大的数学工具，可以用来编程实现有限元分析。
本文旨在介绍如何利用Matlab进行四节点四边形单元悬臂梁的有限元编程，并对其进行分析和展望。
通过本文的研究，我们希望能够为工程实践提供一定的参考和指导，同时也为进一步的研究提供基础。
1.2 文章结构本文主要分为三个部分：引言、正文和结论。
引言部分将介绍文章的背景和目的，明确文章研究的问题和意义。
正文部分包括理论基础、Matlab有限元编程介绍和四节点四边形单元悬臂梁建模三个小节。
其中，理论基础将介绍与悬臂梁有关的理论知识，Matlab有限元编程介绍将详细介绍如何使用Matlab进行有限元分析，最后，四节点四边形单元悬臂梁建模将展示具体的悬臂梁建模过程。
结论部分将对实验结果进行分析与总结，探讨本研究的意义和潜在研究方向。
1.3 目的本文旨在利用Matlab编程实现四节点四边形单元悬臂梁的有限元分析，通过建立合适的数学模型，探索悬臂梁在受力状态下的力学特性。
具体目的包括：1. 建立悬臂梁的有限元数学模型，包括节点、单元和材料参数的设置；2. 实现悬臂梁在不同受力情况下的应力、应变、位移等力学性能的计算；3. 分析悬臂梁受力情况下的应力分布情况，探讨悬臂梁的破坏模式和极限承载能力；4. 验证Matlab编程方法的有效性和准确性，为工程实际中悬臂梁等复杂结构的有限元分析提供参考和借鉴。
通过本文的研究，旨在为工程实践提供可靠的数值计算工具和理论分析方法，为解决工程结构强度和稳定性问题提供一定的指导和参考价值。
2.正文2.1 理论基础在介绍四节点四边形单元悬臂梁的Matlab有限元编程之前，我们首先需要了解一些基本的理论知识。
悬臂梁是一种常见的结构形式，在工程领域中广泛应用于桥梁、建筑物等领域。
平面四节点等参单元分析程序
平面四节点等参单元分析程序(总13页)--本页仅作为文档封面，使用时请直接删除即可----内页可以根据需求调整合适字体及大小--变分原理与有限元大作业平面四节点等参单元分析程序姓名：潘清学号：SQ完成时间：2011-4-26一、概述通常情况下的有限元分析过程是运用可视化分析软件(如ANSYS、ABAQUS、SAP等)进行前处理和后处理，而中间的计算部分一般采用自己编制的程序来运算。
具有较强数值计算和处理能力的Fortran语言是传统有限元计算的首选语言。
随着有限元技术的逐步成熟,它被应用在越来越复杂的问题处理中，但在实际应用中也暴露出一些问题。
有时网格离散化的区域较大,而又限于研究精度的要求，使得划分的网格数目极其庞大，结点数可多达数万个，从而造成计算中要运算的数据量巨大，程序运行的时间较长的弊端，这就延长了问题解决的时间，使得求解效率降低。
因为运行周期长，不利于程序的调试，特别是对于要计算多种运行工况时的情况；同时大数据量处理对计算机的内存和CPU 提出了更高的要求，而在实际应用中，单靠计算机硬件水平的提高来解决问题的能力是有限的。
因此，必须寻找新的编程语言。
随着有限元前后处理的不断发展和完善，以及大型工程分析软件对有限元接口的要求，有限元分析程序不应只满足解题功能，它还应满足软件工程所要求的结构化程序设计条件，能够对存储进行动态分配，以充分利用计算机资源，它还应很容易地与其它软件如CAD 的实体造型，优化设计等接口。
现在可编写工程应用软件的计算机语言较多，其中C语言是一个较为优秀的语言，很容易满足现在有限元分析程序编程的要求。
C语言最初是为操作系统、编译器以及文字处理等编程而发明的。
随着不断完善，它已应用到其它领域，包括工程应用软件的编程。
近年来，C语言已经成为计算机领域最普及的一个编程语言，几乎世界上所有的计算机都装有C 的编译器，从PC机到巨型机到超巨型的并行机，C与所有的硬件和操作系统联系在一起。
双线性四边形等参单元有限元程序
(9)
(10)
(11)
进一步，我们有：
⎛ ∂u ⎞ ⎜ ∂ξ ⎟ ⎜ ⎟ 1 −η ⎡η − 1 0 ⎜ ∂u ⎟ ⎜ ∂η ⎟ 1 ⎢ξ − 1 0 −1 − ξ ⎜ ⎟= ⎢ 0 ⎜ ∂v ⎟ 4 ⎢ 0 η − 1 ⎢ ⎜ ∂ξ ⎟ 0 ⎣ 0 ξ −1 ⎜ ⎟ ⎜ ∂v ⎟ ⎜η ⎟ ⎝ ⎠
4.2 单元刚度矩阵及整体刚度矩阵的集成
(7)
4 节点 4 边形单元刚度矩阵从单元应变能而来：
1 U e = h ∫ [σ ]T [ε ]dA 2 e (8)
h 为结构厚度。 应变矩阵可以表示为位移的函数：
⎛ ∂u ⎞ ⎜ ⎟ ⎛ ε x ⎞ ⎜ ∂x ⎟ ⎜ ⎟ ⎜ ∂v ⎟ [ε ] = ⎜ ε y ⎟ = ⎜ ⎟ ⎜ ⎟ ⎜ ∂y ⎟ ⎝ γ xy ⎠ ⎜ ∂u ∂v ⎟ ⎜ + ⎟ ⎝ ∂y ∂x ⎠ 采用前面的微分转化公式(7)将[ε]表示为局部坐标的微分形式： ⎛ ∂u ⎞ ⎜ ∂ξ ⎟ ⎜ ⎟ ⎜ ∂u ⎟ ⎜ ∂η ⎟ ⎟ [ε ] = [ A] ⎜ ⎜ ∂v ⎟ ⎜ ∂ξ ⎟ ⎜ ⎟ ⎜ ∂v ⎟ ⎜η ⎟ ⎝ ⎠ 这里：
一个典型的计算结果如图 6 所示：
图 6 输出结果
其中 u 为节点位移，u(1,i)为第 i 个节点 x 方向的位移，u(2,i)为第 i 个节点 y 方向的位移。strNote 为节点应力，strNote(1,i)为第 i 个节点的 σx，strNote(2,i) 为 第 i 个节点的 σy， strNote(1,i)为第 i 个节点的 σxy。 StrCent 为单元最佳应力点应力， strCent(1,i)为第 i 个单元的 σx，strCent(1,i)为第 i 个单元的 σy，strCent(1,i)为第 i 个单元的 σxy。 2 程序结构说明 本程序的结构和一般的有限元程序没什么不同，只是没有必要的前处理过 程，需要用户手动输入。程序流程图如图 6 所示： 主程序 main()调用过程如下： 输入数据加载：
abaqus中xfem扩展有限元教程
abaqus中xfem扩展有限元教程一、part模块中的操作：1. 生成一个新的part，取名为plate，本part选取3D deformable solid extrusion类型（如图1）2.通过Rectangle工具画出一长3，高6的矩形。
考虑使用工具栏add-dimension和edit dimension来画出精确长度的模型。
强烈建议此矩形的左上角坐标为（0，3），右下角坐标为（3，-3）（如图2）3. 完成后拉伸此矩形，深度为1.（如图3）4. 生成一个新的part，取名为crack，本part选取3D deformable shell extrusion类型（如图4）5. 生成一条线，此线的左端点坐标为（0，0.08），右端点坐标为（1.5，0.08）6 . 完成后拉伸此线，深度为1.（如图6）7.保存此模型为XFEMtutor（如图7），以后经常保存模型，不再累述。
8. 在part Plate中分别创建4个集合，分别为：all，bottom，top和fixZ，各部分的内容如图8～11所示。
二、Material模块中的操作：1. 创建材料elsa，其弹性参数为E=210GPa，泊松比为0.3（如图12）最大主应力失效准则作为损伤起始的判据，最大主应力为84.4MPa（如图13）损伤演化选取基于能量的、线性软化的、混合模式的指数损伤演化规律，有关参数为G1C= G2C= G3C=42200N/m， =1.（如图14）2.创建一个Solid Homogeneous 的section，名为solid（如图15），此section与材料elsa相联（如图16），并将此section赋给plate part（也就是集合all）（如图17）3.赋予材料取向，分别如图18～21所示。
三、划分网格：网格控制为：Hex型structured（如图22），单元类型为C3D8R（如图23）设置plate各边的网格种子为8，26，36（如图24），各边种子的个数不能改变（如图25）四、装配模块：选中plate和crack两个part，分别生成2个实体（如图26），生成一个参考点，参考点的坐标为（1.5，-3，0）（如图27，28）。
有限元讲义弹性力学平面问题有限单元法(四结点四边形等参元,八结点曲线四边形等参元,问题补充)分析
2.6 四结点四边形单元(The four-node quadrilateral element)前面介绍了四结点的矩形单元其位移函数:xy y x U 4321αααα+++=xy y x V8765αααα+++=为双线性函数，应力，应变在单元内呈线性变化，比常应力三角形单元精度高。
但它对边界要求严格。
本节介绍的四结点四边形等参元，它不但具有较高的精度，而且其网格划分也不受边界的影响。
对任意四边形单元（图见下面）若仍直接采用前面矩形单元的位移函数，在边界上它便不再是线性的（因边界不与x,y 轴一致），这样会使得相邻两单元在公共边界上的位移可能会出现不连续现象（非协调元），而使收敛性受到影响。
可以验证，利用坐标变换就能解决这个问题，即可以通过坐标变换将整体坐标中的四边形（图a ）变换成在局部坐标系中与四边形方向无关的边长为2的正方形。
正方形四个结点i,j,m,p 按反时钟顺序对应四边形的四个结点i j m p 。
正方形的 1-=η 和 1=η 二条边界，分别对应四边形的i ，j 边界和p,m 边界；ξ=-1和ξ=+1分别对应四边形的i ，p 边界和j ，m 边界。
如果用二组直线等分四边形的四个边界线段，使四边形绘成一个非正交网格，那么该非正交网格在正方形上对应着一个等距离的规则网格（见图a, b ）。
当然, 局部坐标上的A 点与整体坐标的A 点对应。
一、四结点四边形等参单元的形函数及坐标变换由于可以将整体坐标下的四边形单元变换成局部坐标下的正方形单元，对于这种正方形单元，自然仍取形函数为： ξηαηαξαα2321+++=U ξηαηαξαα8765+++=V引入边界条件，即可得位移函数：∑=ijmpi i U N Ui ijmpi V N V ∑==写成矩阵形式：{}{}[]{}ee p i p i ed N d N N N N V U f =⎥⎦⎤⎢⎣⎡=⎭⎬⎫⎩⎨⎧=000 式中形函数: ()()()ηηξξηξi i i N ++=1141, ()p m j i ,,, 按照等参元的定义,我们将坐标变换式亦取为: p p m m j j i i i ijmpi x N x N x N x N x N x +++==∑p p m m j j i i i ijmpi y N y N y N y N y N y +++==∑ ()162-- 式中形函数N 与位移函数中的完全一致。
c++ 二维四节点四边形等参元刚度矩阵
C++ 二维四节点四边形等参元刚度矩阵在计算机辅助工程领域，C++语言被广泛应用于有限元分析（FEA）和计算力学等方面。
而在有限元分析中，等参元（Isoparametric Element）是一种常用的元素类型，用于对复杂的结构和材料进行建模和分析。
本文将围绕C++语言下的二维四节点四边形等参元刚度矩阵展开讨论。
1. 了解四节点四边形等参元在开始讨论C++下的四节点四边形等参元刚度矩阵之前，首先应该对四节点四边形等参元有所了解。
四节点四边形等参元是指在二维空间中，以四个节点构成的四边形元素，同时该元素的几何形状和内部分布的自由度均由相同的形函数进行描述的元素。
在有限元分析中，等参元的应用可以大大简化模型建立的复杂度，并提高计算的精度。
2. 实现C++中的四节点四边形等参元在C++语言中，实现四节点四边形等参元需要考虑节点的坐标、材料属性、边界条件等因素，并结合有限元理论中的导出公式进行编码。
通过C++语言的面向对象特性，可以将节点、单元、材料等抽象为对象，以便更好地管理和组织相关数据和函数。
3. 深入分析四节点四边形等参元的刚度矩阵四节点四边形等参元的刚度矩阵是描述该元素对外加载的响应性能，是有限元分析中至关重要的一部分。
通过推导理论公式，并结合C++语言进行具体实现，我们可以得到该元素的刚度矩阵。
在这一过程中，需要考虑矩阵装配的方法、高效的数据结构、数值计算的稳定性等因素。
4. 个人观点和总结通过对C++下的二维四节点四边形等参元刚度矩阵的探讨，我深刻地意识到了有限元分析与计算力学领域的复杂性和重要性。
C++语言作为一种高效、灵活的编程语言，为工程领域的数值计算提供了强大的支持。
对于工程师和研究人员来说，深入理解和掌握有限元分析的原理和实现方法，将有助于提升对复杂结构和材料行为的认识和预测能力。
对于C++语言下的四节点四边形等参元刚度矩阵，需要我们不仅要掌握有限元分析的理论知识，还需要熟练掌握C++语言的编程技巧和工程应用。
第九章 平面广义四节点等参元 GQ4
18.24 23.00 23.24 23.06 23.43 23.78 23.02 23.04 23.69 23.67 23.82 23.9
0.2113 0.2209 0.2287 0.2221 0.23337 0.2226 0.222 0.2661 0.2261 0.2393 0.2377 0.2360
式中 N e 为单元 e 的插值函数矩阵 值函数对应的广义形函数形式为
De 为单元 e 的广义自由度向量 一阶广义节点插 y , (i = 1,2,3,4) 0 x2 0 y2 0
1 0 x 0 y N ei = ϕ ei 0 1 0 x 0 1 0 x 0 y N =ϕ 0 1 0 x 0
9.3.3
剪切自锁考查
MacNeal 细长梁问题
弹性模量为 E = 106 泊松比为
计算简图见图 9-3
材料参数选为
ν = 0.3 纯
弯和端部受剪两种工况 计算网格及工况如图 9-3 中(a) 格计算结果列于表 9-3 中
(b)和(c)
不同工况下各种网
0.5 50 0.2 1 工况 1 0.5 50 工况 2
( ) (E )(B )
i T e g
(9-10)
j e g
式中 n g 为单元内高斯点个数 取值
t 表示材料的厚度
下标 g 表示该表达式在高斯点处的
W g 为高斯点积分权系数 具体的数值实现步骤如下
(1) 按照所选用的高斯积分阶次 (2) 按单元节点循环 i. 形成该单元的所有高斯点局部坐标 (ξ i ,η i )
vC
网格 a
ε A (max) ε B (min)
vC
网格 b
ε A (max) ε B (min)
基于MATLAB的四结点等参元的编程
基于MATLAB的四结点等参元的编程
栾小兵;李自林;朱斌
【期刊名称】《天津城市建设学院学报》
【年(卷),期】2006(012)004
【摘要】基于有限元理论,用MATLAB语言编制了四结点等参元计算机程序.计算了悬臂双向板在均布荷载与集中力共同作用下的位移、应力和应变,并与大型通用有限元软件ANSYS计算结果进行了比较,二者吻合较好,验证了本程序编制的正确性.利用MATLAB编程对结构进行分析计算比FORTRAN、C语言简便,具有精度高的优点.
【总页数】4页(P259-262)
【作者】栾小兵;李自林;朱斌
【作者单位】天津城市建设学院,土木工程系,天津,300384;天津城市建设学院,土木工程系,天津,300384;天津城市建设学院,土木工程系,天津,300384
【正文语种】中文
【中图分类】TP314
【相关文献】
1.四结点等参元XFEM程序设计及在裂纹问题中的应用 [J], 苏毅;王生楠;闫晓中
2.机械密封环应力场的四结点等参元分析 [J], 皇甫晶;周启慧;魏昭成
3.简明有效的四结点多变量协调等参元 [J], 史宝军;鹿晓阳
4.平面四结点多变量非协调等参元 [J], 史宝军;袁明武;等
5.基于“求逆犯规”修正方法的平面四结点非协调等参元 [J], 鹿晓阳
因版权原因，仅展示原文概要，查看原文内容请购买。
四边形四节点等参元matlab程序
悬臂钢梁，尺寸如图一所示；v=0.3。
h=1，E=2.1e11.图一悬臂钢梁图二单元划分与结点编号附录：%---------------四边形四节点等参元matlab计算程序----------------------------% 2013年% 13级建筑与土木工程Bruce% B15-405% 本程序只输出了节点位移、单元中心点的应力%******************************************************************* % 变量说明% E v h% 弹性模量泊松比厚度% NPOIN NELEM NVFIX NNODE NFPOIN% 总结点数, 单元数, 约束结点个数, 单元节点数,受力结点数% COORD LNODS% 结构节点整体坐标数组, 单元定义数组,% FPOIN FORCE FIXED% 结点力数组，总体荷载向量, 约束信息数组% HK DISP% 总体刚度矩阵,结点位移向量%******************************clear allformat short eFP1=fopen(' sjd.txt','rt'); %打开数据文件%%读入控制数据E=fscanf(FP1,'%f',1); %弹性模量v=fscanf(FP1,'%f',1); % 泊松比h=fscanf(FP1,'%f',1); %厚度NELEM=fscanf(FP1,'%d',1); %单元数NPOIN=fscanf(FP1,'%d',1); % 总结点数NNODE=fscanf(FP1,'%d',1); %单元节点数NFPOIN=fscanf(FP1,'%d',1); %受力结点数NVFIX=fscanf(FP1,'%d',1); %约束结点个数LNODS=fscanf(FP1,'%f',[NNODE,NELEM])'; % 单元定义：单元结点号（逆时针）COORD=fscanf(FP1,'%f',[2,NPOIN])'; % 结点号x,y坐标(整体坐标下) FPOIN=fscanf(FP1,'%f',[3,NFPOIN])';% 节点力：结点号、X方向力(向右正)，Y方向力(向上正)FIXED=fscanf(FP1,'%d',[3,NVFIX])';%约束信息数组(n,3) n:受约束节点数目, (n,1):约束点号%(n,2)与(n,3)分别为约束点x方向和y方向的约束情况,受约束为1否则为0% 刚度矩阵的生成%计算刚度矩阵，并对约束条件进行处理Ke=zeros(2*NNODE,2*NNODE); % 单元刚度矩阵并清零HK=zeros(2*NPOIN,2*NPOIN); % 张成总刚矩阵并清零%调用子程序生成单元刚度矩阵for m=1:NELEM %m为单元号Ke=K(E,v,h,...COORD(LNODS(m,1),1),COORD(LNODS(m,1),2),...COORD(LNODS(m,2),1),COORD(LNODS(m,2),2),...COORD(LNODS(m,3),1),COORD(LNODS(m,3),2),...COORD(LNODS(m,4),1),COORD(LNODS(m,4),2)); %调用单元刚度矩阵a=LNODS(m,:); %临时向量,用来记录当前单元的节点编号%对总刚度矩阵的处理for j=1:4for k=1:4HK((a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)=HK((a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)+...Ke(j*2-1:j*2,k*2-1:k*2);endendend %—————————————————————————————————%对荷载向量进行处理FORCE=zeros(2*NPOIN,1); % 张成总荷载向量并清零for i=1:NFPOINb1=FPOIN(i,1)*2-1;b2=FPOIN(i,1)*2; %FPION(i,1)为作用点FORCE(b1)=FPOIN(i,2); %FPION(i,2)为x方向的节点力FORCE(b2)=FPOIN(i,3); %FPION(i,3)为y方向的节点力end %—————————————————————————————————%将约束信息加入总刚,总荷载for i=1:NVFIXif FIXED(i,2)==1c1=2*FIXED(i,1)-1;HK(c1,:)=0; %将一约束序号处的总刚列向量清0HK(:,c1)=0; %将一约束序号处的总刚行向量清0HK(c1,c1)=1; %将行列交叉处的元素置为1FORCE(c1)=0;endif FIXED(i,3)==1c2=2*FIXED(i,1);HK(c2,:)=0;HK(:,c2)=0;HK(c2,c2)=1;FORCE(c2)=0;endendDISP=HK\FORCE %计算节点位移向量%———————————求解单元应力————————————————stress=zeros(3,NELEM);for m=1:NELEMu(1:8)=0;d=LNODS(m,:); %临时向量,用来记录当前单元的节点编号for i=1:NNODEu(i*2-1:i*2)=DISP(d(i)*2-1:d(i)*2);%从总位移向量中取出当前单元的节点位移endD=(E/(1-v*v))*[1 v 0;v 1 0;0 0 (1-v)/2];%弹性矩阵%形成应变矩阵BMBM=zeros(3,8);for i=1:NNODEJ=Jacobi(COORD(LNODS(m,1),1),COORD(LNODS(m,1),2),...COORD(LNODS(m,2),1),COORD(LNODS(m,2),2),...COORD(LNODS(m,3),1),COORD(LNODS(m,3),2),...COORD(LNODS(m,4),1),COORD(LNODS(m,4),2),0,0);[N_s,N_t]=DHS(0,0);B1i=J(2,2)*N_s(i)-J(1,2)*N_t(i);B2i=-J(2,1)*N_s(i)+J(1,1)*N_t(i);BM(1:3,2*i-1:2*i)=[B1i 0;0 B2i;B2i B1i]/det(J);endstressm=D*BM*u';stress(:,m)=stressm;endstress %输出应力function Ke=K(E,v,h,x1,y1,x2,y2,x3,y3,x4,y4)%=========单元刚度矩阵=============== % E 弹性模量% v 泊松比% h 厚度% x1,y1,x2,y2,x3,y3,x4,y4 为4个角结点的坐标%矩阵尺寸为8 x 8Ke=zeros(8,8);D=(E/(1-v*v))*[1 v 0;v 1 0;0 0 (1-v)/2];%弹性矩阵%高斯积分采用3 x 3 个积分点W1=5/9;W2=8/9;W3=5/9; %加权系数W=[W1 W2 W3];r=15^(1/2)/5;x=[-r 0 r];%积分点for i=1:3for j=1:3B=eleB(x1,y1,x2,y2,x3,y3,x4,y4,x(i),x(j));J=Jacobi(x1,y1,x2,y2,x3,y3,x4,y4,x(i),x(j));Ke=Ke+W(i)*W(j)*B'*D*B*det(J)*h;endendfunction B=eleB(x1,y1,x2,y2,x3,y3,x4,y4,s,t)%调用导函数[N_s,N_t]=DHS(s,t);%求Jacobi矩阵J=Jacobi(x1,y1,x2,y2,x3,y3,x4,y4,s,t);%求应变矩阵BB=zeros(3,8);for i=1:4B1i=J(2,2)*N_s(i)-J(1,2)*N_t(i);B2i=-J(2,1)*N_s(i)+J(1,1)*N_t(i);B(1:3,2*i-1:2*i)=[B1i 0;0 B2i;B2i B1i];endB=B/det(J);function J=Jacobi(x1,y1,x2,y2,x3,y3,x4,y4,s,t) %-------Jacobi-----------%单元坐标x=[x1 x2 x3 x4];y=[y1 y2 y3 y4];%%调用形函数对局部坐标的导数[N_s,N_t]=DHS(s,t);%求Jacobi矩阵的行列式的值x_s=0;y_s=0;x_t=0;y_t=0;for i=1:4x_s=x_s+N_s(i)*x(i);y_s=y_s+N_s(i)*y(i);x_t=x_t+N_t(i)*x(i);y_t=y_t+N_t(i)*y(i);endJ=[x_s y_s;x_t y_t];function N=shape(s,t)%ξ,ηN(1)=(1-s)*(1-t)/4;N(2)=(1+s)*(1-t)/4;N(3)=(1+s)*(1+t)/4;N(4)=(1-s)*(1+t)/4;function [N_s,N_t]=DHS(s,t)%形函数求导%ξ,ηN_s(1)=-1/4*(1-t);N_s(2)=1/4*(1-t);N_s(3)=1/4*(1+t);N_s(4)=-1/4*(1+t);N_t(1)=-1/4*(1-s);N_t(2)=-1/4*(1+s);N_t(3)=1/4*(1+s);N_t(4)=1/4*(1-s);sjd.txt 文件数据2.1E11 0.3 1 5 12 4 1 21 2 8 72 3 9 83 4 10 94 5 11 105 6 12 110.0 0.01.0 0.02.0 0.03.0 0.04.0 0.05.0 0.00.0 1.01.0 1.02.0 1.03.0 1.04.0 1.05.0 1.06 0 -100001 1 17 1 1。
等参单元
等参变换
例二：考察实际单元为矩形单元(2a*2b)情形，定 义如前 （3.10） 式中，（x0,y0)为局部坐标原点，由上式得 （3.11） 重新组合 有 （3.12）
等参变换
将式（3.12）与（3.6）对照，即有 同理可得
这说明矩形单元只是四节点四边形单元的一个特 例，通过这个特例可以加深对实际单元到基本单 元变换的理解。
Ni , x E S D B Ni , x i i 2 1 1 1 Ni , y 2 Ni , y Ni , y 1 1 Ni , x 2
四节点四边形下等参元矩阵
下面求解单元刚矩[K] 目标：四节点四边形单元刚度矩阵[K] 问题：将[K]中的面积微元dxdy用dξ dη 进行代换
（3.29） 其中[K]e可划分成四行四列的子矩阵
K
e
B D B tdxdy
T
K
e
k11 k 21 k31 k41
k12 k22 k32 k42
k13 k23 k33 k43
k14 k24 k34 k44
写成矩阵形式 Ni , x, Ni , x, 其中 x x x, x,
y, Ni , x N y, i, y
y y,
（3.15）
y y,
四节点四边形下等参元矩阵
（3.9）
等参变换
把以参数 ξ 代表的 x方程和 y方程消去 ξ ， 则得 x , y 所组成的直线方程 y=kx+b 这表明ξ η 平面上的水平边界经过变换成为了xy平 面上的斜边界，这也是双线性函数的特点体现,即所 以母单元上的四个边都可以通过映射在x y坐标面上 得出一个任意四边形。
最新平面四边形4结点等参有限单元法
有限元程序设计平面四边形4结点等参有限单元法程序设计1、程序功能及特点a.该程序采用四边形4节点等参单元，能解决弹性力学的平面应力应变问题。
b.前处理采用网格自动划分技术，自动生成单元及结点信息。
b.能计算受集中力、自重体力、分布面力和静水压力的作用。
c.计算结点的位移和单元中心点的应力分量及其主应力。
d.后处理采取整体应力磨平求得各个结点的应力分量。
e.算例计算结果与ANSYS计算结果比较，并给出误差分析。
f.程序采用Visual Fortran 5.0编制而成。
2、程序流程及图框图2-１程序流程图图2-２子程序框图其中，各子程序的主要功能为：INPUT――输入原始数据HUAFEN――自动网格划分，形成COOR(2,NP),X,Y的坐标值与单元信息CBAND――形成主元素序号指示矩阵MA（*）SKO――形成整体刚度矩阵[K]CONCR――计算集中力引起的等效结点荷载{R}eBODYR――计算自重体力引起的等效结点荷载{R}eFACER――计算分布面力引起的等效结点荷载{R}eDECOP――支配方程LU三角分解FOBA――LU分解直接解法中的回代过程OUTDISP――输出结点位移分量STRESS――计算单元应力分量OUTSTRE――输出单元应力分量STIF――计算单元刚度矩阵FDNX――计算形函数对整体坐标的导数TiiyNxN⎥⎦⎤⎢⎣⎡∂∂∂∂，=i1,2,3,4。
FUN8――计算形函数及雅可比矩阵[J]SFUN ――应力磨平-单元下的‘K’=NCN‘SCN――应力磨平-单元下的右端项系数‘CN‘SUMSKN――应力磨平-单元下的右端项集成到总体的‘P‘SUMSTRS――应力磨平-单元下的集成到总体的‘K‘GAUSTRSS――高斯消元求磨平后的应力3、输入数据及变量说明当程序开始运行时，按屏幕提示，键入数据文件的名字。
在运行程序之前，根据程序中INPUT需要的数据输入建立一个存放原始数据的文件，这个文件的名字为INDAT.DAT。
平面四节点等参单元matlab实现
计算力学报告平面四节点等参单元学生姓名：**学号：********一、问题描述及分析在无限大平面内有一个小圆孔。
孔内有一集中力p，试求用有限元法编程和用ANSYS软件求出各点应力分量和位移分量，并比较二者结果。
根据圣维南原理建立半径为10mm的大圆，设小圆孔的半径a=0.5mm，在远离大圆边界的地方模型是比较精确的。
由于作用在小圆孔上的力引起的位移随距离的衰减非常快，所以可以把大圆边界条件设为位移为零。
二、有限元划分描述在划分单元时，单元数量比较多，于是我采取了使用ansys软件建模自动划分单元网格的方法。
具体操作如下：打开ansys，在单元类型中选择solid->Quad 4 node 182单元；建立类半径为0.5外半径为10的圆环；使用mashtool中的智能划分和将单元退化成三角形单元；使用工具栏中List中的Nodes和Elements 选项将节点和单元数据导出并导入Excle中，总共得到了207个单元和229个节点。
如下图：图1三、有限元程序及求解程序求解使用了matlab语言。
具体如下：程序：clcclearE=2e11; %弹性模量NU=0.3; %泊松比t=0.1; %厚度X=xlsread('D:\data','nodes'); %读取节点坐标elem=xlsread('D:\data','elements'); %读取单元编号w=[1,2,3,4,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28]; %有位移约束的节点n=size(X,1); %节点数m=size(elem,1); %单元数K=zeros(2*n); %初始总体刚度矩阵for i=1:msyms Ks Et x y I1 I2 a b B; %定义可能存在的变量e=[1,1;-1,1;-1,-1;1,-1];for j=1:4 %形函数N(j)=0.25*(1+e(j,1)*Ks)*(1+e(j,2)*Et);endx=0;y=0;for j=1:4 %标准母单元映射到真实单元x=x+N(j)*X(elem(i,j),1);y=y+N(j)*X(elem(i,j),2);endJ1=jacobian([x;y],[Ks;Et]); %雅克比矩阵及其转置J=J1';for j=1:4I1=diff(N(j),Ks); %形函数分别对Ks和Et的偏导数I2=diff(N(j),Et);C=(J^-1)*[I1;I2];a=C(1); %形函数对x，y的偏导数b=C(2);B(1,2*j-1)=a; %组成B阵B(1,2*j)=0;B(2,2*j-1)=0;B(2,2*j)=b;B(3,2*j-1)=b;B(3,2*j)=a;endD=(E/(1-NU*NU))*[1,NU,0;NU,1,0;0,0,(1-NU)/2]; %D阵k=zeros(8,8);Kss=[-0.906179,-0.538469,0,0.538469,0.906179]; %5*5高斯积分点Ett=[-0.906179,-0.538469,0,0.538469,0.906179];H=[0.236926,0.478628,0.568888,0.478628,0.236926];%高斯积分权系数for j=1:5 %高斯积分求单元刚度阵for l=1:5Ks=Kss(j);Et=Ett(l);B=subs(B);J=subs(J);k=k+H(j)*H(l)*B'*D*B*det(J);endendG=zeros(8,2*n); %初始总刚变换矩阵G(1,2*elem(i,1)-1)=1; %总刚变换矩阵G(2,2*elem(i,1))=1;G(3,2*elem(i,2)-1)=1;G(4,2*elem(i,2))=1;G(5,2*elem(i,3)-1)=1;G(6,2*elem(i,3))=1;G(7,2*elem(i,4)-1)=1;G(8,2*elem(i,4))=1;K=K+G'*k*G; %总体刚度矩阵合成endKK=K;b=size(w,1);for i=1:bK(2*w(i)-1,2*w(i)-1)=1e20;K(2*w(i),2*w(i))=1e20;end %置大数法f=zeros(2*n,1); %初始载荷矩阵f(10)=-10e3; %加载荷10kNU=K\f; %节点位移for i=1:m %将每个单元各个节点位移集合u(:,i)=[U(2*elem(i,1)-1);U(2*elem(i,1));U(2*elem(i,2)-1);U(2*elem(i,2));U(2*ele m(i,3)-1);U(2*elem(i,3));U(2*elem(i,4)-1);U(2*elem(i,4))];endfor i=1:m %求单元应力syms Ks Et x y I1 I2 a b B;e=[1,1;-1,1;-1,-1;1,-1];for j=1:4N(j)=0.25*(1+e(j,1)*Ks)*(1+e(j,2)*Et);endx=0;y=0;for j=1:4x=x+N(j)*X(elem(i,j),1);y=y+N(j)*X(elem(i,j),2);endJ1=jacobian([x;y],[Ks;Et]);J=J1';for j=1:4I1=diff(N(j),Ks);I2=diff(N(j),Et);C=(J^-1)*[I1;I2];a=C(1);b=C(2);B(1,2*j-1)=a;B(1,2*j)=0;B(2,2*j-1)=0;B(2,2*j)=b;B(3,2*j-1)=b;B(3,2*j)=a; %以上同前面部分为得到B阵endD=(E/(1-NU*NU))*[1,NU,0;NU,1,0;0,0,(1-NU)/2];w=D*B*u(:,i);w1=subs(w,{Ks,Et},{1,1}); %求单元上各节点的应力sigma1(:,i)=double(w1);w2=subs(w,{Ks,Et},{-1,1});sigma2(:,i)=double(w2);w3=subs(w,{Ks,Et},{-1,-1});sigma3(:,i)=double(w3);w4=subs(w,{Ks,Et},{1,-1});sigma4(:,i)=double(w4);endc=[24,29,47,58,78,79,137,149,160,166,186]'; %如截图选取圆半径方向的节点号d=[166,167,168,169,170,171,172,173,174,175,176,177,178,179,180,181,182,183,184 ,185]';%圆周方向选择的节点号e=size(c,1);f=size(d,1);sigmar=zeros(3,e); %r相同角度不同节点应力矩阵sigmat=zeros(3,f); %角度不同r不同节点应力矩阵msigmar=zeros(1,e); %半径方向节点处的mises应力msigmat=zeros(1,f); %圆周方向节点处的mises应力for i=1:e %绕节点平均g=0;for j=1:m %判断节点在单元的那个位置并加上相应的应力值if elem(j,1)-c(i)==0g=g+1;sigmar(:,i)=sigmar(:,i)+sigma1(:,j);endif elem(j,2)-c(i)==0g=g+1;sigmar(:,i)=sigmar(:,i)+sigma2(:,j);endif elem(j,3)-c(i)==0g=g+1;sigmar(:,i)=sigmar(:,i)+sigma3(:,j);endif elem(j,4)-c(i)==0g=g+1;sigmar(:,i)=sigmar(:,i)+sigma4(:,j);endendsigmar(:,i)=sigmar(:,i)/g; %求应力绕节点平均msigmar(:,i)=(0.5*((sigmar(1,i)-sigmar(2,i))^2+sigmar(1,i)^2+sigmar(2,i)^2+6*(si gmar(3,i))^2))^0.5; %求节点处的mises应力endmsigmar %mises应力for i=1:f %同上g=0;for j=1:mif elem(j,1)-d(i)==0g=g+1;sigmat(:,i)=sigmat(:,i)+sigma1(:,j);endif elem(j,2)-d(i)==0g=g+1;sigmat(:,i)=sigmat(:,i)+sigma2(:,j);endif elem(j,3)-d(i)==0g=g+1;sigmat(:,i)=sigmat(:,i)+sigma3(:,j);endif elem(j,4)-d(i)==0g=g+1;sigmat(:,i)=sigmat(:,i)+sigma4(:,j);endendsigmat(:,i)=sigmat(:,i)/g;msigmat(:,i)=(0.5*((sigmat(1,i)-sigmat(2,i))^2+sigmat(1,i)^2+sigmat(2,i)^2+6*(si gmat(3,i))^2))^0.5;endmsigmat四、计算结果：1.ANSYS软件计算结果计算结果分别罗列圆周方向的单元和半径方向的单元位移和应力。
重力坝开裂过程扩展有限元数值模拟
重力坝开裂过程扩展有限元数值模拟靳旭;董羽蕙【摘要】扩展有限元法(XFEM)是一种求解不连续问题的数值方法.它继承了常规有限元法(CFEM)的所有优点,在模拟裂纹扩展、界面、复杂流体等不连续问题时特别有效,近十多年得到了快速发展.介绍了XFEM的基本原理,给出了进行混凝土裂纹扩展分析的方法.利用XFEM模拟混凝土重力坝裂纹扩展,通过对比有、无裂纹情况下的重力坝应力分布,分析裂纹存在对重力坝应力场分布的影响；分析裂纹扩展受网格疏密程度的影响；计算在不同岩基弹性模量下裂纹的扩展方向.%Extended finite element method(XFEM)is a numerical solution for analyzing discontimuity problem . It inherited all the advantages of the conventional finite element method (CFEM) , in the simulation of crack extension , interface, complex fluid and other discontinuities are particularly effective , in the past decade it has been rapid development. The basic theory of XFEM in introduced and the method of analyzing concrete fracture is presented. The XFEM is utilized to simulate the crack propagation in concrete gravity dam. By the contrast of stress distribution under no crack and crack circumstance of gravity dam the discipline of stress field distribution is analyzed; It is also used for influence of mesh density to crack propagation and is calculated the crack propagation direction in batholith elastic modulus.【期刊名称】《科学技术与工程》【年(卷),期】2012(012)033【总页数】6页(P9100-9104,9109)【关键词】重力坝;扩展有限元法;裂纹扩展;网格疏密;弹性模量【作者】靳旭;董羽蕙【作者单位】昆明理工大学建筑工程学院,昆明650500;昆明理工大学建筑工程学院,昆明650500【正文语种】中文【中图分类】TV313;TV642.3实际工程中，无论采用多么严格的裂缝控制措施，混凝土结构仍然会带裂缝工作。
1基于4节点四边形单元的矩形薄板分析(Quad2D4Node).doc
F=o--o-- 000000002 2(4-187)1基于4节点四边形单元的矩形薄板分析(Quad2D4Node)如图4-21所示的一个薄平板，在右端部受集屮力F 作用，其屮的参数为：E = lxlO 7Pa, ;z=l/3,r=0.lm, F=lxl05No 基于MATLAB 平台,按平面应力问题计算各个节点位移、支朋反力以及单元的应力。
(a)问题描述 (b)有限元分析模型图4・21右端部受集中力作用的薄平板解答：对该问题进行有限元分析的过程如下。
(1) 结构的离散化与编号将结构离散为二个4节点矩形单元，单元编号及节点编号如图4-21(b)所示，连接关系 见表4・6,节点的几何坐标见表4・7,载荷F 按静力等效原则向节点1, 2移置。
表4-6结构的单元连接关系单元号节点号I 3 5 6 421 3 4 2表节点的坐标节点节点坐标/mXy 1 2 i 2 2 0 3 1 1 4 10 5 01 6节点位移列阵q = [u A v, u 2 v 2 w 3 v 3 u 4 v 4 u 5 v 5 u 6 v 6](4-186)节点外载列阵P=F+R=R= 0约束的支反力列阵R=[0 0 0 0 0 0 0 0 /?x5 R y5 R X 6 R y6J (4-188)总的节点载荷列阵Rd 心 6 (4 89)心 5其中，(R\5，R 、,5)和(心6,心6)分别为节点5和节点6的两个方向的支反力。
(2) 计算各单元的刚度矩阵(以国际标准单位)首先在MATLAB 环境下，输入弹性模量E 、泊松比NU 、薄板厚度h 和平面应力问题 性质指示参数ID,然后针对单元1和单元2,分别两次调用函数Quad2D4Node_Stiffness, 就可以得到单元的刚度矩阵kl(8X8)和k2(8X8)o» E=le7; » NU=l/3; »1=0」； »ID= 1; » kl = Quad2D4Node_Stiffness(E,NU,t, 1,1,0,1,0,0,1,0, ID); »k2=Quad2D4Node_Stiffness(E,NU,t, 2,1,1,1,1,0,2,0, ID);(3) 建立整体刚度方程由于该结构共有6个节点，则总共的自由度数为12,因此，结构总的刚度矩阵为KK(12 X 12),先对KK 清零，然后两次调用函数Quad2D4Node_Assembly 进行刚度矩阵的组装。
二、平面四边形4结点等参有限单元法程序
-1- -2- 平面四边形4结点等参有限单元法程序设计1、程序功能及特点a.该程序采用四边形4节点等参单元能解决弹性力学的平面应力应变问题。
b.前处理采用网格自动划分技术自动生成单元及结点信息。
b.能计算受集中力、自重体力、分布面力和静水压力的作用。
c.计算结点的位移和单元中心点的应力分量及其主应力。
d.后处理采取整体应力磨平求得各个结点的应力分量。
e.算例计算结果与ANSYS计算结果比较并给出误差分析。
f.程序采用Visual Fortran 5.0编制而成。
2、程序流程及图框启动输入原始数据自动划分网格形成MA计算NNHMX形成整体刚度矩阵K形成荷载列向量RLU分解KLU回代并求得结点位移输入结点位移计算单元应力及主应力等整体应力磨平结点应力停机图2- 程序流程图-3- MAINPROGRAMINPUTHUAFENCBANDSKOSTIFFDNXFUN8CONCRBODYRFA CERDECOPFOBASTRESSGAUSSSTRESSSUMSSUMSTRSOUTDISTRE 图2-子程序框图其中各子程序的主要功能为INPUT――输入原始数据HUAFEN――自动网格划分形成COOR2NPXY的坐标值与单元信息CBAND――形成主元素序号指示矩阵MA SKO――形成整体刚度矩阵K CONCR――计算集中力引起的等效结点荷载Re BODYR――计算自重体力引起的等效结点荷载Re FACER――计算分布面力引起的等效结点荷载Re DECOP――支配方程LU三角分解FOBA――LU分解直接解法中的回代过程OUTDISP――输出结点位移分量STRESS――计算单元应力分量OUTSTRE――输出单元应力分量STIF――计算单元刚度矩阵FDNX――计算形函数对整体坐标的导数TiiyNxNi1234。
FUN8――计算形函数及雅可比矩阵J SFUN ――应力磨平-单元下的…K‟NCN… SCN――应力磨平-单元下的右端项系数…CN… SUMSKN――应力磨平-单元下的右端项集成到总体的…P… -4- SUMSTRS――应力磨平-单元下的集成到总体的…K… GAUSTRSS――高斯消元求磨平后的应力3、输入数据及变量说明当程序开始运行时按屏幕提示键入数据文件的名字。
matlab四节点四边形等参元的刚度矩阵计算程序
-1.2150 0.1440 2.1968 -0.8933 0.5212 -0.0151 -1.5030 0.7645
0.0616 -0.0419 -0.8933 1.8399 0.0673 -0.5250 0.7645 -1.2729
x2=2;y2=0;
x3=2.25;y3=1.5;
x4=1.25;y4=1;
材料参数：
E=30e12;
NU=0.3;
h=1;
ID=1;
Ai=1
Aj=1;
3、程序（由附件给出）
4、计算结果及讨论（需有图表来说明）
计算结果如下：
k =
1.0e+012 *
1.4619 0.2211 -1.2150 0.0616 -0.3716 -0.4236 0.1248 0.1409
-0.3716 -0.4236 0.5212 0.0673 1.1645 0.2763 -1.3141 0.0800
-0.4236 -0.2582 -0.0151 -0.5250 0.2763 0.9061 0.1624 -0.1229
0.1248 0.0585 -1.5030 0.7645 -1.3141 0.1624 2.6923 -0.9854
%输入平面问题性质参数ID（1为平面应力，2为平面应变）
%输出单元刚度矩阵
%-------------------------------------------------------------------
symsst;
a=[-(1-t)*x1+(1-t)*x2+(1+t)*x3-(1+t)*x4]/4;
平面四边形四节点等参单元Fortran源程序复习过程
平面四边形四节点等参单元F o r t r a n源程序C ************************************************C * FINITE ELEMENT PROGRAM *C * FOR Two DIMENSIONAL ELASticity PROBLEM *C * WITH 4 NODE *C ************************************************PROGRAM ELASTICITYcharacter*32 dat,cchDIMENSION SK(80000),COOR(2,300),AE(4,11),MEL(5,200),& WG(4),JR(2,300),MA(600),R(600),iew(30),STRE(3,200)COMMON /CMN1/ NP,NE,NM,NRCOMMON /CMN2/ N,MX,NHCOMMON /CMN3/ RF(8),SKE(8,8),NN(8)WRITE(*,*)'PLEASE ENTER INPUT FILE NAME'READ(*,'(A)')DATOPEN(4,FILE=dat,STATUS='OLD')OPEN(7,FILE='OUT',STATUS='UNKNOWN')READ(4,*)NP,NE,NM,NRWRITE(7,'(A,I6)')'NUMBER OF NODE---------------------NP=',npWRITE(7,'(A,I6)')'NUMBER OF ELEMENT------------------NE=',ne WRITE(7,'(A,I6)')'NUMBER OF MATERIAL-----------------NM=',nm WRITE(7,'(A,I6)')'NUMBER OF surporting---------------NC=',Nr CALL INPUT (JR,COOR,AE,MEL)CALL CBAND (MA,JR,MEL)DO I=1,NHSK(I)=0.0enddoCALL SK0(SK,MEL,COOR,JR,MA,AE)do I=1,NR(I)=0.0enddopause 'aaa'stopREAD(4,*)NCP,NBE,izWRITE(*,'(5i8)')NCP,NBE,izWRITE(7,'(5i8)')NCP,NBE,izIF(NCP.GT.0)CALL CONCR(NCP,R,JR)IF(NBE.GT.0) CALL BODYR(NBE,R,MEL,COOR,JR,AE)IF(iz.GT.0)thendo jj=1,izREAD (4,*)Js,nse,(WG(I),I=1,4)read(4,*)(iew(m),m=1,nse)CALL FACER(iew,NSE,R,MEL,COOR,JR,WG)enddoendifCALL DECOP (SK,MA)CALL FOBA (SK,MA,R)CALL OUTDISP(NP,R,JR)CALL STRESS (COOR,MEL,JR,AE,R,STRE)WRITE(7,'(A)')' PROGRAM SAFF HAS BEEN ENDED'WRITE(*,'(A)')' PROGRAM SAFF HAS BEEN ENDED'STOPc RETURNENDC *********************************************SUBROUTINE INPUT (JR,COOR,AE,MEL)DIMENSION JR(2,*),COOR(2,*),AE(4,*),MEL(5,*)COMMON /CMN1/ NP,NE,NM,NRCOMMON /CMN2/ N,MX,NHDO 70 I=1,NPREAD(4,*) IP,X,YCOOR(1,IP)=XCOOR(2,IP)=Y70 CONTINUEDO 11 J=1,NEREAD(4,*)NEE,NME,(MEL(I,NEE),I=1,4)MEL(5,NEE)=NME11 CONTINUEDO 10 I=1,NPDO 10 J=1,210 JR(J,I)=1DO 20 I=1,NRREAD(4,*) IP,IX,IYJR(1,IP)=IXJR(2,IP)=IY20 CONTINUEN=0DO 30 I=1,NPDO 30 J=1,2IF (JR(J,I)) 30,30,2525 N=N+1JR(J,I)=N30 CONTINUEDO 55 J=1,NMREAD (4,*)JJ,(AE(I,JJ),I=1,4)WRITE(*,910) JJ,(AE(I,JJ),I=1,4)55 CONTINUE910 FORMAT (/20X,'MATERIAL PROPERTIES'/(3X,I5,4(1x,E8.3))) RETURNENDC **********************************************SUBROUTINE CBAND (MA,JR,MEL)DIMENSION MA(*),JR(2,*),MEL(5,*),NN(8)COMMON /CMN1/ NP,NE,NM,NRCOMMON /CMN2/ N,MX,NHDO 65 I=1,N65 MA(I)=0DO 90 IE=1,NEDO 75 K=1,4IEK=MEL(K,IE)DO 95 M=1,2JJ=2*(K-1)+MNN(JJ)=JR(M,IEK)95 CONTINUE75 CONTINUEL=NDO 80 I=1,2*4NNI=NN(I)IF(NNI.EQ.0) GO TO 80IF(NNI.LT.L) L=NNI80 CONTINUEDO 85 M=1,2*4JP=NN(M)IF(JP.EQ.0) GO TO 85JPL=JP-L+1IF(JPL.GT.MA(JP)) MA(JP)=JPL85 CONTINUE90 CONTINUEMX=0MA(1)=1DO 10 I=2,NIF(MA(I).GT.MX) MX=MA(I)MA(I)=MA(I)+MA(I-1)10 CONTINUENH=MA(N)WRITE(7,'(A,I8)')'TOTAL DEGREES OF FREEDOM-----------N= ',N WRITE(7,'(A,I8)')'MAX-SEMI-BANDWIDTH-----------------MX=',MX WRITE(7,'(A,I8)')'TOTAL-STORAGE----------------------NH=',NH500 FORMAT (/5X,'FREEDOM N='*,I5,3X,'SEMI-BANDWI. MX=',I5,3X,* 'STORAGE NH=',I7)RETURNENDC **********************************************SUBROUTINE SK0(SK,MEL,COOR,JR,MA,AE)DIMENSION SK(*),MEL(5,*),COOR(2,*),JR(2,*),MA(*), * AE(4,*),XYZ(2,4),iven(4)COMMON /CMN1/ NP,NE,NM,NRCOMMON /CMN2/ N,MX,NHCOMMON /CMN3/ RF(8),SKE(8,8),NN(8)COMMON /CMN4/ NEE,NMECOMMON /GAUSS/ RSTG(3),H(3)H(1)=0.5555555555555560H(2)=0.8888888888888890H(3)=H(1)RSTG(1)=-0.7745966692414830RSTG(2)=0.00RSTG(3)=-RSTG(1)DO 10 IE=1,NENEE=IENME=MEL(5,IE)DO 75 K=1,4IEK=MEL(K,IE)iven(k)=IEKDO 95 M=1,2JJ=2*(K-1)+MNN(JJ)=JR(M,IEK)95 XYZ(M,K)=COOR(M,IEK)75 CONTINUECALL STIF(XYZ,AE,iven)DO 60 I=1,8DO 60 J=1,8II=NN(I)JJ=NN(J)IF ((JJ.EQ.0).OR.(II.LT.JJ)) GO TO 60JN=MA(II)-(II-JJ)SK(JN)=SK(JN)+SKE(I,J)60 CONTINUE70 CONTINUEwrite(7,1111) ((ske(i,j),j=1,8),i=1,8)1111 format(2x,8f12.2)10 CONTINUERETURNENDC *********************************************SUBROUTINE STIF(XYZ,AE,iven)DIMENSION AE(4,*),DNX(2,4),XYZ(2,*),iven(*),* RJAC(2,2)COMMON /CMN1/ NP,NE,NM,NRCOMMON /CMN2/ N,MX,NHCOMMON /CMN3/ RF(8),SKE(8,8),NN(8)COMMON /CMN4/ NEE,NMECOMMON /GAUSS/ RSTG(3),H(3)DO 40 I=1,8RF(I)=0.00DO 30 J=1,8SKE(I,J)=0.0030 CONTINUE40 CONTINUEE=AE(1,NME)U=AE(2,NME)GAMA=AE(3,NME)D1=E*(1.00-U)/((1.00+U)*(1.00-2.00*U))D2=E*U/((1.00+U)*(1.00-2.00*U))D3=E*0.50/(1.00+U)DO 120 I=1,4II=2*(I-1)I1=II+1I2=II+2DO 115 J=1,4JJ=2*(J-1)J1=JJ+1J2=JJ+2DXX=0DXY=0DYX=0DYY=0DO 99 IS=1,3S=RSTG(IS)SH=H(IS)DO 98 IR=1,3R=RSTG(IR)RH=H(IR)CALL FDNX (XYZ,DNX,DET,R,S,RJAC,iven,NEE) DNIX=DNX(1,I)DNIY=DNX(2,I)DNJX=DNX(1,J)DNJY=DNX(2,J)DXX=DXX+DNIX*DNJX*DET*RH*SHDXY=DXY+DNIX*DNJY*DET*RH*SHDYX=DYX+DNIY*DNJX*DET*RH*SHDYY=DYY+DNIY*DNJY*DET*RH*SH98 CONTINUE99 CONTINUESKE(I1,J1)=DXX*D1+DYY*D3SKE(I2,J2)=DYY*D1+DXX*D3SKE(I1,J2)=DXY*D2+DYX*D3SKE(I2,J1)=DYX*D2+DXY*D3115 CONTINUE120 CONTINUERETURNENDC ********************************************* SUBROUTINE CONCR(NCP,R,JR)DIMENSION R(*),JR(2,*),XYZ(2)DO 100 I=1,NCPREAD (4,*) IP,PX,PYXYZ(1)=PXXYZ(2)=PYDO 95 J=1,2L=JR(J,IP)IF(L.EQ.0) GO TO 95R(L)=R(L)+XYZ(J)95 CONTINUE100 CONTINUERETURNENDC ********************************************** SUBROUTINE BODYR(NBE,R,MEL,COOR,JR,AE) DIMENSION R(*),MEL(5,*),COOR(2,*),JR(2,*),& AE(4,*),XYZ(2,4),iven(4)COMMON /CMN1/ NP,NE,NM,NRCOMMON /CMN2/ N,MX,NHCOMMON /CMN3/ RF(8),SKE(8,8),NN(8)COMMON /CMN5/ FUN(4),PN(2,4),XJAC(2,2)COMMON /GAUSS/ RSTG(3),H(3)H(1)=1.0H(2)=1.0RSTG(1)=-0.5773502691896260RSTG(2)=-RSTG(1)DO 10 IE=1,NBEDO I=1,8RF(I)=0.00ENDDOc READ(4,*)NEENEE=ieNME=MEL(5,NEE)GAMA=AE(3,NME)DO 75 K=1,4IEK=MEL(K,NEE)iven(k)=iekDO 95 M=1,2JJ=2*(K-1)+MNN(JJ)=JR(M,IEK)95 XYZ(M,K)=COOR(M,IEK)75 CONTINUEDO 99 IS=1,2S=RSTG(IS)SH=H(IS)DO 98 IR=1,2RR=RSTG(IR)RH=H(IR)CALL FUN8 (XYZ,RR,S,DET)DO 30 I=1,4J=2*IRF(J)=RF(J)-FUN(I)*RH*SH*DET*GAMA30 CONTINUE98 CONTINUE99 CONTINUECALL ASLOAD (R)10 CONTINUERETURNENDC *********************************************SUBROUTINE FACER(iew,NSE,R,MEL,COOR,JR,WG) DIMENSION R(*),MEL(5,*),COOR(2,*),JR(2,*),wg(*) * ,XYZ(2,4),iew(*),PR(2)COMMON /CMN1/ NP,NE,NM,NRCOMMON /CMN2/ N,MX,NHCOMMON /CMN3/ RF(8),SKE(8,8),NN(8)COMMON /CMN4/ NEE,NMECOMMON /GAUSS/ RSTG(3),H(3)H(1)=1.0H(2)=1.0RSTG(1)=-0.5773502691896260RSTG(2)=-RSTG(1)nwf=0nnf=0ir=wg(1)+0.1if(ir.eq.1)nwf=1if(ir.eq.2)nnf=1DO 510 IE=1,NSEDO I=1,8RF(I)=0.00ENDDOnee=iew(ie)DO 575 K=1,4IEK=MEL(K,NEE)DO 595 M=1,2JJ=2*(K-1)+MNN(JJ)=JR(M,IEK)595 XYZ(M,K)=COOR(M,IEK)575 CONTINUEIF(NWF.EQ.1) thenGAMA=WG(2)Z0=WG(3)NSU=WG(4)+0.1CALL SURLOD (NSU,XYZ,PR,Z0,GAMA,1)endifIF(NNF.EQ.1) thenq=WG(2)NSU=WG(4)+0.1do j=1,2PR(J)=qenddoCALL SURLOD (NSU,XYZ,PR,Z0,GAMA,2)endifCALL ASLOAD (R)510 CONTINUERETURNENDC *********************************************SUBROUTINE SURLOD (NSU,XYZ,PR,Z0,GAMA,NSI)DIMENSION XYZ(2,*),RST(3),PR(2),KCRD(4),KFACE(2,4), & FVAL(4),NODES(2),FACT(4)COMMON /CMN1/ NP,NE,NM,NRCOMMON /CMN2/ N,MX,NHCOMMON /CMN3/ RF(8),SKE(8,8),NN(8)COMMON /CMN4/ NEE,NMECOMMON /CMN5/ FUN(4),PN(2,4),XJAC(2,2)COMMON /GAUSS/ RSTG(3),H(3)DATA KCRD/1,1,2,2/DATA KFACE/1, 4,* 2, 3,* 1, 2,* 4, 3/DATA FVAL/-1.00,1.00,-1.00,1.00/FACT(1)=1.0FACT(2)=-1.0FACT(3)=-1.0FACT(4)=1.0FACTNUS=FACT(NSU)DO I=1,2J=KFACE(I,NSU)NODES(I)=JENDDOIF (NSI.EQ.1) THENDO I=1,2J=NODES(I)Z=Z0-XYZ(2,J)PR(I)=0.00IF (Z.GT.0.00) PR(I)=Z*GAMAENDDOENDIFML=KCRD(NSU)IF(ML.EQ.1)MM=2IF(ML.EQ.2)MM=1RST(ML)=FVAL(NSU)DO 70 LX=1,2RST(MM)=RSTG(LX)CALL FUN8 (XYZ,RST(1),RST(2),DET) PXYZ=0.00DO 25 I=1,2J=NODES(I)PXYZ=PXYZ+FUN(J)*PR(I)25 CONTINUEA1=XJAC(MM,2)A2=-XJAC(MM,1)30 DO 60 I=1,2J=NODES(I)K2=2*JK1=K2-1Q=PXYZ*FUN(J)*H(LX)*FACTNUSRF(K1)=RF(K1)+Q*A1RF(K2)=RF(K2)+Q*A260 CONTINUE70 CONTINUEENDC ********************************************* SUBROUTINE ASLOAD (R)DIMENSION R(*)COMMON /CMN1/ NP,NE,NM,NRCOMMON /CMN3/ RF(8),SKE(8,8),NN(8)DO 20 I=1,8L=NN(I)IF (L.EQ.0) GO TO 20R(L)=R(L)+RF(I)20 CONTINUERETURNENDC *********************************************** SUBROUTINE DECOP (SK,MA)DIMENSION SK(*),MA(*)COMMON /CMN2/ N,MX,NHDO 50 I=2,NL=I-MA(I)+MA(I-1)+1K=I-1L1=L+1IF (L1.GT.K) GO TO 30DO 20 J=L1,KIJ=MA(I)-I+JM=J-MA(J)+MA(J-1)+1IF (L.GT.M) M=LMP=J-1IF (M.GT.MP) GO TO 20DO 10 LP=M,MPIP=MA(I)-I+LPJP=MA(J)-J+LPSK(IJ)=SK(IJ)-SK(IP)*SK(JP)10 CONTINUE20 CONTINUE30 IF (L.GT.K) GO TO 50DO 40 LP=L,KIP=MA(I)-I+LPLPP=MA(LP)SK(IP)=SK(IP)/SK(LPP)II=MA(I)SK(II)=SK(II)-SK(IP)*SK(IP)*SK(LPP)40 CONTINUE50 CONTINUEENDC *************************************************************SUBROUTINE FOBA (SK,MA,R)DIMENSION SK(*),MA(*),R(*)COMMON /CMN2/ N,MX,NHDO 10 I=2,NL=I-MA(I)+MA(I-1)+1K=I-1IF (L.GT.K) GO TO 10DO 5 LP=L,KIP=MA(I)-I+LPR(I)=R(I)-SK(IP)*R(LP)5 CONTINUE10 CONTINUEDO 20 I=1,NII=MA(I)45 R(I)=R(I)/SK(II)20 CONTINUEDO 30 J1=2,NI=2+N-J1L=I-MA(I)+MA(I-1)+1K=I-1IF (L.GT.K) GO TO 30DO 25 J=L,KIJ=MA(I)-I+J55 R(J)=R(J)-SK(IJ)*R(I)25 CONTINUE30 CONTINUERETURNENDC ***************************************************************** SUBROUTINE STRESS(COOR,MEL,JR,AE,R,STRE)DIMENSION XYZ(2,4),DNX(2,4),AE(4,*),STRE(3,*),& COOR(2,*),MEL(5,*),JR(2,*),RJAC(2,2),SIG(3),& B(3,8),R(*),iven(4)COMMON /CMN1/ NP,NE,NM,NRCOMMON /CMN3/ RF(8),SKE(8,8),NN(8)COMMON /CMN5/ FUN(4),PN(2,4),XJAC(2,2)DO 106 IE=1,NENME=MEL(5,IE)DO 300 K=1,4IEK=MEL(K,IE)DO 310 M=1,2310 XYZ(M,K)=COOR(M,IEK)DO 320 M=1,2JRR=2*(K-1)+M320 NN(JRR)=JR(M,IEK)300 CONTINUEE=AE(1,NME)U=AE(2,NME)D1=E*(1.00-U)/((1.00+U)*(1.00-2.00*U))D2=E*U/((1.00+U)*(1.00-2.00*U))D3=0.50*E/(1.00+U)SS=0.0RR=0.0CALL FDNX (XYZ,DNX,DET,RR,SS,RJAC,iven,IE) DO 30 I=1,4II=2*(I-1)J1=II+1J2=II+2BI=DNX(1,I)CI=DNX(2,I)B(1,J1)=BIB(2,J1)=0.B(3,J1)=CIB(1,J2)=0.B(2,J2)=CIB(3,J2)=BI30 CONTINUEDO 55 II=1,3SIG(II)=0.0055 CONTINUEDO 70 K=1,8NA=NN(K)IF (NA.EQ.0) GO TO 70DO 60 L=1,3SIG(L)=SIG(L)+B(L,K)*R(NA)60 CONTINUE70 CONTINUESX=D1*SIG(1)+D2*SIG(2)SY=D2*SIG(1)+D1*SIG(2)SXY=D3*SIG(3)STRE(1,IE)=SXSTRE(2,IE)=SYSTRE(3,IE)=SXY106 CONTINUECALL OUTSTRE(NE,STRE)RETURNENDC *********************************************SUBROUTINE FDNX (XYZ,DNX,DET,R,S,RJAC,iven,NEE) DIMENSION XYZ(2,*),DNX(2,*),RJAC(2,2),iven(*)COMMON /CMN5/ FUN(4),PN(2,4),XJAC(2,2)CALL FUN8 (XYZ,R,S,DET)IF (DET.LT.1.0E-5)THENWRITE(7,600) NEE,R,S,detWRITE(7,*) (iven(m),m=1,4)STOPENDIFREC=1.00/DETRJAC(1,1)=REC*XJAC(2,2)RJAC(2,2)=REC*XJAC(1,1)RJAC(2,1)=-REC*XJAC(2,1)RJAC(1,2)=-REC*XJAC(1,2)DO 30 K=1,4DO 20 I=1,2DNX(I,K)=0.DO 25 M=1,2DNX(I,K)=DNX(I,K)+RJAC(I,M)*PN(M,K)25 CONTINUE20 CONTINUE30 CONTINUE600 FORMAT (1X,'ERR0R*** NEGTIVE OR ZERO '* 'JACOBIAN DETERMINANT FOR '* 'ELEMENT'/'ELE.=',I5,' R=',F10.5,6X,'S=',F10.5,* 'det=',f12.5)RETURNENDC *********************************************SUBROUTINE FUN8 (XYZ,R,S,DET)DIMENSION XYZ(2,*),XI(4),ETA(4)COMMON /CMN5/ FUN(4),PN(2,4),XJAC(2,2)DATA XI/-1.0,1.0,1.0,-1.0/DATA ETA/-1.0,-1.0,1.0,1.0/DO 10 I=1,4G1=(1.0+XI(I)*R)G2=(1.0+ETA(I)*S)FUN(I)=0.25*G1*G2PN(1,I)=0.25*XI(I)*G2PN(2,I)=0.25*ETA(I)*G110 CONTINUEDO 80 I=1,2DO 75 J=1,2DET=0.00DO 70 K=1,4DET=DET+PN(I,K)*XYZ(J,K)70 CONTINUEXJAC(I,J)=DET75 CONTINUE80 CONTINUEDET=XJAC(1,1)*XJAC(2,2)* -XJAC(2,1)*XJAC(1,2)RETURNENDC ********************************************** SUBROUTINE OUTDISP(NP,R,JR)DIMENSION R(*),JR(2,*),U(2)WRITE(*,650)WRITE(7,650)DO I=1,NPDO M=1,2L=JR(M,I)IF(L.EQ.0)U(M)=0.0IF(L.GT.0)U(M)=R(L)ENDDOWRITE(*,'(5X,I5,10X,2E14.3)') I,UWRITE(7,'(5X,I5,10X,2E14.3)') I,UENDDO650 FORMAT(/25X,'NODAL DISPLACEMENTS'/8X, * 'NODE',13X,'X-COMP.',8X,'Y-COMP.')RETURNENDC ********************************************** SUBROUTINE OUTSTRE(NE,STRE)DIMENSION STRE(3,*),ST(6)WRITE(*,700)WRITE(7,700)DO IE=1,NESX=STRE(1,IE)SY=STRE(2,IE)SXY=STRE(3,IE)ST(1)=SXST(2)=SYST(3)=SXYH1=SX+SYH2=SQRT((SX-SY)*(SX-SY)+4.0*SXY*SXY)ST(4)=(H1+H2)/2.0ST(5)=(H1-H2)/2.0IF(ABS(SXY).LT.1.0E-4)THENIF (SX.GT.SY) ST(6)=0.0IF (SX.LE.SY) ST(6)=90.0ELSEST(6)=ATAN((ST(4)-SX)/SXY)*57.29578ENDIFWRITE(*,'(6X,I4,3X,6F11.3)') IE,STWRITE(7,'(6X,I4,3X,6F11.3)') IE,STENDDO700 FORMAT(/30X,'ELEMENT STRESSES'/5X,*'ELEMENT',5X,'X-STRESS',3X,'Y-STRESS',*2X,'XY-STRESS',1X,'MAX-STRESS',1X,*'MIN-STRESS',4X,'ANGLE'/)RETURNENDC *********************************************。
有限单元法及程序设计
有限单元法及程序设计有限单元法（Finite Element Method，FEM）是一种用于数值分析和计算的方法，广泛应用于工程和科学领域。
它通过将连续问题离散化成有限个小单元，并在每个小单元上建立数学模型来近似求解问题。
本文将介绍有限单元法的基本原理、步骤以及程序设计方面的注意事项。
一、有限单元法基本原理有限单元法的基本原理是将连续的物理区域划分为有限个离散的小单元，每个小单元内的场量近似表示为一些插值函数的线性组合。
通过对这些小单元进行逐个求解，最终得到整个问题的近似解。
有限单元法的核心思想是利用局部性原则，将整个问题分解成多个小问题。
每个小问题只涉及到相邻的单元，在确定了边界条件和材料特性后，可以进行独立的求解。
最后通过组合各个小问题的解，得到整个问题的解。
二、有限单元法步骤有限单元法的求解过程主要包括几个基本步骤，具体如下：1. 离散化：将连续的物理区域划分为有限的小单元。
常用的小单元形状包括三角形、四边形、六边形等。
2. 建立数学模型：在每个小单元上建立数学模型，通常使用插值函数来近似表示物理量。
插值函数的选择对求解结果的准确性和效率有重要影响。
3. 形成总体方程：根据物理规律和边界条件，利用适当的数学方法推导出总体方程。
常见的总体方程包括稳定性方程、运动方程等。
4. 矩阵装配：将每个小单元的局部方程装配成整个系统的总体方程。
这一步骤常常需要对单元进行编号和排序，以便正确地装配矩阵。
5. 边界条件处理：根据实际问题的边界条件，对总体方程进行修正。
边界条件的处理通常包括施加约束和设定边界值。
6. 求解方程：通过数值方法，如有限差分法或有限元法，求解总体方程。
常用的求解方法包括直接法和迭代法。
7. 后处理：对求解结果进行计算和分析，以获得实际问题的有用信息。
后处理包括输出位移、应力、应变等字段，以及进行可视化展示。
三、程序设计注意事项在进行有限单元法的程序设计时，需要充分考虑以下几个方面的注意事项：1. 算法选择：根据问题的特点和求解需求，选择合适的有限单元类型、插值函数和数值解法。
有限元-4-等参
表5-1 三结点三角形单元与四结点矩形单元比较
单元类型
优点
缺点
三 结 点 三 角 适应复杂形状，
形单元
单元大小过渡方便
计算精度低
四结点矩形 单元
单元内的应力、应变 是线性变化的，计算 精度较高
不能适应曲线边界和 非正交的直线边界
❖ 如果任意形状的四边形四结点单元采用矩形单元的位 移模式，则在公共边界上不满足位移连续性条件。为 了既能得到较高的计算精度，又能适应复杂的边界形
❖ 参照矩形单元，四结点正方形单元的位 移模式为，
❖ u N1u1 N2u2 N3u3 N4u4
❖ v N1v1 N2v2 N3v3 N4v4 （4-4） ❖ 其中，
❖
N1
1 4
(1 )(1 )
❖ ❖ ❖
N2
1 (1 )(1 )
4
N3
1 (1 )(1 )
4
（4-5）
❖
N4
1 (1 )(1 )
4
❖ 四个结点的坐标为(i ,i ) ，定义新的变量，
❖ 0 i ，0 i
(i=1，2，3，4)
❖ 形态函数表示为，
❖
Ni
1 4
(1
0
)(1
0
)
（i=1，2，3，4）
(4-6) （4-7）
把及 作为任意四边形单元的局部坐标，把（4-4）的位移模 式和（4-7）的形态函数用于任意形状的四边单元，可得：
1）在四个结点处可以得到结点的位移；
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
［1 ］ 不连 续 问 题 。 XFEM 是 基 于 单 位 分 解 的 方 法 （ PUM） 对单元的形函数加以改进， 从而考虑所研究
2012 年 第 31 卷
12 月 第 12 期
机械科学与技术 Mechanical Science and Technology for Aerospace Engineering
December Vol． 31
2012 No． 12
四结点等参元 XFEM 程序设计及 在裂纹问题中的应用
苏
苏 毅
毅， 王生楠， 闫晓中
4 4
∫ σ·δεdΩ = ∫ t·δudΓ + ∫ b·δudΩ
Ω Γt Ω
（ 3）
b 分别为分布体力、 式中： t， 分布面力。 位移模式构 造后， 和常规有限元方法一样， 将式 （ 2 ） 代入式 （ 3 ） 由虚功原理推导扩展有限元的支配方程 Kδ = F （ 4）
式中： δ 是结点未知向量； K 和 F 分别为总体刚度矩 K 和 F 先逐个单元计算， 阵和总体荷载列阵， 再按常 单元层次上的 k 和 f 分别为 K I 单元 规步骤进行组装。 劲度矩阵 k ［ k ij］ = B， B， B］ ［B ，
XFEM Programming for Four Nodes Isoparametric Element and Its Application in Crack Problems
Su Yi，Wang Shengnan，Yan Xiaozhong
（ School of Aeronautics，Northwestern Polytechnical University，Xi'an 710072 ）
( )
( )
( )
( )
b j2 i ) 为改进结点两个方向的
T
只有第 附加自由度。当结点为常规单元的结点时， 一项； 当结点为被裂纹贯穿单元的结点时 ， 有第一项 和第二项； 当结点为裂纹尖端所在单元的结点时 ， 有 第一项和第三项； 结点同时处于裂纹贯穿单元和裂 纹尖端所在单元， 应优先属于裂纹尖端所在单元， 加 强方式选择裂纹尖端所在单元结点加强方式 。 2 ） 虚功方程和支配方程
T
为结点与 Heaviside 函数有关的加强自由度 ； 第三 项只有在点 x 处在裂尖单元 、 单元结点需要改进时 当第三项有意义时 β i = 才有意义 ； β i 为指示函数 ， 1， ； （ x ） x 否则为零 φ j 为点 的改进函数 ｛ φ（ r ， r sin θ ， r cos θ ， r sin θ sin θ ， θ） ｝ = ｛ 槡 槡 槡 2 2 2 r cos θ sin θ ｝ ( b j1 i 槡 2
（ 西北工业大学 航空学院， 西安 710072 ）
要： 虽然扩展有限元法( XFEM) 在处理裂纹等这种强不连续性 问题 时 理 论 上 是 成 功 的， 得到很 快的发展和广泛的应用， 但在实际应用中， 尚存在许多技术问题如网格密度等值得研究。为了验证 摘 和提高 XFEM 在 计算 裂 纹 应力 强 度因子 上 的 有 效性， 针 对 四 结 点 等 参 元 推 导 了 XFEM 的 相 应 公 式， 编写了用于计算含裂 纹板 裂 纹 尖端 应力 强 度因子的 完 整 的 Matlab 程 序。 针 对 典 型 含裂 纹 平 板， 采用本文编写的程序计算了裂纹尖端应力强度因子， 与采用传统有限元法的结果进行了对比分 XFEM 在 计算 裂 纹 尖端 应力 强 度 析， 并进一步研究了网格参数， 对 XFEM 结果 的 影响。 结果 表明， 因子上有很好的计算精度， 但其计算结果对网格密度较为敏感， 在实际应用中应当引起重视。 关 键 词： 扩展有限元; 裂纹尖端; 应力强度因子; 四结点等参元; 网格密度 中图分类号： O346. 1 文献标识码： A 8728 （ 2012 ） 12194906 文章编号： 1003-
v i 为点 x 所 式中 ： N i （ x ） 为标准有限元形函数 ； u i 、 4 （ u 、 v 在单元 个结点位移分量 i i 并非最终的结点 单 位移值 ） ； 第二项只有在点 x 处在裂纹贯穿单元 、 元结点需要改进时才有意义 ； α i 为指示函数 ， 当第 二项有意义时 α i = 1 ， 否则为零 ； H 为跳跃函数 ， 在 局部坐标系下定义为 H = +1 {－ 1 y ＞ 0 y ＜ 0 ； （ a1i ， a2i ）
x =
Ni xi ， y ∑ i =1
=
Ni yi ∑ i =1
（ 1）
XFEM 中四结点等参单元的位 关于裂纹问题， 移模式为
[u ] = ∑ N （ x） [ v ] + ∑ α N （ x） H [ a v
4
ui
i
4
a1i
2i
i
i
i
i =1 4
i =1
]
+
∑ β i N i （ x） φj
i =1
1 2 3 4 Ω T 48 ×3
［ D］ B1 ， B2 ， B3 ， B4］ 3 ×3［ 3 ×48 tdxdy =
uu k14 au k14 ua k14 aa k14 b ja k14 ub j k14 ab j k14 b jb j k14 j =1 ～4 ， ub j k44 ab j k44 b jb j k44
[ ]
b j1i b 2i
j
，j = 1 ～ 4
（ 2）
k11 au k11 kbju 11 uu k41 au k41 au k41
uu
ua k11 aa k11 b ja k11
ub j k11 ab j k11 b jb j k11
… … … …
问题的不连续、 奇异性等特性。 XFEM 所使用的网 格与结构内部几何或物理界面无关， 从而克服了裂 纹尖端等高应力和变形集中区内网格划分的困难 ， 使得在模拟裂纹生长过程中无需对网格进行重新划
1950
机械科学与技术
第 31 卷
分。自 XFEM 问世以来， 在国际上得到了很快的发 ［2 ～ 6 ］ 。 文献［ 7］ 展和广泛的应用 综述了 XFEM 在静 态 和 扩 展 裂 纹 问 题 中 的 应 用， 并与广义有限元 （ GFEM） 进行了比较， 8］ 结果符合的很好。文献［ 采 用 XFEM 研究了双材料界面裂纹问题的应力强度 9］综述了 XFEM 的基本思想、 因子的计算。 文献［ 实施步骤及其应用， 初步展望了该领域发展需要研 究的课题。 尽管 XFEM 在处理裂纹等这种强不连续性问 题时理论上是成功的， 但在方法的实现和应用上， 尚 存在许多技术如网格密度、 结点加强区域范围等问 题值得研究。为了验证扩展有限元法在计算裂纹应 力强度因子上的有效性， 在已有 XFEM 理论的基础 上， 推导了四结点等参元相应的计算公式 ， 编写了完 对 XFEM 整的 Matlab 程序。 针对典型含裂纹平板， 的应用进行了研究， 并进一步研究了其计算结果对 网格密度的敏感性。 1 XFEM 中四结点等参元的计算 SIF
Abstract： It was successful that extended finite element method （ XFEM ） was theoretically applied to the strong discontinuity problems such as cracks so that the development and application of the method are being increased rapidly． However，in practical applications，there are still many technical issues such as the mesh density to be studied． In order to verify and improve the effectiveness of XFEM on the calculation of crack stress intensity factor， the corresponding formula of XFEM with fournode isoparametric element was derived and a complete Matlab code was also edited aiming at calculating the crack tip stress intensity factor of a plate with crack． Based upon the program，the crack tip stress intensity factors for a typical plate with crack were calculated and the results by XFEM were compared with those by the traditional finite element method （ FEM ） ． The effect of the mesh parameters on the XFEM results was further studied． The study shows that XFEM has a very good accuracy in the calculation of crack tip stress intensity factor，but the results by XFEM are sensitive to the density of the mesh，this should be paid more attention to practical applications． Key words： XFEM； crack tip； stress intensity factor； four nodes isoparametric element； mesh density； discontinuity problem； cracks 1999 年， 以美国西北大学 Belytschko 教授为代 表的研究组首先提出用扩展有限元 （ XFEM ） 来解决