第九章平面广义四节点等参元 GQ4

合集下载

有限元讲义弹性力学平面问题有限单元法(四结点四边形等参元,八结点曲线四边形等参元,问题补充)分析

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 与位移函数中的完全一致。

有限元 2-4弹性力学平面问题(2.6四结点四边形等参元,2.7八结点曲线四边形等参元,2.8几个问题补充 )

将(2-6-6)代入即可获得[B], [B]是ξ,η的函数。
有限单元法
土木工程学院
P-17/73
三、单元刚度矩阵
获得[B]后,便可由单刚的一般表达式:
[K ] = t ∫∫ [B ] [D ][B ]dxdy
T
求出四结点四边形的单元刚度矩阵。在按上述公式作积分运算时，必须把面积元 dxdy变换成dξdη,图a上的面积元abdc的面积等于 → → 矢量 ac 与矢量 ab 的矢量积的模，即微元为
→
→
I × I = J× J = 0
→
I× J =1
→
则有:
→ → ⎞ ⎛ ∂x → →⎞ ⎛ ∂x ∂y ∂y dA = ⎜ ⎟ ⎜ ∂η dη I + ∂η dη J ⎟ ⎟×⎜ ⎜ ∂ξ dξ I + ∂ξ dξ J ⎟ ⎠ ⎠ ⎝ ⎝ → → ⎛ ∂x ∂y ∂y ∂x ⎞ =⎜ ⎟dξdη I × J = J dξdη ⎜ ∂ξ ∂η − ∂ξ ∂η ⎟ ⎠ ⎝
∂N i ∂N i ∂x ∂N i ∂y = + ∂ξ ∂x ∂ξ ∂y ∂ξ
∂N i ∂N i ∂x ∂N i ∂y = + ∂x ∂η ∂y ∂η ∂η
(2 − 6 − 3)
有限单元法
土木工程学院
P-12/73
或写成：
⎧ ∂N i ⎫ ⎧ ∂N i ⎫ ⎪ ⎪ ⎪ ∂ξ ⎪ ⎪ ⎪ ∂x ⎪ ⎪ ⎨ ∂N ⎬ = [J ]⎨ ∂N ⎬ ⎪ i⎪ ⎪ i⎪ ⎪ ⎪ ⎩ ∂y ⎪ ⎭ ⎩ ∂η ⎪ ⎭
引入边界条件，即可得位移函数：
U = ∑ N iU i
ijmp
V = ∑ N iVi
ijmp
有限单元法
土木工程学院

c++ 二维四节点四边形等参元刚度矩阵

C++ 二维四节点四边形等参元刚度矩阵在计算机辅助工程领域，C++语言被广泛应用于有限元分析（FEA）和计算力学等方面。

而在有限元分析中，等参元（Isoparametric Element）是一种常用的元素类型，用于对复杂的结构和材料进行建模和分析。

本文将围绕C++语言下的二维四节点四边形等参元刚度矩阵展开讨论。

1. 了解四节点四边形等参元在开始讨论C++下的四节点四边形等参元刚度矩阵之前，首先应该对四节点四边形等参元有所了解。

四节点四边形等参元是指在二维空间中，以四个节点构成的四边形元素，同时该元素的几何形状和内部分布的自由度均由相同的形函数进行描述的元素。

在有限元分析中，等参元的应用可以大大简化模型建立的复杂度，并提高计算的精度。

2. 实现C++中的四节点四边形等参元在C++语言中，实现四节点四边形等参元需要考虑节点的坐标、材料属性、边界条件等因素，并结合有限元理论中的导出公式进行编码。

通过C++语言的面向对象特性，可以将节点、单元、材料等抽象为对象，以便更好地管理和组织相关数据和函数。

3. 深入分析四节点四边形等参元的刚度矩阵四节点四边形等参元的刚度矩阵是描述该元素对外加载的响应性能，是有限元分析中至关重要的一部分。

通过推导理论公式，并结合C++语言进行具体实现，我们可以得到该元素的刚度矩阵。

在这一过程中，需要考虑矩阵装配的方法、高效的数据结构、数值计算的稳定性等因素。

4. 个人观点和总结通过对C++下的二维四节点四边形等参元刚度矩阵的探讨，我深刻地意识到了有限元分析与计算力学领域的复杂性和重要性。

