结构力学课程设计word.docx
- 文档编号:25801150
- 上传时间:2023-06-14
- 格式:DOCX
- 页数:34
- 大小:125.62KB
结构力学课程设计word.docx
《结构力学课程设计word.docx》由会员分享,可在线阅读,更多相关《结构力学课程设计word.docx(34页珍藏版)》请在冰豆网上搜索。
结构力学课程设计word
结构力学课程设计
专业:
班级:
姓名:
学号:
指导老师:
日期:
2015年7月5日
目录
前言1
问题一:
3
问题描述:
3
程序说明:
3
全选主元高斯约当消去法:
3
全选主元高斯约当消去法的程序及注解如下:
4
运行结果:
6
问题二:
6
问题描述:
6
方法一:
追赶法7
程序说明:
7
追赶法带型的计算程序及注解:
7
运行结果:
9
总结与思考:
9
方法二:
列选主元高斯消去法算带型问题10
程序说明:
10
列选主元高斯消去法算带型计算程序及注解:
10
运行结果:
12
反思与对比(收获):
12
问题三:
13
问题描述:
13
程序框图:
14
程序特点:
14
1.主要变量:
15
2.子例行子程序哑元信息:
15
3.文件管理:
16
4.数据文件格式:
16
源程序:
17
输入数据如下(input.txt):
23
输出数据如下(output.txt):
23
程序运行后输出数据结果如下(需要手动打开output.txt文件):
24
总结与收获:
25
参考文献:
26
前言:
经过这学期的学习与积累,对结构力学这门课程有所收获,结构力学这门课程对我们学习飞行器设计与专业的学生来说,那就是手足的关系,因为我感觉任何航空、航天器都离不开结构的设计,只要有结构就牵涉到结构力学的分析与计算,因为航空器在空中飞行要遇到很多“挫折”,结构力学就是来分析这些个“挫折”下,看航空器能不能经受得了。
结构力学课程从内容上讲,主要涉及机构的几何组成分析,求解静定、超静定结构内力的虚功原理。
具体分析问题的方法包括力法、位移法等。
但对于复杂结构来讲,简单的手算的方法过于繁琐。
因此,由于课程设计偏重于利用Fortran语言编写有限元子程序来完成复杂结构的内力计算,我就恶补了好几天的与Fortran有关的知识,下面就现学现卖的计算了王老师给的三个问题,肯定有不妥之处,希望读者纠错。
问题一:
一、利用全选主元的高斯约当(Gauss-Joadan)消去法求解如下方程组,并给出详细
的程序注解和说明:
问题描述:
这是一个五元线性方程组,需要采用全选主元高斯约当消去法来求解。
程序说明:
全选主元高斯约当消去法:
AGJDN(A,B,N,M,L,JS)
A——双精度实型二维数组,体积为N×N,输入参数。
存放方程组的系数矩阵,返回时将被破坏。
B——双精度实型二维数组,体积为N×M,输入兼输出参数。
调用时存放M组常数向量;返回M组解向量。
N——整型变量,输入参数,方程组阶数。
M——整型变量,输入参数。
方程组右端常数向量的组数。
L——整型变量,输出参数。
若返回L=0,说明方程组系数矩阵奇异,求解失败;若L≠0,表示正常返回。
JS——整型一维数组,长度为N。
本子程序的工作数组。
全选主元高斯约当消去法的程序及注解如下:
子程序:
SUBROUTINEAGJDN(A,B,N,M,L,JS)
DIMENSIONA(N,N),B(N,M),JS(N)
DOUBLEPRECISIONA,B,D
L=1
DO100K=1,N
Q=0.0
DO10I=K,N
DO10J=K,N
IF(ABS(A(I,J)).GT.Q)THEN
Q=ABS(A(I,J))
JS(K)=J
IS=I
ENDIF
10CONTINUE
IF(Q+1.0.EQ.1.0)THEN
WRITE(*,20)
L=0
RETURN
ENDIF
20FORMAT(1X,'FAIL')
DO30J=K,N
D=A(K,J)
A(K,J)=A(IS,J)
A(IS,J)=D
30CONTINUE
DO40J=1,M
D=B(K,J)
B(K,J)=B(IS,J)
B(IS,J)=D
40CONTINUE
DO50I=1,N
D=A(I,K)
A(I,K)=A(I,JS(K))
A(I,JS(K))=D
50CONTINUE
DO60J=K+1,N
60A(K,J)=A(K,J)/A(K,K)
DO70J=1,M
70B(K,J)=B(K,J)/A(K,K)
DO90I=1,N
IF(I.NE.K)THEN
DO80J=K+1,N
80A(I,J)=A(I,J)-A(I,K)*A(K,J)
DO85J=1,M
85B(I,J)=B(I,J)-A(I,K)*B(K,J)
ENDIF
90CONTINUE
100CONTINUE
DO110K=N,1,-1
DO110J=1,M
D=B(K,J)
B(K,J)=B(JS(K),J)
B(JS(K),J)=D
110CONTINUE
RETURN
END
主程序:
DIMENSIONA(5,5),B(5,1),JS(5)
DOUBLEPRECISIONA,B
DATAA/5.0,7.0,6.0,5.0,1.0,7.0,10.0,8.0,7.0,2.0,6.0,8.0,10.0,9.0,3.0,5.0,7.0,9.0,10.0,4.0,1.0,2.0,3.0,4.0,5.0/
DATAB/24.0,34.0,35.0,36.0,15.0/
N=5
M=1
CALLAGJDN(A,B,N,M,L,JS)
IF(L.NE.0)THEN
WRITE(*,10)((B(I,J),I=1,5),J=1,1)
ENDIF
10FORMAT(1X,4D15.6)
END
运行结果:
问题二:
2、利用追赶法求解如下方程组,并给出详细的程序注解和说明:
问题描述:
这是一个五条对角线的带型的方程组,问题要求需要用追赶法来解。
方法一:
追赶法
程序说明:
pendag(a,b,c,d,e,r,u,n)
n——整型变量,输入参数,方程阶数
a——n个元素的一维实型数组,输入参数,存放系数矩阵的下次对角元素下面的元素(
...
),其中
任意
b——n个元素的一维实型数组,输入参数,存放系数矩阵的下次对角元素(
...
),其中
任意
c——n个元素的一维实型数组,输入参数,存放系数矩阵的对角元素
d——n个元素的一维实型数组,输入参数,存放系数矩阵的上次对角元素(
...
),其中
任意
e——n个元素的一维实型数组,输入参数,存放系数矩阵的上次对角元素上面的元素(
...,
),其中
任意
r——n个元素的一维实型数组,输入参数,存放方程的右端向量
u——n个元素的一维实型数组,输出参数,输出方程的解向量
追赶法带型的计算程序及注解:
SUBROUTINEpendag(a,b,c,d,e,r,u,n)
PARAMETER(nmax=100)
REALa(n),b(n),c(n),d(n),e(n),r(n),u(n),w(nmax),beta(nmax),alpha(nmax),cg(nmax),h(nmax)
INTEGERk,n
w
(1)=c
(1)
Beta
(1)=0.0
Beta
(2)=d
(1)/w
(1)
alpha
(1)=0.0
alpha
(2)=e
(1)/w
(1)
alpha(n)=0.0
alpha(n+1)=0.0
Dok=2,n
Cg(k)=b(k)-a(k)*beta(k-1)
W(k)=c(k)-a(k)*alpha(k-1)-cg(k)*beta(k)
If(w(k).eq.0.0)pause'w(k)=0.0inpendag'
Beta(k+1)=(d(k)-cg(k)*alpha(k))/w(k)
alpha(k+1)=e(k)/w(k)
Enddo
H
(1)=0.0
H
(2)=r
(1)/w
(1)
Dok=2,n
H(k+1)=(r(k)-a(k)*h(k-1)-cg(k)*h(k))/w(k)
Enddo
U(n)=h(n+1)
U(n-1)=h(n)-beta(n)*u(n)
Dok=n-2,1,-1
U(k)=h(k+1)-beta(k+1)*u(k+1)-alpha(k+1)*u(k+2)
Enddo
EndSUBROUTINEpendag
PROGRAMD1R4
!
DriverprogramforroutinePENDAG
PARAMETER(N=8)
DIMENSIONA1(N,N),A(N),B(N),C(N),D(N),E(N),R(N),U(N),X(N)
DATAA1/3.,-2.,1.,0.,0.,0.,0.,0.,-4.,-5.,3.,2.,&
0.,0.,0.,0.,1.,6.,-1.,5.,-3.,0.,0.,0.,0.,1.,2.,-5.,1.,6.,0.,0.,&
0.,0.,-3.,6.,-1.,1.,-4.,0.,0.,0.,0.,-1.,2.,-3.,1.,5.,0.,0.,0.,0.,&
-5.,2.,-1.,1.,0.,0.,0.,0.,0.,-9.,2.,-7./
DATAR/13.,-6.,-31.,64.,-20.,-22.,-29.,7./
Print*,'已知的方程的右端向量'
DoI=1,N
WRITE(*,'(1X,3F12.6)')R(I)
ENDDO
DOI=3,N
A(I)=A1(I,I-2)
ENDDO
DOI=2,N
B(I)=A1(I,I-1)
ENDDO
DOI=1,N-1
D(I)=A1(I,I+1)
ENDDO
DOI=1,N-2
E(I)=A1(I,I+2)
ENDDO
DOI=1,N
C(I)=A1(I,I)
ENDDO
CallPENDAG(A,B,C,D,E,R,U,N)
WRITE(*,*)
Print*,'计算出的方程的解'
DOI=1,N
WRITE(*,'(1X,3F12.6)')U(I)
ENDDO
!
将计算的解B乘以系数矩阵,以检验计算结果的正确
DOL=1,N
X(L)=0.
DOJ=1,N
X(L)=X(L)+A1(L,J)*U(J)
ENDDO
ENDDO
WRITE(*,*)
Print*,'计算出的解乘以系数矩阵的结果'
DOI=1,N
WRITE(*,'(1X,3F12.6)')X(I)
ENDDO
END
运行结果:
总结与思考:
用完了追赶法,我又偶然发现这个也可以用第一题类似的算法——列选主元高斯消去法算带型问题(方法二):
方法二:
列选主元高斯消去法算带型问题
程序说明:
ABAND(B,D,N,L,IL,M,IT)
B——双精度实型二维数据,体积为N×B,输入参数。
存放带型矩阵A中带区内的元素,该数组在返回时被破坏。
D——双精度实型二维数组,体积为N×M,输入兼输出参数。
调用时存放方程组右端的M组常数向量;返回M组解向量。
N——整型变量,输入参数。
方程组的阶数。
L——整型变量,输入参数。
带型矩阵A的半带宽。
IL——整型变量,输入参数。
存放带型矩阵A的带宽,要求IL=2*L+1.
M——整型变量,输入参数。
方程组右端常数向量的组数。
IT——整型变量,输出参数。
若返回IT<0,说明输入参数IL与L的关系不对;若IT=0,说明系数矩阵A奇异,求解失败;若IT>0,表示正常返回。
列选主元高斯消去法算带型计算程序及注解:
子程序:
SUBROUTINEABAND(B,D,N,L,IL,M,IT)
DIMENSIONB(N,IL),D(N,M)
DOUBLEPRECISIONB,D,T
IT=1
IF(IL.NE.2*L+1)THEN
IT=-1
WRITE(*,20)
RETURN
ENDIF
LS=L+1
DO100K=1,N-1
P=0.0
DO10I=K,LS
IF(ABS(B(I,1)).GT.P)THEN
P=ABS(B(I,1))
IS=I
ENDIF
10CONTINUE
IF(P+1.0.EQ.1.0)THEN
IT=0
WRITE(*,20)
RETURN
ENDIF
20FORMAT(1X,'***FAIL***')
DO30J=1,M
T=D(K,J)
D(K,J)=D(IS,J)
D(IS,J)=T
30CONTINUE
DO40J=1,IL
T=B(K,J)
B(K,J)=B(IS,J)
B(IS,J)=T
40CONTINUE
DO50J=1,M
50D(K,J)=D(K,J)/B(K,1)
DO60J=2,IL
60B(K,J)=B(K,J)/B(K,1)
DO90I=K+1,LS
T=B(I,1)
DO70J=1,M
70D(I,J)=D(I,J)-T*D(K,J)
DO80J=2,IL
80B(I,J-1)=B(I,J)-T*B(K,J)
B(I,IL)=0.0
90CONTINUE
IF(LS.NE.N)LS=LS+1
100CONTINUE
IF(ABS(B(N,1))+1.0.EQ.1.0)THEN
IT=0
WRITE(*,20)
RETURN
ENDIF
DO110J=1,M
110D(N,J)=D(N,J)/B(N,1)
JS=2
DO150I=N-1,1,-1
DO120K=1,M
DO120J=2,JS
120D(I,K)=D(I,K)-B(I,J)*D(I+J-1,K)
IF(JS.NE.IL)JS=JS+1
150CONTINUE
RETURN
END
主程序:
DIMENSIONB(8,5),D(8,1)
DOUBLEPRECISIONB,D
DATAB/3.0,-2.0,1.0,2.0,-3.0,6.0,-4.0,5.0,-4.0,-5.0,3.0,5.0,5*1.0,6.0,-1.0,-5.0,-1.0,-3.0,-1.0,-7.0,0.0,1.0,2.0,6.0,3*2.0,3*0.0,-3.0,-1.0,-5.0,-9.0,2*0.0/
DATAD/13.0,-6.0,-31.0,64.0,-20.0,-22.0,-29.0,7.0/
CALLABAND(B,D,8,2,5,1,IT)
IF(IT.GT.0)THEN
WRITE(*,10)((D(I,J),J=1,1),I=1,8)
ENDIF
10FORMAT(1X,3D15.6)
END
运行结果:
反思与对比(收获):
对比发现二者的计算结果有点出入,咨询过老师后,才知道这是由于双精度的原因,计算机计算存在点误差,属正常现象。
问题三:
三、编写Fortran完整程序,完成对图示刚架结构的节点力、约束力的求解,已知各
杆E=30Mpa,A=0.18m2,I=5.4×10−3m4:
问题描述:
这是一个梁单元问题,需要用Fortran程序来计算总刚矩阵,最终求出想要的结果。
程序框图:
程序特点:
问题类型:
可用于计算结构力学的平面刚架问题
单元类型:
直接利用杆梁单元
载荷类型:
节点载荷及非节点载荷,其中非节点载荷包括均布荷载和垂直于杆件的集中荷载
材料性质:
所有杆单元几何性质相同,且由相同的均匀材料组成
方程求解:
结构刚度矩阵采用满阵存放,Gauss消元过程采用《数值分析》中的列主元素消去法
输入文件:
按先处理法的要求,由手工生成输入数据文件
1.主要变量:
ne:
单元个数nj:
结点个数n:
自由度e:
弹性模量(单位:
KN/m2)
a:
杆截面积zi:
惯性矩np:
结点荷载个数nf:
非结点荷载个数
x(nj):
存放结点的x轴坐标y(nj):
存放结点的y轴坐标
ij(ne,2):
存放单元结点编号,其中ij(nj,1)存放起始结点编号,ij(nj,2)存放终止结点编号
jn(nj,3):
存放结点位移编号,以组成单元定位数组
pj(np,3):
存放结点荷载信息,其中pj(np,1)存放结点荷载作用结点号,pj(np,2)存放荷载方向代码(1—x方向;2—y方向;3—转角),pj(np,3)存放荷载大小
pf(ne,4):
存放非结点荷载信息,其中pf(ne,1)存放荷载作用单元号,pf(ne,2)存放荷载代码(1—均布荷载,2—垂直集中荷载),pf(ne,3)存放荷载大小,pf(ne,4)荷载作用距离(均布荷载,集中荷载均以单元起始结点为计算起始位置)。
2.子例行子程序哑元信息:
第一部分:
基本部分
I.subroutinelsc(Length&Sin&Cos):
输入哑元:
m(单元号),nj,ne,x,y,ij
输出哑元:
bl(杆件长度),si(正弦值),co(余弦值)
II.subroutineelv(ElementLocationVector):
输入哑元:
m,ne,nj,ij,jn
输出哑元:
lv(单元定位数组)
III.subroutineesm(ElementStiffnessMatrix):
输入哑元:
e,a,zi,bl,si,co
输出哑元:
ek(整体坐标系下的单刚矩阵)
IV.subroutineeff(ElementFixed-endForces)
输入哑元:
i,pf,nf,bl
输出哑元:
fo(局部坐标系下单元固端力)
第二部分:
主程序直接调用部分
I.subroutinetsm(TotalStiffnessMatrix计算总刚矩阵)
输入哑元:
ne,nj,n,e,x,y,ij,a,zi,jn
输出哑元:
tk
II.subroutinejlp(JointLoadVector计算结点荷载)
输入哑元:
ne,nj,n,np,nf,x,y,ij,jn,pj,pf
输出哑元:
p(结点荷载列矩阵)
III.subroutinegauss(带列主元素消去的高斯法)
输入(输出)哑元:
tk,p,n;(注意,算出位移后,直接存储到结点荷载列矩阵)
IV.subroutinemvn(Member-endforcesofelements计算各单元的杆端力)
输入哑元:
ne,nj,n,nf,e,x,y,ij,a,zi,jn,pf,p
3.文件管理:
源程序文件:
pff.for
程序需读入的数据文件:
input.txt
程序输出的数据文件:
output
4.数据文件格式:
【输入文件格式】:
栏目
格式说明
实际需输入的数据
基本模型数据
第1行,每两个数之间用“,”号隔开
单元个数,结点个数,总自由度,弹性模量,杆截面积,惯性矩,结点荷载个数,非结点荷载个数
结点位置信息
第2行,每两个数之间用“,”号隔开
依次输入各结点的坐标(x,y)
单元结点信息
每输入一个单元换行(回车),两个数之间用“,”号隔开
依次输入各单元的起点结点号和终点结点号
结点约束信息
每输入一个结点换行(回车),两个数之间用“,”号隔开
按先处理法要求,输入各结点编号
结点荷载信息
每个结点荷载成一行,每两个数之间用“,”号隔开
每行依次输入荷载作用的结点号,荷载方向代码,荷载大小(参考“主要变量”的叙述)
非结点荷载信息
每个非结点荷载成一行,每两个数之间用“,”号隔开
每行依次输入荷载作用单元,荷载代码,荷载大小,荷载作用”长度”(参考“主要向量”的叙述)
【输出文件格式】:
1.第1部分:
每行数据依次为:
结点号,结点x方向位移,结点y方向位移,结点转角位移
2.第2部分:
每行数据依次为:
单元号,
,
,
,
,
,
源程序:
programPFF
implicitnone
realtk(100,100),x(50),y(50),p(100),pj(50,3),pf(50,4)
integerij(50,2),jn(50,3)
integerne,nj,n,np,nf
reale,a,zi
open(1,file="input.txt",status="old")
open(2,file="output.txt",status="old")
read(1,*)ne,nj,n,e,a,zi,np,nf
callinput(ne,nj,x,y,ij,jn,np,nf,pj,pf)
calltsm(ne,nj,n,e,x,y,ij,a,zi,jn,tk)
calljlp(ne,nj,n,np,nf,x,y,ij,jn,pj,pf,p)
callgauss(tk,p,n)
callmvn(ne,nj,n,nf,e,x,y,ij,a,zi,jn,pf,p)
end
subroutineinput(ne,nj,x,y,ij,jn,np,nf,pj,pf)
dimensionx(nj),y(nj),ij(ne,2),jn(nj,3),pj(np,3),pf(nf,4)
read(1,*)(x(i),y(i),i=1,nj)
read(1,*)(ij(i,1),ij(i,2),i=1,ne)
read(1,*)((jn(i,j),j=1,3),i=1,nj)
if(np>0)read(1,*)((pj(i,j),j=1,3),i=1,np)
if(nf>0)read(1,*)((pf(i,j),j=1,4),i=1,nf)
end
subroutinetsm(ne,nj,n,e,x,y,ij,a,zi,jn,tk)
dimensionx(nj),y(nj),ij(ne,2),jn(nj,3),tk(n,n),ek(6,6),lv(6)
doi=1,n
doj=1,n
tk(i,j)=0
enddo
enddo
dom=1,ne
calllsc(m,ne,nj,x,y,ij,bl,si,co)
callesm(e,a,zi,bl,si,co,ek)
callelv(m,ne,nj,ij,jn,lv)
dol=1,6
i=lv(l)
if(i/=0)then
dok=1,6
j=lv(k)
if(j/=0)tk(i,j)=tk(i,j)+ek(l,k)
enddo
endif
enddo
enddo
end
subroutinelsc(m,ne,nj,x,y,ij,bl,si,co)
dimensionx(nj),y(nj),ij(ne,2)
i=ij(m,1)
j=ij(m,2)
dx=x(j)-x(i)
dy=y(j)-y(i)
bl=sqrt(dx*dx+dy*dy)
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 结构 力学 课程设计 word