有限元分析半开卷资料(共10页).docx

精选优质文档-倾情为你奉上有限元分析半开卷资料基本概念：位移a=u1v1u2v2unvnT；等效结点荷载R=R1xR1yR2xR2yRnxRnyTu=uv=Nae; Fe=kae; =xyxyT=Lu=LNae=Bae; =xyxyT=D=DBae=Sae第二章 平面弹性力学问题1.位移模式与收敛性条件三角形单元的位移模式：u=Niui+Njuj+Nmumv=Nivi+Njvj+Nmvm式中Ni=ai+bix+ciy2A i,j,m；ai=xjym-xmyj;bi=yj-ym; ci=-(xj-xm), A为单元面积A=121xiyi1xjyj1xmym (为了使面积A不成为负值，规定结点I, j, m的次序按逆时针转向)矩形单元的位移模式：u=Niui+Njuj+Nmum+Npupv=Nivi+Njvj+Nmvm+Npvp式中Ni=141+ixa1+iyb; i=xixi; i=yiyi (i,j,m,p)位移模式需满足的条件：a 位移模式必须能反映单元的刚体位移（与本单元的变形无关的位移） b. 位移模式必须能反映单元的常量应变（与位置坐标无关的应变） c. 位移模式应当尽可能反映位移的连续性（相邻单元之间位移的连续性和单元内部位移也是连续的）以三角形单元为例说明，a1,a4,a5-a32反映了刚体移动和刚体转动，a2,a6,a3+a5反映常量应变2.形函数及其性质形函数：即插值基函数，反应单元的位移形态，因而也称为位移的形态函数，简称形函数(NiNjNm)形函数的性质：（2） 在节点i上Ni=1，在其他节点上Ni=0，该性质可导出形函数在三角形单元上的积分和在某边界上的积分为eNidxdy=13A,ijNids=12lije表示对单元积分，ij表示对单元的ij边线积分（3） 在单元中，任意点形函数之和等于1（4） 形函数的值在0-1间变化3.面积坐标及三角形高次形函数的构造面积坐标：三角形单元中，任一点P(x,y)与其3个角点相连形成3个子三角形，其位置可以用三个比值来确定，即面积坐标。面积坐标只限于用在一个三角形单元内，因而是一种局部坐标。在平行于jm边的一根直线上的所有点，都具有相同的Li坐标，而且这个坐标就等于“该直线至jm边的距离”与“结点i至jm边的距离”的比值用直角坐标表示面积坐标的关系式：Li=(ai+bix+ciy)2A式中ai=xjym-xmyj,bi=yj-ym,ci=-xj+xm用面积坐标表示直角坐标的关系式：x=xiLi+xjLj+xmLmy=yiLi+yjLj+ymLm六结点三角形单元的位移模式：u=Niui+Njuj+Nmum+N1u1+N2u2+N3u3v=Nivi+Njvj+Nmvm+N1v1+N2v2+N3v3Ni=Li2Li-1 (i,j,m),N1=4LjLm (1,2,3)4.有限元支配方程的推导（结构力学法/变分原理）根据平衡条件，各环绕单元对该结点作用的结点力之和应等于由各环绕单元移置而来的结点荷载之和，即eFi=eRi，再将结点力公式代入，变可得有限元的支配方程Ka=R,式中K为整体刚度矩阵，a为整体结点位移列阵，R为整体结点荷载列阵变分原理导出有限元支配方程证明：用最小势能原理推导出有限元的求解方程=12TDtdxdy-uTftdxdy-SuTfds式中，t是平面弹性体的厚度，f是体积力，f是物体表面的面力离散成有限网格，其中=Bae, u=Nae代入上式可得=12e(ae)TeBTDBtdxdyae-e(ae)TeNTftdxdy-e(ae)TseNTftds利用ae=Cea, 则=12aTKa-aTR, 其中K=eCeTkCe;R=eCeTRe;k=eBTDBtdxdy; Re=eNTftdxdy+seNTftds根据最小势能原理=0，即a=0, 这样就得到有限元的求解方程Ka=R虚功原理建立有限元的支配方程假设单元发生了虚位移，其相应虚位移为u=uvT而该单元上各结点的相应虚位移为ae=uiviujvjumvmT按照静力等效原理，即结点荷载与原荷载在上述虚位移上的虚功相等，有(ae)TRe=uTP将u=Nae 代入，得(ae)TRe=（ae）TNTP由于虚位移是任意的，得Re=NTP5.荷载列阵：单元到整体体力引起的等效结点荷载：Re=-bb-aaNTftdxdy面力引起的等效结点荷载：Re=-bbNTftdy常见分布荷载产生的等效结点荷载1.单元自重 Re=-13gtA010 101T 2.在ij边界上受x方向均布力q作用，边界长度为l，Re=12qlt101 00 0T 3.在ij边界上受三角形分布荷载作用，边界长度为l，Re=12ql23013 00 0T 整体结点荷载列阵：确定每个单元的结点荷载列阵，然后根据各个单元的结点局部编码与整体编码的对应关系，将单元的结点荷载列阵中每个子矩阵叠加到R的相应位置上。6.刚度矩阵：单元到整体单元刚度矩阵（平面应变问题）k=-bb-aaBTDBtdxdy线性位移模式下，有k=BTDBt=kiikijkimkjikjjkjmkmikmjkmmFe=UiViUjVjUmVm=kiixxkiixykiiyxkiiyyuiviujvjumvm kij表示j结点对i结点的刚度贡献kiixx+kjixx+kmixx=0;kiiyx+kjiyx+kmiyx=0单元刚度矩阵的力学意义：单元刚度矩阵中任一个元素（如kijyx）表示当j结点x方向发生单位位移时，在i结点y方向产生的结点力单元刚度矩阵的性质：对称性，奇异性，主元素恒正，单元均匀放大或缩小不会改变刚度矩阵的数值，单元水平或竖向移动不会改变刚度矩阵的数值平面应力问题 krs=Et41-v2Abrbs+1-v2CrCsvbrbs+1-v2CrCsvCrCs+1-v2brbsCrCs+1+v2brbs (r=I,j,m s=I,j,m)平面应变问题krs=E1-vt41+v1-2vAbrbs+1-2v2(1-v)crcsv1-vbrcs+1-2v2(1-v)crbsv1-vcrbs+1-2v2(1-v)brcscrcs+1-2v2(1-v)brbs整体单元建立步骤：将K全部充零，逐个单元地建立单元 的刚度矩阵，然后根据单元结点的局部编码与整体编码的关系，将单元的刚度矩阵中每一个子矩阵叠加到K中的相应位置上。对所有的单元全部完成上述叠加步骤，就形成了整体刚度矩阵。7.简单问题的有限元具体计算（利用已算好的刚度矩阵）8.计算结果的整理与分析计算成果主要包括两个方面：即位方面和应力方面。主要进行对应力方面的计算成果讨论。把计算出来的常量应力作为单元形心处的应力，拉应力用箭头表示，压应力用平头表示。由计算成果推出结构物内某一点的接近实际的应力，必须通过某种平均计算，通常可采用绕结点平均法或二单元平均法。绕结点平均法：把环绕某一结点的各单元中的常量应力加以平均，用来表征该结点处的应力。（在内结点处具有较好的表征性，但在边界结点处表征性较差，因此边界结点处的应力宜用内结点应力插值推算）二单元平均法：把两个相邻单元中的常量应力加以平均用来表征公共边中点处的应力。PS：只容许对厚度及弹性常数都相同的单元进行平均计算9.网格划分的注意事项a. 应根据所需精度，在合理的计算时间内来决定单元的大小 b.为了便于整理和分析应力成果，往往采用直角三角形的单元 c.应合理利用结构的对称性，减少计算量 d.为了使地基弹性对结构物中应力的影响能反映出来，必须把和结构物相连的那一部分地基也取为弹性体，和结构物一起作为计算对象e. 可进行二次计算在以下情况下应将单元尺寸减小：a.计算对象的厚度或其弹性有突变之处，且应吧突变线作为单元的界线，b.计算对象受季度突变的分布荷载或受集中荷载 c.结构物具有凹槽或孔洞第三章 平面等参有限元1. 等参单元的概念将位移模式与坐标变换式具有相同的形式，即形函数相同，参数个数相同的单元称为等参单元。2.单元形函数的构造方法整体坐标系为（x,y），在每个单元上建立局部坐标系(,),只需雅克比行列式大于0即可保证坐标一一对应的实现。J=xyxy>0坐标转换式：x=N1x1+N2x2+N3x3+N4x4 y=N1y+N2y2+N3y3+N4y4Ni=14(1+i)(1+i)3.位移模式及收敛性要求单元位移模式：u=N1u1+N2u2+N3u3+N4u4 v=N1v1+N2v2+N3v3+N4v4收敛性要求：完备性和连续性以4结点四边形等参单元为例说明完备性，i=14Ni=1，i=14Nixi=x，i=14Niyi=yi=14Ni=14(1-)1-+141+1-+141+1+141-1+=121-+121+=1对于具有m各结点的等参单元，为了位移模式满足完备性要求，形函数必须满足i=1mNi=1连续性：只需考察任意两单元交界面上位移是否连续即可。为了保证整体坐标与局部坐标的一一对应关系，单元不能歪斜，单元的各边长不能等于0。4.单元刚度矩阵和荷载列阵的计算Bi=Nix00NiyNiyNix; Si=E1-v2NixvNiyvNixNiy1-v2Niy1-v2Nix 需根据复合函数的求导法则来求出各形函数对整体坐标的导数Ni=Nixx+Niyy; Ni=Nixx+Niyy 即可解得Nix和Niy单元刚度矩阵：k=eBTDBtdA=-11-11BTDBtJdd体力引起的单元结点荷载Re=eNTftdA=-11-11NTftJdd面力引起的单元结点荷载，在=±1；Re=-11Ntft(x)2+(y)2d;在=±1；Re=-11Ntft(x)2+(y)2d以平行四边形为例，常见荷载下的单元结点荷载a. 在自重f=0-g;Re=-14gtATb. 在=1;Re=12atqTc. 在=1上法向分布面力的集度为q, Re=-11NTy-xqtd5.高斯数值积分数值积分一般有两类方法，一类是等间距数值积分，如辛普森方法等。另一类是不等间距数值积分，如高斯数值积分法。高斯积分法对积分点的位置进行了优化处理，所以精度较高。例：求两点高斯积分的积分点坐标和积分权系数二次多项式：P=(-1)(-2)积分点位置由下式确定-11iPd=0 i=0,1； P=-1-2*(-n)当i=0时，-11-1-2d=23+212=0当i=1时，-11-1-2d=-231+2=0联立方程后可解得1=-0.577；2=0.577积分权系数为：Hi=-11lin-1()d; lin-1=-1-2*-i-1-i+1*(-n)i-1i-2*i-i-1i-i+1*(i-n)可解得H1=H2=1第四章 空间弹性力学问题1.单元形函数的构造方法V为四面体的体积V=161xiyizi1xjyjzj11xmxpymypzmzp 式中 ai=xjyjzjxmymzmxpypzp; bi=-1yjzj1ymzm1ypzp; ci=-xj1zjxm1zmxp1zp; ai=-xjyj1xmym1xpyp12.位移模式及收敛性要求位移模式：u=Niui+Njuj+Nmum+Npup ; v=Nivi+Njvj+Nmvm+Npvpw=Niwi+Njwj+Nmwm+Npwp; Ni=16Vai+bix+ciy+diz (i,j,m,p)3.单元刚度矩阵和荷载列阵的计算krs=E(1-v)361+v1-2vV*brbs+A2(crcs+drds)A1brcs+A2crbsA1brbs+A2drbsA1crbs+A2brcscrcs+A2(brbs+drds)A1crds+A2drcsA1drbs+A2brdsA1drcs+A2crdsdrds+A2(crcs+brbs)若单元受自重作用，体积力f=00-gT;Re=-14gV1T,表面单元重量平均分到4个结点。若单元某边界面如ijm面x方向受线性分力作用，设结点i的面力集度为q，结点就，p集度为0，则面力矢量f=Niq00T，Re=13qAijm0000T4.高斯数值积分第六章 有限元程序设计1.主元素序号指示矩阵MA由子程序CBAND来实现（1） 将MA（*）全部充零（2） 对全部单元循环，利用结构自由度序号矩阵JR（2，*），计算每个单元的结点自由度序号数组NN（8）单元i，NNx=单元ii，NNx=单元iii，NNx=单元iv，NNx=（3） 计算单元自由度序号中最小的序号，赋给L（4） 计算与单元各自由度对应行的半带宽。把各单元计算得到的半带宽中的最大值赋给MA（*）数组, MA*=9（5） 对数组MA（*）逐个元素累加，便得到所需的主元素序号指示相量，此时MA*=2.半带宽的计算K中每行从第一个非零元素到对角线元素的元素个数称为该行的半带宽。半带宽是衡量结点编号好坏的一个指标，最优的结点编号可以使半带宽最小半带宽等于该行行号（即自由度序号）减去与该自由度号关联的最小自由度号再加1，即JPL=JP-L+1,JPL为半带宽，JP为自由度序号，L为某一单元最小自由度序号。3.一维变带宽的存储将按一维变带宽存储的K所需的存储量记为NH，存储有效元素的一维数组记为SK（NH），结构总自由度数，亦即需要建立的方程个数记为N。则对于本例来讲，N=12，NH=66。通过指示向量MA（*）便可将一维数组SK（*）与二维矩阵K的元素一一对应起来一维数组与二维矩阵K有如下对应关系（1） K中第I行的主元素就是SKMA(I)（2） KI,J=SKMAI-(I-J)4.结点自由度序号矩阵JRJR2,12=5.体力引起的荷载列阵6.面力引起的荷载列阵7.荷载列阵的集成8.刚度矩阵的集成由子程序SKO来实现（1） 采用3个积分点的高斯数值积分来计算单元刚度矩阵，对高斯点局部坐标和权系数赋值（2） 在对单元循环的循环体内形成单元自由度数组NN（8）和单元结点坐标值XYZ（2,4）。NN（8）是为了将单元刚度矩阵的元素累加到整体刚度矩阵的需要，XYZ（2,4）是为了计算单元刚度矩阵的需要。（3） 调用子程序STIF计算单元的刚度矩阵SKE（8,8）（4） 将每个单元的刚度矩阵累加到一维刚度矩阵数组SK（*）中去。9.读懂源程序10在源程序基础上修改某些功能三角形单元有关公式1. 位移模式（三角形单元）：u=Niui+Njuj+Nmumv=Nivi+Njvj+Nmvm式中Ni=ai+bix+ciy2A i,j,m；ai=xjym-xmyj;bi=yj-ym; ci=-(xj-xm), A为单元面积A=121xiyi1xjyj1xmym (为了使面积A不成为负值，规定结点I, j, m的次序按逆时针转向)2. 应变转换矩阵B=12Abi0bj0bm00cicibi0cjcjbj0cmcmbm 3. 应力转换矩阵平面应变问题Si=E2-v2Abivcivbici1-2v2(1-v)ci1-2v2(1-v)bi 对于平面应力问题，要把E换成E1-v2, v换成v1-v;Si=E(1-v)21+v1-2vAbiv1-vciv1-vbici1-2v2(1-v)ci1-2v2(1-v)bi4. 单元刚度矩阵平面应力问题 krs=Et41-v2Abrbs+1-v2CrCsvbrbs+1-v2CrCsvCrCs+1-v2brbsCrCs+1+v2brbs (r=I,j,m s=I,j,m)平面应变问题krs=E1-vt41+v1-2vAbrbs+1-2v2(1-v)crcsv1-vbrcs+1-2v2(1-v)crbsv1-vcrbs+1-2v2(1-v)brcscrcs+1-2v2(1-v)brbs5. 等效结点荷载体力引起的等效结点荷载：Re=-bb-aaNTftdxdy面力引起的等效结点荷载：Re=-bbNTftdy常见分布荷载产生的等效结点荷载1.单元自重 Re=-13gtA010 101T 2.在ij边界上受x方向均布力q作用，边界长度为l，Re=12qlt101 00 0T 3.在ij边界上受三角形分布荷载作用，边界长度为l，Re=12ql23013 00 0T 4结点四边形等参单元的有关主要公式1. 位移模式u=uv=Nae；其中N=N10N20N30N400N10N20N30N4;ae=u1v1u2v2u3v3u4v4T;Ni=14(1+i)(1+i) 2. 坐标变换x=N1x1+N2x2+N3x3+N4x4 y=N1y+N2y2+N3y3+N4y43. 应变公式=xyxyT=Bae 其中B=B1B2B3B4; Bi=Nix00NiyNiyNix; NixNiy=J-1NiNi4. 应力公式=D; x=D1x+D2y;y=D2x+D1y; xy=D3xy 其中D=D1D20D2D1000D3; D1=1-vE(1+v)(1-2v); D2=vE(1+v)(1-2v); D3=E2(1+v)5. 单元刚度矩阵k=-11-11BTDBJdd; k=k11k12k13k14k21k22k23k24k31k41k32k42k33k43k34k44kij=-11-11Nix0Niy0NiyNixD1D20D2D1000D3Nix00NiyNiyNixJdd 令每一子块矩阵为kij=kiikijkjikjj 其中kii=-11-11(D1NixNjx+D3NiyNjy)Jdd ;kij=-11-11(D2NixNjy+D3NiyNjx)Jddkji=-11-11(D2NiyNjx+D3NixNjy)Jdd; kjj=-11-11(D1NiyNjy+D3NixNjx)Jdd6. 单元等效结点荷载（1） 集中力P=PxPyT; Re=NTP; Rix=NiPx;Riy=NiPy （2） 自重体力f=0-gT; Re=-11-11NT0-gJdd ; Rix=0; Riy=-g-11-11NiJdd （3） 分布面力f=fxfyT 在=±1边界，Re=-11Nf(x)2+(y)2d在=±1边界，Re=-11Nf(x)2+(y)2d如果面力为法向压力，则f=-qlqm(i) 在=-1边界 l=-y(x)2+(y)2, m=x(x)2+(y)2 ; Re=-11NTy-xqd(ii) 在=1边界 l=y(x)2+(y)2, m=-x(x)2+(y)2; Re=-11NTy-xqd(iii) 在=-1边界l=y(x)2+(y)2, m=-x(x)2+(y)2；Re=-11NTy-xqd子程序的主要功能专心-专注-专业INPUT输入原始数据CBAND形成主元素序号指示矩阵MA（600）SKO形成整体刚度矩阵KCONCR计算集中力引起的等效结点荷载ReBODYR计算自重体力引起的等效结点荷载ReFACER计算分布面力引起的等效结点荷载ReDECOP支配方程Crout直接解法中的分解和前代计算FOBACrout直接解法中的回代计算OUTDISP输出结点位移分量STRESS计算单元应力分量OUTSTRE输出单元应力分量STIF计算单元刚度矩阵FDNX计算形函数对整体坐标的导数NixNiyTFUN8计算形函数及雅克比矩阵J程序中的主要变量及其意义NP结点总数NE单元总数NM材料类型数NR受约束的结点数N自由度总数MX最大半带宽NH整体劲度矩阵按一维存储的总容量NCP受集中力作用的结点的数目NBE要考虑自重作用的单元的数目IZ不同类型面力的批数程序中的主要数组及其意义，括号中的数字表示该数组所设定的最大容量SK(800)按一维存储的整体刚度矩阵COOR(2,300)结点坐标AE(4,11)材料常数，每种材料4各常数，它们分别为弹性模量、泊松比、容量和单元厚度MEL(5，200)单元信息，每个单元需要5个数据，前4个为该单元的4各结点号码，第5个为该单元的材料类型号WG(4)存放面力信息，对于某种类型的面力，这4个数分别为Ir, ,Z0,NSUIr=1受水压力作用Ir=2受均布法向压力作用 水的容重Z0最高水位所对应的y坐标值（y轴向上）NSU该批受面力作用的单元的单元号JR(2，300)结点自由度序号矩阵MA(600)主元素序号指示矩阵R(600)开始存放等效结点荷载，求解方程以后存放结点位移Iew(30)存放受面力作用的单元的单元号STRE(3，200)单元应力分量RF(8)单元等效结点荷载SKE(8,8)单元刚度矩阵NN(8)单元结点自由度序号矩阵RSTG(3)高斯积分点的坐标值H(3)高斯积分点的权系数对于3个积分点的数值积分格式，它们为RSTG(1)=-0.RSTG(2)=0RSTG(3)=0.H(1)=0.H(2)=0.8888H(3)=0.对于2个积分点的数值积分格式，它们为RSTG(1)=-0.RSTG(2)=0.57750H(1)=1H(2)=1FUN(4)4个函数PN(2,4)形函数对局部坐标，的偏导数DNX(2,4)形函数对整体坐标x,y的偏导数XJAC(2,2)雅克比矩阵J