C++语言作为一种高效、灵活的编程语言，为工程领域的数值计算提供了强大的支持。

对于工程师和研究人员来说，深入理解和掌握有限元分析的原理和实现方法，将有助于提升对复杂结构和材料行为的认识和预测能力。

对于C++语言下的四节点四边形等参元刚度矩阵，需要我们不仅要掌握有限元分析的理论知识，还需要熟练掌握C++语言的编程技巧和工程应用。

平面四边形4结点等参有限单元法

平⾯四边形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――计算形函数对整体坐标的导数TiNxN，=i1,2,3,4。

FUN8――计算形函数及雅可⽐矩阵[J]SFUN ――应⼒磨平-单元下的‘K’=NCN‘SCN――应⼒磨平-单元下的右端项系数‘CN‘SUMSKN――应⼒磨平-单元下的右端项集成到总体的‘P‘SUMSTRS――应⼒磨平-单元下的集成到总体的‘K‘GAUSTRSS――⾼斯消元求磨平后的应⼒3、输⼊数据及变量说明当程序开始运⾏时，按屏幕提⽰，键⼊数据⽂件的名字。

在运⾏程序之前，根据程序中INPUT需要的数据输⼊建⽴⼀个存放原始数据的⽂件，这个⽂件的名字为INDAT.DAT。

四边形四节点等参元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。

数学结构有限元法

x1
y1
x
2
x y
N1 0
0 N1
N2 0
0 N2
N3 0
0 N3
N4 0
0 N4
y2 x3ຫໍສະໝຸດ y3x4 y4
N1
1 4
1 1
N2
1 4
1 1
N3
1 4
1 1
N4
1 4
1 1
用结点的坐标值 x1, y1, x2 , y2 , x3 , y3 , x4 , y4 插值表示出单元内的坐标
x, y
与单元分析中常用的结点位移插值一样， N1, N2 , N3, N4 也可称为形状函数，Ni , 称为几何形状函数。
1, 1 代入上式，则得 x x2 , y y2 变成图(b)中相应线
两个单元的等百分线也一一对应
1 的直线，通过式(6-1)变换之后即是 x y 平面上2－3直线
2．等参单元的收敛性
§6.2 八结点曲边等参单元
N1
1 4
1 1
1
8
8
u Ni ,ui , v Ni ,vi
i1
i1
N5
1 2
1 2
1
8
8
x Ni ,xi , y Ni ,yi
i1
i1
每一条边都是一条二次曲线。如令 1得
x
x2 2
1
x3 2
1
第六章等参数单元
§6.1 平面四结点等参元 6.1.1 坐标变换与等参单元
实单元
母单元
随单元形状而不同的局部坐标系，称为单元的自然坐标系
正方形的四个边对应于实际单元的边界，四个顶点也一一对应于四个结点；正方形内任一点 P(,)都对应于实际单元内的一个点

有限元课件第4讲等参元和高斯积分

关于坐标系
直角坐标系( x , y , z)
极坐标(r,) ，2维球坐标系(r,θ, ) 柱坐标系 (, , z)
自然坐标系
自然坐标系:
➢选轨迹上任一点O为原点 ➢用轨迹长度S 描写质点位置
m
OS
n
➢质点沿切线前进方向的单位矢量为切向单位矢量(tangential unit vector)
➢质点与切向正交且指向轨迹曲线凹侧的单位矢量为法向单位矢量(normal unit vector)
U e 1 (x( )) (x( ))dV 1 x2 Ee (x( )) (x( ))Aedx
2 e
2 x1
U e 1 1 EeB( )qeB( )qe Ae (le / 2)d
2 1
U e
1 qeT [
1
(l e
/
T
2)B
( )Ee AeB( )d ]qe
2
1
U e 1 qeT Keqe 2
x(,) N(,)xe
u(x(,), y(,)) u(,) N(,)qe
N1
1 4
(1
)(1 )
N2
1 4
(1
)(1 )
N3
1 4
(1
)(1 )
N4
1 4
(1 )(1)
ε(x(,), y(,)) u(,) N(,)qe B(,)qe
x
B(
, )
0
y
0
x 1 2 3 4 N1x1 N2x2 N3x3 N4x4
y
1
2
3
4
N1 y1
N2
y2
N3
y3
N4
y4
N1

