第九章 平面广义四节点等参元 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
最新平面四边形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。
- 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
大连理工大学博士学位论文