算例计算如图1所示平面刚架各结点的位移和各梁的内力及支座反力.docx
《算例计算如图1所示平面刚架各结点的位移和各梁的内力及支座反力.docx》由会员分享，可在线阅读，更多相关《算例计算如图1所示平面刚架各结点的位移和各梁的内力及支座反力.docx（18页珍藏版）》请在冰豆网上搜索。
算例计算如图1所示平面刚架各结点的位移和各梁的内力及支座反力
算例
计算如图-1所示平面刚架各结点的位移和各梁的内力及支座反力。
已知：
图-1受分布力和集中力的平面刚架
分析：
首先建立有限元模型，即定义结点坐标，定义单元的结点号和材料特性，定义约束条件，给定结点力等。
把杆件的连接点和集中力的作用点取为结点，并按1～7编号，其中5号结点就是集中力作用点，如下图所示。
为了确定结点的坐标，我们要建立一个整体坐标系，其原点为结点1，水平向右为x轴正方向，竖直向上为y轴正方向。
这样就能根据框架的尺寸确定7个结点的坐标。
然后将这7个结点两两组合成6个单元，并按①～⑥编号，示于图-2中。
该刚架有3个结点被固定，每个结点有3个自由度，因此共有9个自由度被约束。
另外还有1个集中力，两个单元有分布荷载。
这些信息即构成了有限元模型，可编制一个函数程序PlaneFrameModel来指定。
图-2单元划分图
有限元模型生成函数
functionPlaneFrameModel
%定义平面杆系的有限元模型
%输入参数：
%无
%返回值：
%无
%说明：
%该函数定义平面杆系的有限元模型数据：
%gNode-------节点定义
%gElement----单元定义
%gMaterial---材料定义，包括弹性模量，梁的截面积和梁的抗弯惯性矩
%gBC1--------约束条件
%gNF---------集中力
%gDF---------分布力
globalgNodegElementgMaterialgBC1gNFgDF
%节点坐标
%xy
gNode=[0.0,0.0%节点1
0.0,4.0%节点2
3.0,0.0%节点3
3.0,4.0%节点4
4.5,4.0%节点5
6.0,0.0%节点6
6.0,4.0];%节点7
%单元定义
%节点1节点2材料号
gElement=[1,2,1%单元1
2,4,1%单元2
3,4,1%单元3
4,5,1%单元4
5,7,1%单元5
6,7,1];%单元6
%材料性质
%弹性模量抗弯惯性矩截面积
gMaterial=[2.1e11,2.0e-4,1.0e-2];%材料1
%第一类约束条件
%节点号自由度号约束值
gBC1=[1,1,0.0
1,2,0.0
1,3,0.0
3,1,0.0
3,2,0.0
3,3,0.0
6,1,0.0
6,2,0.0
6,3,0.0];
%集中力
%节点号自由度号集中力值
gNF=[5,2,-80e3];
%分布载荷（线性分布）
%单元号节点1载荷值节点2载荷值自由度号
gDF=[1-30e302
2-15e3-15e32];
return
有了有限元模型数据，下一步的工作就是求解。
其具体流程如下图所示
图-3有限元程序流程图
计算分析主程序
MATLAB源程序代码
functionplane-beam
%本程序为采用平面梁单元计算平面刚架的变形和内力
%输入参数：
无
%输出结果：
节点位移和单元节点力
PlaneFrameModel;%定义有限元模型
SolveModel;%求解有限元模型
DisplayResults;%显示计算结果
return;
functionPlaneFrameModel
%定义平面杆系的有限元模型
%输入参数：
%无
%返回值：
%无
%说明：
%该函数定义平面杆系的有限元模型数据：
%gNode-------节点定义
%gElement----单元定义
%gMaterial---材料定义，包括弹性模量，梁的截面积和梁的抗弯惯性矩
%gBC1--------约束条件
%gNF---------集中力
%gDF---------分布力
globalgNodegElementgMaterialgBC1gNFgDF
%节点坐标
%xy
gNode=[0.0,0.0%节点1
0.0,4.0%节点2
3.0,0.0%节点3
3.0,4.0%节点4
4.5,4.0%节点5
6.0,0.0%节点6
6.0,4.0];%节点7
%单元定义
%节点1节点2材料号
gElement=[1,2,1%单元1
2,4,1%单元2
3,4,1%单元3
4,5,1%单元4
5,7,1%单元5
6,7,1];%单元6
%材料性质
%弹性模量抗弯惯性矩截面积
gMaterial=[2.1e11,2.0e-4,1.0e-2];%材料1
%第一类约束条件
%节点号自由度号约束值
gBC1=[1,1,0.0
1,2,0.0
1,3,0.0
3,1,0.0
3,2,0.0
3,3,0.0
6,1,0.0
6,2,0.0
6,3,0.0];
%集中力
%节点号自由度号集中力值
gNF=[5,2,-80e3];
%分布载荷（线性分布）
%单元号节点1载荷值节点2载荷值自由度号
gDF=[1-30e302
2-15e3-15e32];
return
functionSolveModel
%求解有限元模型
%输入参数：
%无
%返回值：
%无
%说明：
%该函数求解有限元模型，过程如下
%1.计算单元刚度矩阵，集成整体刚度矩阵
%2.计算单元的等效节点力，集成整体节点力向量
%3.处理约束条件，修改整体刚度矩阵和节点力向量
%4.求解方程组，得到整体节点位移向量
globalgNodegElementgMaterialgBC1gNFgDFgKgDelta
%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_number
k=StiffnessMatrix（ie,1）;
AssembleStiffnessMatrix（ie,k）;
end
%step3.把集中力直接集成到整体节点力向量中
[nf_number,dummy]=size（gNF）;
forinf=1:
1:
nf_number
n=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_number
enf=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_number
n=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;
end
%step6.求解方程组，得到节点位移向量
gDelta=gK\f;
return
functionDisplayResults
%显示计算结果
%输入参数：
%无
%返回值：
%无
globalgNodegElementgMaterialgBC1gNFgDFgKgDelta
fprintf（'节点位移\n'）;
fprintf（'节点号x方向位移y方向位移转角\n'）;
[node_number,dummy]=size（gNode）;
fori=1:
node_number
fprintf（'%6d%16.8e%16.8e%16.8e\n',...
i,gDelta（（i-1）*3+1）,gDelta（（i-1）*3+2）,gDelta（（i-1）*3+3））;
end
fprintf（'\n\n节点力\n'）;
fprintf（'轴力剪力弯矩\n'）;
[element_number,dummy]=size（gElement）;
forie=1:
element_number
enf=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））;
end
return
functionk=StiffnessMatrix（ie,icoord）
%计算单元刚度矩阵
%输入参数:
%ie-------单元号
%icoord--坐标系参数，可以是下面两个之一
%1----整体坐标系
%2----局部坐标系
%返回值:
%k----根据icoord的值，相应坐标系下的刚度矩阵
globalgNodegElementgMaterial
k=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==1
T=TransformMatrix（ie）;
k=T*k*transpose（T）;
end
return
functionAssembleStiffnessMatrix（ie,k）
%把单元刚度矩阵集成到整体刚度矩阵
%输入参数:
%ie---单元号
%k---单元刚度矩阵
%返回值:
%无
globalgElementgK
fori=1:
1:
2
forj=1:
1:
2
forp=1:
1:
3
forq=1:
1:
3
m=（i-1）*3+p;
n=（j-1）*3+q;
M=（gElement（ie,i）-1）*3+p;
N=（gElement（ie,j）-1）*3+q;
gK（M,N）=gK（M,N）+k（m,n）;
end
end
end
end
return
functionenf=EquivalentNodeForce（ie,p1,p2,idof）
%计算线性分布荷载的等效节点力
%输入参数:
%ie-----单元号
%p1-----第一个节点上的分布力集度值
%p2-----第二个节点上的分布力集度值
%idof---分布力的种类，它可以是下面几种
%1---分布轴向力
%2---分布横向力
%3---分布弯矩
%返回值:
%enf-----整体坐标系下等效节点力向量
globalgElementgNode
enf=zeros（6,1）;%定义6x1的等效节点力向量
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=sqrt（（xj-xi）^2+（yj-yi）^2）;
switchidof
case1%分布轴向力
enf
（1）=（2*p1+p2）*L/6;
enf（4）=（p1+2*p2）*L/6;
case2%分布横向力
enf
（2）=（7*p1+3*p2）*L/20;
enf（3）=（3*p1+2*p2）*L^2/60;
enf（5）=（3*p1+7*p2）*L/20;
enf（6）=-（2*p1+3*p2）*L^2/60;
case3%分布弯矩
enf
（2）=-（p1+p2）/2;
enf（3）=（p1-p2）*L/12;
enf（5）=（p1+p2）/2;
enf（6）=-（p1-p2）*L/12;
otherwise
disp（sprintf（'分布力的种类错误，单元号:
%d',ie））;
end
T=TransformMatrix（ie）;%计算单元的转换矩阵
enf=T*enf;%把等效节点力转换到整体坐标下
return
functionT=TransformMatrix（ie）
%计算单元的坐标转换矩阵（局部坐标->整体坐标）
%输入参数
%ie-----节点号
%返回值
%T-------从局部坐标到整体坐标的坐标转换矩阵
globalgElementgNode
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=sqrt（（xj-xi）^2+（yj-yi）^2）;
c=（xj-xi）/L;
s=（yj-yi）/L;
T=[c-s0000
sc0000
001000
000c-s0
000sc0
000001];
return
functionenf=ElementNodeForce（ie）
%计算单元的节点力
%输入参数
%ie-----节点号
%返回值
%enf-----单元局部坐标系下的节点力
globalgElementgNodegDeltagDF
i=gElement（ie,1）;
j=gElement（ie,2）;
de=zeros（6,1）;
de（1:
3）=gDelta（（i-1）*3+1:
（i-1）*3+3）;
de（4:
6）=gDelta（（j-1）*3+1:
（j-1）*3+3）;
k=StiffnessMatrix（ie,1）;
enf=k*de;
[df_number,dummy]=size（gDF）;
foridf=1:
1:
df_number
ifie==gDF（idf,1）
enf=enf-EquivalentNodeForce（gDF（idf,1）,...
gDF（idf,2）,gDF（idf,3）,gDF（idf,4））;
break;
end
end
T=TransformMatrix（ie）;
enf=transpose（T）*enf;
return
程序计算机如果如下
节点位移
节点号x方向位移y方向位移转角
17.19450697e-019-3.10908802e-020-2.78411083e-019
27.88387267e-004-3.10908802e-005-3.44682851e-005
32.66615896e-019-1.29980350e-019-1.62999824e-019
47.70766801e-004-1.29980350e-004-2.52075453e-004
57.63456283e-004-5.67794228e-0042.15592934e-005
61.29964769e-018-7.70240078e-020-4.19430144e-019
77.56145765e-004-7.70240078e-0052.71750964e-004
节点力
轴力剪力弯矩
单元号1节点号11.63227121e+0044.76656742e+0043.56932655e+004
节点号2-1.63227121e+0041.23343258e+004-5.03056852e+003
单元号2节点号21.23343258e+0041.63227121e+0045.03056852e+003
节点号4-1.23343258e+0042.86772879e+004-2.35624322e+004
单元号3节点号36.82396838e+0042.09960018e+0036.84599262e+003
节点号4-6.82396838e+004-2.09960018e+0031.55240811e+003
单元号4节点号41.02347256e+0043.95623959e+0042.20100241e+004
节点号5-1.02347256e+004-3.95623959e+0043.73335698e+004
单元号5节点号51.02347256e+004-4.04376041e+004-3.73335698e+004
节点号7-1.02347256e+0044.04376041e+004-2.33228363e+004
单元号6节点号64.04376041e+0041.02347256e+0041.76160660e+004
节点号7-4.04376041e+004-1.02347256e+0042.33228363e+004
ANSYS计算结果对比
*****POST1NODALDEGREEOFFREEDOMLISTING*****
THEFOLLOWINGDEGREEOFFREEDOMRESULTSAREINGLOBALCOORDINATES
NODEUXUYUZROTXROTYROTZ
10.00000.00000.00000.00000.00000.0000
20.78839E-03-0.31091E-040.00000.00000.0000-0.34468E-04
30.00000.00000.00000.00000.00000.0000
40.77077E-03-0.12998E-030.00000.00000.0000-0.25208E-03
50.76346E-03-0.56779E-030.00000.00000.00000.21559E-04
60.00000.00000.00000.00000.00000.0000
70.75615E-03-0.77024E-040.00000.00000.00000.27175E-03
*****POST1NODALTOTALFORCESUMMATION*****
THEFOLLOWINGX
升级会员 