平面四节点等参单元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软件计算结果计算结果分别罗列圆周方向的单元和半径方向的单元位移和应力。

平面四节点等参单元分析程序

平面四节点等参单元分析程序(总27页)本页仅作为文档页封面，使用时可以删除This document is for reference only-rar21year.March变分原理与有限元大作业平面四节点等参单元分析程序姓名：潘清学号：SQ10018014033完成时间：2011-4-26一、概述通常情况下的有限元分析过程是运用可视化分析软件(如ANSYS、ABAQUS、SAP等)进行前处理和后处理，而中间的计算部分一般采用自己编制的程序来运算。

具有较强数值计算和处理能力的Fortran语言是传统有限元计算的首选语言。

随着有限元技术的逐步成熟,它被应用在越来越复杂的问题处理中，但在实际应用中也暴露出一些问题。

有时网格离散化的区域较大,而又限于研究精度的要求，使得划分的网格数目极其庞大，结点数可多达数万个，从而造成计算中要运算的数据量巨大，程序运行的时间较长的弊端，这就延长了问题解决的时间，使得求解效率降低。

因为运行周期长，不利于程序的调试，特别是对于要计算多种运行工况时的情况；同时大数据量处理对计算机的内存和CPU 提出了更高的要求，而在实际应用中，单靠计算机硬件水平的提高来解决问题的能力是有限的。

因此，必须寻找新的编程语言。

随着有限元前后处理的不断发展和完善，以及大型工程分析软件对有限元接口的要求，有限元分析程序不应只满足解题功能，它还应满足软件工程所要求的结构化程序设计条件，能够对存储进行动态分配，以充分利用计算机资源，它还应很容易地与其它软件如CAD 的实体造型，优化设计等接口。

现在可编写工程应用软件的计算机语言较多，其中C语言是一个较为优秀的语言，很容易满足现在有限元分析程序编程的要求。

C语言最初是为操作系统、编译器以及文字处理等编程而发明的。

随着不断完善，它已应用到其它领域，包括工程应用软件的编程。

近年来，C语言已经成为计算机领域最普及的一个编程语言，几乎世界上所有的计算机都装有C 的编译器，从PC机到巨型机到超巨型的并行机，C与所有的硬件和操作系统联系在一起。

第4章平面问题的等数参单元√ 有限元

