计算如图所示平面刚架各结点的位移和各梁的内力及支

算例计算如图-1所示平面刚架各结点的位移和各梁的内力及支座反力。已知：4422210,2.010,1.010EGPaImAm30kN/m30kN/m80kN图-1受分布力和集中力的平面刚架分析：首先建立有限元模型，即定义结点坐标，定义单元的结点号和材料特性，定义约束条件，给定结点力等。把杆件的连接点和集中力的作用点取为结点，并按1～7编号，其中5号结点就是集中力作用点，如下图所示。为了确定结点的坐标，我们要建立一个整体坐标系，其原点为结点1，水平向右为x轴正方向，竖直向上为y轴正方向。这样就能根据框架的尺寸确定7个结点的坐标。然后将这7个结点两两组合成6个单元，并按①～⑥编号，示于图-2中。该刚架有3个结点被固定，每个结点有3个自由度，因此共有9个自由度被约束。另外还有1个集中力，两个单元有分布荷载。这些信息即构成了有限元模型，可编制一个函数程序PlaneFrameModel来指定。文档收集自网络，仅用于个人学习30kN/m30kN/m80kN图-2单元划分图-1-/16有限元模型生成函数functionPlaneFrameModel%定义平面杆系的有限元模型%输入参数：%无%返回值：%无%说明：%该函数定义平面杆系的有限元模型数据：%gNode-------节点定义%gElement----单元定义%gMaterial---材料定义，包括弹性模量，梁的截面积和梁的抗弯惯性矩%gBC1--------约束条件%gNF---------集中力%gDF---------分布力globalgNodegElementgMaterialgBC1gNFgDF%节点坐标%xygNode=[0.0,0.0%节点10.0,4.0%节点23.0,0.0%节点33.0,4.0%节点44.5,4.0%节点56.0,0.0%节点66.0,4.0];%节点7%单元定义%节点1节点2材料号gElement=[1,2,1%单元12,4,1%单元23,4,1%单元34,5,1%单元45,7,1%单元56,7,1];%单元6%材料性质%弹性模量抗弯惯性矩截面积gMaterial=[2.1e11,2.0e-4,1.0e-2];%材料1文档收集自网络，仅用于个人学习%第一类约束条件%节点号自由度号约束值-2-/16gBC1=[1,1,0.01,2,0.01,3,0.03,1,0.03,2,0.03,3,0.06,1,0.06,2,0.06,3,0.0];%集中力%节点号自由度号集中力值gNF=[5,2,-80e3];%分布载荷（线性分布）%单元号节点1载荷值节点2载荷值自由度号gDF=[1-30e3022-15e3-15e32];return有了有限元模型数据，下一步的工作就是求解。其具体流程如下图所示文档收集自网络，仅用于个人学习图-3有限元程序流程图建立有限元模型（PlaneFrameModel）求解有限元模型（SolveModel）显示计算结果(DisplayResults）结束计算单元刚度矩阵（StiffnessMatrix）集成整体刚度矩阵（AssembleStiffnessMatrix）计算等效结点力（EquivalentNodeForce）处理约束条件求解方程组开始-3-/16计算分析主程序MATLAB源程序代码functionplane-beam%本程序为采用平面梁单元计算平面刚架的变形和内力%输入参数：无%输出结果：节点位移和单元节点力PlaneFrameModel;%定义有限元模型SolveModel;%求解有限元模型DisplayResults;%显示计算结果return;functionPlaneFrameModel%定义平面杆系的有限元模型%输入参数：%无%返回值：%无%说明：%该函数定义平面杆系的有限元模型数据：%gNode-------节点定义%gElement----单元定义%gMaterial---材料定义，包括弹性模量，梁的截面积和梁的抗弯惯性矩%gBC1--------约束条件%gNF---------集中力%gDF---------分布力globalgNodegElementgMaterialgBC1gNFgDF%节点坐标%xygNode=[0.0,0.0%节点10.0,4.0%节点23.0,0.0%节点33.0,4.0%节点44.5,4.0%节点56.0,0.0%节点66.0,4.0];%节点7%单元定义%节点1节点2材料号gElement=[1,2,1%单元1-4-/162,4,1%单元23,4,1%单元34,5,1%单元45,7,1%单元56,7,1];%单元6%材料性质%弹性模量抗弯惯性矩截面积gMaterial=[2.1e11,2.0e-4,1.0e-2];%材料1文档收集自网络，仅用于个人学习%第一类约束条件%节点号自由度号约束值gBC1=[1,1,0.01,2,0.01,3,0.03,1,0.03,2,0.03,3,0.06,1,0.06,2,0.06,3,0.0];%集中力%节点号自由度号集中力值gNF=[5,2,-80e3];%分布载荷（线性分布）%单元号节点1载荷值节点2载荷值自由度号gDF=[1-30e3022-15e3-15e32];returnfunctionSolveModel%求解有限元模型%输入参数：%无%返回值：%无%说明：%该函数求解有限元模型，过程如下%1.计算单元刚度矩阵，集成整体刚度矩阵%2.计算单元的等效节点力，集成整体节点力向量%3.处理约束条件，修改整体刚度矩阵和节点力向量%4.求解方程组，得到整体节点位移向量-5-/16globalgNodegElementgMaterialgBC1gNFgDFgKgDelta文档收集自网络，仅用于个人学习%step1.定义整体刚度矩阵和节点力向量[node_number,dummy]=size(gNode);gK=sparse(node_number*3,node_number*3);f=sparse(node_number*3,1);%step2.计算单元刚度矩阵，并集成到整体刚度矩阵中[element_number,dummy]=size(gElement);forie=1:1:element_numberk=StiffnessMatrix(ie,1);AssembleStiffnessMatrix(ie,k);end%step3.把集中力直接集成到整体节点力向量中[nf_number,dummy]=size(gNF);forinf=1:1:nf_numbern=gNF(inf,1);d=gNF(inf,2);f((n-1)*3+d)=gNF(inf,3);end%step4.计算分布力的等效节点力，并集成到整体节点力向量中[df_number,dummy]=size(gDF);foridf=1:1:df_numberenf=EquivalentNodeForce(gDF(idf,1),gDF(idf,2),gDF(idf,3),gDF(idf,4));文档收集自网络，仅用于个人学习i=gElement(gDF(idf,1),1);j=gElement(gDF(idf,1),2);f((i-1)*3+1:(i-1)*3+3)=f((i-1)*3+1:(i-1)*3+3)+enf(1:3);文档收集自网络，仅用于个人学习f((j-1)*3+1:(j-1)*3+3)=f((j-1)*3+1:(j-1)*3+3)+enf(4:6);文档收集自网络，仅用于个人学习end%step5.处理约束条件，修改刚度矩阵和节点力向量。采用乘大数法[bc_number,dummy]=size(gBC1);foribc=1:1:bc_numbern=gBC1(ibc,1);d=gBC1(ibc,2);m=(n-1)*3+d;f(m)=gBC1(ibc,3)*gK(m,m)*1e15;gK(m,m)=gK(m,m)*1e15;-6-/16end%step6.求解方程组，得到节点位移向量gDelta=gK\f;returnfunctionDisplayResults%显示计算结果%输入参数：%无%返回值：%无globalgNodegElementgMaterialgBC1gNFgDFgKgDelta文档收集自网络，仅用于个人学习fprintf('节点位移\n');fprintf('节点号x方向位移y方向位移转角\n');文档收集自网络，仅用于个人学习[node_number,dummy]=size(gNode);fori=1:node_numberfprintf('%6d%16.8e%16.8e%16.8e\n',...文档收集自网络，仅用于个人学习i,gDelta((i-1)*3+1),gDelta((i-1)*3+2),gDelta((i-1)*3+3));文档收集自网络，仅用于个人学习endfprintf('\n\n节点力\n');fprintf('轴力剪力弯矩\n');文档收集自网络，仅用于个人学习[element_number,dummy]=size(gElement);forie=1:element_numberenf=ElementNodeForce(ie);fprintf('单元号%6d节点号%6d%16.8e%16.8e%16.8e\n',...文档收集自网络，仅用于个人学习ie,gElement(ie,1),enf(1),enf(2),enf(3));fprintf('节点号%6d%16.8e%16.8e%16.8e\n',...文档收集自网络，仅用于个人学习gElement(ie,2),enf(4),enf(5),enf(6));endreturnfunctionk=StiffnessMatrix(ie,icoord)%计算单元刚度矩阵%输入参数:%ie-------单元号-7-/16%icoord--坐标系参数，可以是下面两个之一%1----整体坐标系%2----局部坐标系%返回值:%k----根据icoord的值，相应坐标系下的刚度矩阵globalgNodegElementgMaterialk=zeros(6,6);E=gMaterial(gElement(ie,3),1);I=gMaterial(gElement(ie,3),2);A=gMaterial(gElement(ie,3),3);xi=gNode(gElement(ie,1),1);yi=gNode(gElement(ie,1),2);xj=gNode(gElement(ie,2),1);yj=gNode(gElement(ie,2),2);L=((xj-xi)^2+(yj-yi)^2)^(1/2);k=[E*A/L00-E*A/L00文档收集自网络，仅用于个人学习012*E*I/L^36*E*I/L^20-12*E*I/L^36*E*I/L^2文档收集自网络，仅用于个人学习06*E*I/L^24*E*I/L0-6*E*I/L^22*E*I/L文档收集自网络，仅用于个人学习-E*A/L00E*A/L00文档收集自网络，仅用于个人学习0-12*E*I/L^3-6*E*I/L^2012*E*I/L^3-6*E*I/L^2文档收集自网络，仅用于个人学习06*E*I/L^22*E*I/L0-6*E*I/L^24*E*I/L];文档收集自网络，仅用于个人学习ificoord==1T=TransformMatrix(ie);k=T*k*transpose(T);endreturnfunctionAssembleStiffnessMatrix(ie,k)%把单元刚度矩阵集成到整体刚度矩阵%输入参数:%ie---单元号%k---单元刚度矩阵%返回值:%无globalgElementgKfori=1:1:2forj=1:1:2forp=1:1:3-8-/16forq=1:1:3m=(i-1)*3+p;n=(j-1)*3+q;M=(gElement(ie,i)-1)