平面三角形单元有限元程序设计

1、平面三角形单元有限元程序设计1如图1所示，一个厚度均匀的三角形薄板，在顶点作用沿板厚方向均匀分布的竖向载荷。已知：P=150N/m,E=200GPa,p=0.25,t=0.1m,忽略自重。试计算薄板的位移及应力分布。要求：1.编写有限元计算机程序，计算节点位移及单元应力。（划分三角形单元，单元数不得少丁30个）；2.采用有限元软件分析该问题（有限元软件网格与程序设计网格必须一致），详细给出有限元软件每一步的操作过程，并将结果与程序计算结果进行对比（任选取三个点，对比位移值）；3.提交程序编写过程的详细报告及计算机程序；4.所有同学参加答辩，并演示有限元计算程序。有限元法中三节点三角形分析结构的

2、步骤如下：1）整理原始数据，如材料性质、荷载条件、约束条件等，离散结构并进行单元编码、结点编码、结点位移编码、选取坐标系。2单元分析，建立单元刚度矩阵。整体分析，建立总刚矩阵。建立整体结构的等效节点荷载和总荷载矩阵边界条件处理。6)解方程，求出节点位移。7)求出各单元的单元应力。8)计算结果整理。一、程序设计网格划分如图，将薄板如图划分为6行，IP2)3)4)5)并建立坐标系，则34节点编单元编号刚度矩阵的集成建立与总刚度矩阵等维数的空矩阵，已变单元刚度矩阵的集成。由单元分析已知节点、单元的排布规律，继而通过循环计算求得每个单元对应的节点序号。通过循环逐个计算：（1）每个单元对应2种单元刚度矩

3、阵中的哪一种；2）该单元对应总刚度矩阵的那几行哪几列3）将该单元的单元刚度矩阵加入总刚度矩）最外层：逐行计算（2）中间层：该行逐个计算（3）最里层：区分为第奇/偶数个计算66k1565666k2565666kz5656边界约束的处理：划0置1法适用：这种方法适用丁边界节点位移分量为已知(含为0)的各种约束。做法：(1)将总刚矩阵K中相应丁已知位移行主对角线元素置1,其他元素改为零；同时将载荷列阵R中相应元素用已知位移置换。这样，由该方程求得的此位移值一定等丁已知阵的对应行列循环乂分为3层循环：（1蚌k；ke单元刚度的集成：zKkiek；kZ5量。(2)将K中已知位移相应的列的非主对角成元素也置

4、0,以保持K的对称性。当然，在已知位移分量不为零的情况下，这样做就改变了方程左端的数值，为保证方程成立，须在方程右端减去已知位移对该方程的贡献已知位移和相应总刚元素的乘积。若约束为零位移约束时，此步则可省去。特点：(1)经以上处理同样可以消除刚性位移(约束足够的前提下)，去掉未知约束反力。(2)但这种方法不改变方程阶数，利丁存贮。(3)不过，若是要求出约束反力，仍要重新计算各个划去的总刚元素。程序如下：变量说明NNODE单元节点数NPION总结点数NELEM单元数NVFIX受约束边界点数FIXED约束信息数组NFORCE节点力数FORCE节点力数组6COORD结构节点坐标数组LNODS单元定义

5、数组YOUNG弹性模量POISS泊松比THICK厚度B单元应变矩阵（3*6）D单元弹性矩阵（3*3）S单元应力矩阵（3*6）A单元面积ESTIF单元刚度矩阵ASTIF总体刚度矩阵ASLOD总体荷载向量ASDISP节点位移向量ELEDISP单元节点位移向量STRESS单元应力%*咐始化clearformatshorte%设定输出类型clear%青除内存变量NELEM=36欧元个数（单元编码总数）NPION=28%结点个数（结点编码总数）NVFIX=2%s约束边界点数NFORCE=1咐点荷载个数7YOUNG=2e11姗性模量POISS=0.25购白松比THICK=0.1钏度LNODS=123;24

6、5;253;356;478;485;589;596;6910;71112;7128;81213;8139;91314;91410;101415;111617;111712;121718;121813;131819;131914;141920;142015;152021;162223;162317;172324;172418;182425;182519;251926;192620;202627;202721;212728元定义数组（单元结点号）祯目应为单元结点号（编码）、按逆时针顺序输入COORD=00;-0.751.5;0.751.5;-1.53;03;1.53;-2.254.5;-0.754

7、.5;0.754.5;2.254.5;-36;-1.56;06;1.56;36;-3.757.5;-2.257.5;-0.757.5;0.757.5;2.257.5;3.757.5;-4.59;-39;-1.59;09;1.59;39;4.59%吉点坐标数组%坐标：x,y坐标（共NPOIN组）FORCE=10-15%结点力数组（受力结点编号，x方向,y方向）FIXED=2211;2811%勺束信息（约束点，x约束，y约束）%有约束为1,无约8束为09*%处成单元刚度矩阵并组成总体刚度矩阵ASTIF=zeros(2*NPION,2*NPION);度矩阵并置0%*fori=1:NELEMD=1PO