2
2
变换到母单元上的单元刚度矩阵为
[ K ]e h [ B]T [ D][ B] J dd
1 1
(10' )
这个式子中包含雅可比矩阵的逆矩阵，使上式的解析积分相当困难，因此，通常用数值积分法进行积分。
（见沈永欢：实用数学手册，科学出版社；王勖成：有限单元法，清华大学出版社）
４单元的等效荷载体积力的等效荷载考虑单元内任一点的体积力为
i 1,2,3,4 (11) j 5,6,7,8

实际单元中任一点坐标可用下式表示
x N i ( , ) xi，y N i ( , ) yi
i 1 i 1
8
8
２单元刚度矩阵几何方程为 u x x x v { } y 0 y xy v u x y y 其中
u 1 2 3 4 2 5 6 2 7 2 8 2 2 2 2 2 v 9 10 11 12 13 14 15 16
(8)
应变可写为
{ } [ B]{ }e
其中应变矩阵为
(6' )
N r x 0 N r y N s x 0 N s y 0 N s y N s x
N i 0 x x [ B] 0 [N ] 0 y N i y x y
1 形函数与坐标变换在母单元上，取位移函数为 u 1 2 3 4 v 5 6 7 y 8
按照上一章求插值函数的方法，得

平面四边形四节点等参单元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 *********************************************。

有限元讲义等参单元

显然该位移模式在ξ,η坐标系下是双线性位移模式，在x，y坐标系下不是双线性位移模式。容易验证，上述形函数满足形函数性质。
• 等参单元的提出对于有限元法在工程实践中的应用具有重要意义。
§6.2 平面四节点等参单元 1、局部坐标系与位移模式
下图为一个4节点任意四边形单元，单元有8个自由度。将矩形单元放松为4节点任意四边形单元将带来许多好处。
§6.1 引言
• 解决上述矛盾的出路就是突破矩形单元和六面体单元几何上的限
下图为一个4节点制任，意四使边形其单成元，为单元平有面8个自任由意度。四边形和空间任意六面体单元，如果再增
2 平面四节点等参单元
由于等参单元刚加度矩边阵中积分间式中节被积点函，数很还难导可出以解析成表达为式曲，因边此等四参边单元形的计和算曲都采面用数六值面积分体求积高分精的近度似值单，元有限。元中对四
有限元讲义等参单元
§6.1 引言
• 回顾前面的各种二、三维单元，这些单元受到两个方面的限制：
➢ 第一是单元的精度，显然单元的节点数越多，单元精度越高。因此在这一点上，矩形单元优于3节点三角形单元，六面体单元优于四面体单元；
➢ 第二是单元几何上的限制，矩形和六面体（长方体）单元要求单元的边（面）平行于坐标轴（面），因此都不能模拟任意形状和方位的结构。此外，线性单元都是直线边界，处理曲边界几何体误差较大。
• 建立了局部坐标系后，在ξ-η平面内单元就是一个边长为2的正方形。
§6.2 平面四节点等参单元
• 该局部坐标系使得在x-y平面上的任意四边形与ξ-η平面上的正方形之间形成了1-1对应的映射。正方形的4个顶点对应任意四边形单元的四个节点； 4条边对应任意四边形单元的4条边；正方形内任一点p(ξ,η)对应于任意四边形内一点p(x,y)。

有限元分析第5章平面矩形单元与平面等参单元

f

1 1

x3 y3

f

1 1

x3 y3

f

1

1

求出待定系数，得
xN1x1N2x2N3x3N4x4 yN1y1N2y2N3y3N4y4
其中： N 1
1 1 1
r,si, j,m,k
平面矩形单元小结
1. 优点：矩形单元的应力、应变为一次线性函数，精度要比三角形三节点高；
2. 不足：实际问题很难用4节点矩形单元划分，特别是边界适合性不强；
3. 问题：能否构造一种任意四边形单元，则在提高精度得前提下，边界适应性还强？
等参单元
二、平面等参单元
等参单元（iso-parametric element）的概念
um
vi i (xi , yi)
vj
e
uj
ui
j (xj , yj)
一、平面矩形单元
b
vk uk
k (xk, yk)
vi
(xi, yi)
i
u
i
a
y, x,
a
vm
m
um
(xm, ym)
vj j
(x , y )
jj
uj
b
矩形单元也是常用的单元之一，由于采用了比常应变三角形单元更高次数的位移模式，故可以更好地反映弹性体的位移状态和应力状态。

J
1

N i

y

x

J

一种基于张量列式的高性能平面4节点杂交_混合单元

文章编号:0258 2724(2000)01 0072 05一种基于张量列式的高性能平面4节点杂交/混合单元杨帆, 王王君, 陈大鹏(西南交通大学工程科学研究院,四川成都610031)摘要:采用张量形式,从Hellinger Reissner 变分原理出发,建立一种列式杂交/混合有限单元的有效方法。

该方法利用了张量的不变性,剔除了坐标系转换的影响,通过应力参数与独立变形模式一一对应的方法,选择并优化应力场,由此建立高性能的杂交/混合单元。

据此方法,建立了一个4节点杂交/混合单元。

数值算例表明,该单元具有不变性、普适性、不含零能模式,性能优良。

关键词:张量;有限元法;杂交/混合有限元中图分类号:O242.21 文献标识码:AAn Efficient 4 Node Plane Hybrid/Mixed Element with Tensorial FormulationYANG Fan, WANG Jun , C HE N Da peng(School of Eng.Science,Southwest Jiaotong University,Chengdu 610031,China)Abstract:The paper is mainly concerned with developing a kind of tensorial formulation scheme for loworder hybrid/mixed finite elements based on Hellinger Reisser variational principle.By making use of the invariance properties of tensors,the disadvantages due to coordinates transfor mation are totally avoided.Also,in matching the derived strain field,whic h is in tensorial form,a proper selection of the optimal stress field is successfully achieved on the basis of the consistency of one stress parameter with one independent deformation mode .Thus,a 4 node plane quadrilateral ele ment is obtained.Numerical simulations reveal the re markable charac teristics of the present elements,such as invariance,robustness,having no zero energy,etc .Key words:tensors;finite ele ment methods;hybrid/mixed finite ele ment有限元法已广泛应用于实际工程问题的计算分析。

有限元四面体及六面体单元

空间问题有限基本元概念分析
7.单元刚度矩阵
4 节点四面体
空间问题有限基本元概念分析
7.单元刚度矩阵
4 节点四面体
空间问题有限基本元概念分析
7.单元刚度矩阵
4 节点四面体
空间问题有限基本元概念分析
8. 4节点四面体单元的常系数应变和应力
4 节点四面体
空间问题有限基本元概念分析
节
位移场可以设定4个待定系数，根据节点个数以及确定位移模式的基本原则（从低阶到高阶的完备性、唯一确定性），选取该单元的位移模式为
点
四
面
(4-104)
体
(4-105)
空间问题有限基本元概念分析
2.单元位移场的表达
4
将式(9-3)代入节点条件(9-4)中，可求取待定系数(ai，bi，ci)，i＝0，1
节
取该单元的位移模式为
点
正
六
(4-117)
面
体
(4-118)
空间问题有限基本元概念分析
3.其它物理参量的表达
8
在得到该单元的形状函数矩阵后，就可以按照有限元分析的标准过程推导相应的几何矩阵、刚度矩阵、节点等效载荷矩阵以及刚度方
节
程，相关情况如下
点
正
(4-119)
六
面
体
(4-120)
(4-121)
1. 单元的几何和节点描述
8
该单元为由8节点组成的正六面体单元(hexahedron element)，每个节点有3 个位移(即3个自由度)，单元的节点及节点位移如图所示
节
点
正
六
面
体
(4-115)

1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。

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)
0th GQ4 QA[4] QA2[4] QR4a[4] QR4b[4] QA4[4] QM6[4] QE2[4] Q8[4] 1st GQ4 2nd GQ4 最佳解[4]
11.73 20.89 21.70 21.27 22.48 23.29 21.05 21.35 23.02 22.77 23.46 23.9
-0.1525 -0.1944 -0.1974 -0.1945 -0.2003 -0.1888 -0.1854 -0.1859 -0.1863 -0.1927 -0.2050 -0.2010
9.3.2
变
计算精度考查
悬臂梁
但只引起微小的应计算网格本算例常用于考查单元的精度
悬臂梁
由于载荷容易引起端部附近单元大幅度的刚体运动计算结果列于表 9-2
0.1640 0.1798 0.1931 0.1840 0.2074 0.1883 0.1773 0.1956 0.2231 0.2327 0.2358 0.2360
-0.1723 -0.1746 -0.1868 -0.1780 -0.2129 -0.1894 -0.1666 -0.1448 -0.1855 -0.2200 -0.1933 -0.2010
0.0926 0.9029 0.9769 0.9926 0.9713 0.9815 1.0000
0.0278 0.7632 0.8844 0.9149 0.7882 0.9630 1.0000
0.0370 0.8409 0.9195 0.9602 0.8233 0.9722 1.0000
9.3.4
表 9-4
Table 9-3 Comparison of vertical displacements at the free end of MacNeal thin beam 单元类型 0 (Q4) QA[4] QA2[4] QA4[4] Q8[5] 1stGQ4 2ndGQ4
th
纯弯
规则梯形平行四边形规则
245
大连理工大学博士学位论文
v. vi.
累加各高斯点处的广义节点单刚阵转到 i 直到循环结束
K ij = ∑ K ijg
g
(b) 将 K ij 组装到总体刚度阵中的相应位置 (3) 转到步骤 (2) 直到循环结束 (4) 完成单元刚度阵的形成与组装
9.3
能
数值算例
从不同角度考查 GQ4 的性
选用文献[4]中用于考查四边形单元的一些典型例题
其中分块子阵 K ij (i, j = 1,2,3,4) 为
K 12 K 22 K 32 K 42
K 13 K 23 K 33 K 42
K 14 K 24 K 34 K 44
(9-8)
K ij = ∫ Bei ( E ) Bej tdxdy
A
( ) ( )
(9-9)
K ij 称为广义节点刚度阵其阶数与所采用的广义节点插值函数的阶次有关对于一
(a) 规则网格
248
第九章
平面广义四节点等参元 GQ4
45
(b) 不规则网格 1
45
(c) 不规则网格 2 图 9-3 Fig.9-3 MacNeal 问题计算网格
Regular and irregular Meshes for MacNeal thin beam
表 9-3
MacNeal 细长梁自由端竖向位移计算结果比较
i, j 分别遍历单元的四个节点
(a) 按高斯点循环在局部坐标下即计算面内形成该高斯点处的传统四节点等参元形函数阵 ϕ e ii. iii. iv. 求高斯点的自然坐标 ( x i , y i ) 在自然坐标系下求广义节点插值函数在高斯点处的取值及导数值按式(9-5)或(9-6)形成高斯点处的几何矩阵 Bij 高斯点处的广义节点单刚 K ijg 并按式(9-10)形成该
∂ ( y 2ϕ ei ) L ∂x
式
9-6 分别为一阶广义节点与二阶广义节点所对应的几何阵形式 9-5 9-6 可以看出在广义四节点等参元中几何阵各分量的求导可 y 的导数另一一部分是传统四节点等参元的形函数 ϕ e 对总体坐标 x y 的导数而且注意到进行即计算面内
部分则是广义节点插值函数对总体坐标 x 数对坐标的导数的计算无需在自然坐标系下
∂ ( yϕ ei ) ∂x 0 ∂ ( yϕ ei ) ∂y
∂ ( yϕ ei ) (i = 1,2,3,4) (9-5) ∂y ∂ ( yϕ ei ) ∂x 0 ∂ ( y 2ϕ ei ) (i = 1,2,3,4) ∂y ∂ ( y 2ϕ ei ) ∂x 0
(9-6)
∂ϕ ei ∂x B ei = 0 ∂ϕ ei ∂y
9-5 从式分为两部分
0 ∂ϕ ei ∂y ∂ϕ ei ∂x
∂ ( xϕ ei ) ∂x 0 ∂ ( xϕ ei ) ∂y
0
∂ ( xϕ ei ) L 0 ∂y ∂ ( xϕ ei ) ∂ ( y 2ϕ ei ) L ∂x ∂y
9.3.1
材料参数
计算精度考查
取弹性模量为 E = 1.0
Cook 变截面斜悬臂梁
计算简图见图 9-1 作用于自由端泊松比为 ν = 1 / 3 厚度为 t = 1
该问题用于考查单元在实际不规则网格情形下的计算精度的荷载 P 为单位荷载计算结果见表 9-1
y
B
16
1
y
B 16
1
A x 48
A x 48
因此是一种特别容易引起误差的结构
如图 9-2 共划分了 5 种不同网格
梁挠度精确解为 4.03
(100,10) P (0,0)
(a)
P
(b)
(75,10) P (25,0)
(c)
247
大连理工大学博士学位论文
P
(d)
(50,10) (16.67,0) (50,0) (83.33,10) P
(e) 图 9-2 悬臂梁单元划分 Fig. 9-2 FEM Meshes and load cases of cantilever
(a) 2 × 2 网格
(b) 4 × 4 网格
图 9-1 变载面斜悬臂梁弯曲计算网格 Fig. 9 FEM meshes for a beam with variable section
246
第九章
平面广义四节点等参元 GQ4
表 9-1 变截面斜悬臂梁弯曲计算结果比较 Table 9-1 Calculated displacements at point A and stresses at point B and C of the beam 单元
广义节点插值函
直接求导便可
9.2.3
Байду номын сангаас
单元刚度阵
用 E 表示弹性矩阵
T A
考虑线弹性问题
则单元刚度的表达式为 (9-7)
K e = ∫ (B e ) ( E )(B e ) tdxdy
具体地有
244
第九章
平面广义四节点等参元 GQ4
K 11 K K e = 21 K 31 K 41
端部受剪
梯形平行四边形
0.0926 0.9100 0.9837 1.0 0.9641 0.9926 0.9963
0.0222 0.7693 0.8948 0.9222 0.7722 0.9926 0.9963
0.0296 0.8793 0.9556 0.9852 0.8811 0.9926 0.9963
单元的畸形敏度分析
表 9-5 分别为纯弯载荷和梁端部集中力作用下梁端点的挠度值与精确值的比弹性模量为 E = 1 500 泊松比为
考查图 9-4 所示的悬臂梁在纯弯载荷和端部集中力载荷作用下单元的畸形敏度值材料参数选为
ν = 0.25 载荷如图中所示
249
大连理工大学博士学位论文