8、ISS0;POISS10;0(1-POISS)/2*YOUNG/(1-POISSA2)%*咐算当前单元的面积A=-det(1COORD(LNODS(i,1),1)COORD(LNODS(i,1),2);1COORD(LNODS(i,2),1)COORD(LNODS(i,2),2);1COORD(LNODS(i,3),1)COORD(LNODS(i,3),2)/2%*%成应变矩阵Bforj=0:2b(j+1)=COORD(LNODS(i,(rem(j+1),3)+1),2)-COORD(LNODS(i,(rem(j+2),3)+1),2);c(j+1)=-COORD(LNODS(i,(rem(j

9、+1),3)+1),1)+COORD(LNODS(i,(rem(j+2),3)+1),1);%生成特定大小总体刚10endB=b(1)0b(2)0b(3)0;0c(1)0c(2)0c(3);c(1)b(1)c(2)b(2)c(3)b(3)/(2*A);B1(:，:,i)=B;%*咻应力矩阵S=D*BS=D*B;ESTIF=B*S*THICK*A;%求解单元刚度矩阵a=LNODS(i,:);%|缶时向量，用来记录当前单元的节点编号forj=1:3fork=1:3ASTIF(a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)=ASTIF(a(j)*2-1):a(j)*2,(a(

10、k)*2-1):a(k)*2)+ESTIF(j*2-1:j*2,k*2-1:k*2);%根据节点编号对应关系将单元刚度分块叠加到总刚咖矩阵中endendend%*%务约束信息加入总体刚度矩阵(对角元素改一法)fori=1:NVFIXifFIXED(i,2)=111ASTIF(:,(FIXED(i,1)*2-1)=0;%H列为零ASTIF(FIXED(i,1)*2-1),:)=0;%H行为零ASTIF(FIXED(i,1)*2-1),(FIXED(i,1)*2-1)=1;%对角元素为1end%*咻成单元刚度矩阵并组成总体刚度矩阵%*ifFIXED(i,3)=1ASTIF(:,FIXED(i,1)

11、*2)=0;n列为零ASTIF(FIXED(i,1)*2,:)=0;%一行为零ASTIF(FIXED(i,1)*2,FIXED(i,1)*2)=1;冽角元素为1endend%*州成荷载向量12fori=1:NFORCEASLOD(FORCE(i,1)*2-1):FORCE(i,1)*2)=FORCE(i,2:3);end%*咻解内力ASDISP=ASTIFASLOD咐算节点位移向量ELEDISP(1:6)=0;%当前单元节点位移向量fori=1:NELEMforj=1:3ELEDISP(j*2-1:j*2)=ASDISP(LNODS(i,j)*2-1:LNODS(i,j)*2);咖出当前单元的

12、节点位移向量endiSTRESS=D*B1(:,:,i)*ELEDISP郦内力end(程序计算结果和有限元软件得出的结果稍有偏差，可能是程序某些地方数据输入时出了问题，还在寻找具体原因)二、有限元软件分析ASLOD(1:2*NPION)=0;%总体荷载向量置零1311PropertyValueUnits归ReferenceTemperatue22cd6H囹IsotropeEta&hdty7DerivefromYoungs二8YounModulus2E+11Pa二gPoissonsFtatio0.2S10EiJkModulusU533E+11Pa11ShearModulus3E+10Pa

13、12H罔FieWVariables设置材料参DetailsotCurtatebody*GraphicsPrpertisADefinitionSupprssfdstiffnessBehaviorFlexib-leCoordinateSystem&efaulttoardinateSystemReferenceTemperatureflyEnvironmentThickness100,mmThi而此ModeManualOffsetTypeMiddle-MaterialAwignmntStruttjraJStelNonlinearEffectsYesV14添加载DetskofEdgeSiznc

14、-Shirg边界旬Project0童|HoddA40-愈Geonei宅SurfaceEody刁/XCooHinabeS尊拒ms.点Gob甘CcordirtesvstentglMesh“，i 旬AHnanglesMethod-?白StaticStructuralfA5)Analysis兄ttirgmE?富Solutiofi(A6)：-?ScLibcninfitnrationUrrojrcrE囹Model(同三、励feoireTj/匚川独熊畸夺；火湖盎踏关竹了去GobaCwriinaTS/steni3-.；(5lMe?h鹿”TnandasMelioc-具魂3网白日StaticStructuralA

15、S):/上MalyseSettngs三?童Solution(A6JGlutenlrferrateScopeSeepingRetiedGeenrSelectiorGeometryiBodyDefinitiafi轴ppres盆NoMethodT沽n。否EemertMidadeNodesUseGlobjlSdtTijT-ScopeScDpirg吊曲，GacmetieleannGecmetT3Ed牌-Mnition5u叩割讲TypeNumbefofDi击砒Nuwer:fOMasiEenavixHardTBKTypeNoBias网格划分DetaisofAllFrianglesMethodtMetnoc15约束16添加载A：staticstnietuulEquSvnlEms;Type;Eqjn/dlent(wcn-rvtises)Skm-Tcp/DottciiiUnit:MP*Tim