结构力学课程设计(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 ）消去法求解如下方程组，并给出详细 的程序注解和说明：

??????????

????????=???????????????????????????????????1536353424543214109753910862781071567554321x x x x x

问题描述：

这是一个五元线性方程组，需要采用全选主元高斯约当消去法来求解。 程序说明：

全选主元高斯约当消去法：

AGJDN(A,B,N,M,L,JS)

A ——双精度实型二维数组，体积为N ×N ，输入参数。存放方程组的系数矩阵，返回时将被破坏。

B ——双精度实型二维数组，体积为N ×M ，输入兼输出参数。调用时存放M 组常数向量；返回M 组解向量。

N ——整型变量，输入参数，方程组阶数。

M ——整型变量，输入参数。方程组右端常数向量的组数。

L ——整型变量，输出参数。若返回L=0，说明方程组系数矩阵奇异，求解失败；若L ≠0，表示正常返回。

JS ——整型一维数组，长度为N 。本子程序的工作数组。

全选主元高斯约当消去法的程序及注解如下：

子程序：

SUBROUTINE AGJDN(A,B,N,M,L,JS)

DIMENSION A(N,N),B(N,M),JS(N)

DOUBLE PRECISION A,B,D

L=1

DO 100 K=1,N

Q=0.0

DO 10 I=K,N

DO 10 J=K,N

IF (ABS(A(I,J)).GT.Q) THEN

Q=ABS(A(I,J))

JS(K)=J

IS=I

END IF

10 CONTINUE

IF (Q+1.0.EQ.1.0) THEN

WRITE(*,20)

L=0

RETURN

END IF

20 FORMA T(1X,'FAIL')

DO 30 J=K,N

D=A(K,J)

A(K,J)=A(IS,J)

A(IS,J)=D

30 CONTINUE

DO 40 J=1,M

D=B(K,J)

B(K,J)=B(IS,J)

B(IS,J)=D

40 CONTINUE

DO 50 I=1,N

D=A(I,K)

A(I,K)=A(I,JS(K))

A(I,JS(K))=D

50 CONTINUE

DO 60 J=K+1,N

60 A(K,J)=A(K,J)/A(K,K)

DO 70 J=1,M

70 B(K,J)=B(K,J)/A(K,K)

DO 90 I=1,N

IF (I.NE.K) THEN

DO 80 J=K+1,N

80 A(I,J)=A(I,J)-A(I,K)*A(K,J)

DO 85 J=1,M

85 B(I,J)=B(I,J)-A(I,K)*B(K,J)

END IF

90 CONTINUE

100 CONTINUE

DO 110 K=N,1,-1

DO 110 J=1,M

D=B(K,J)

B(K,J)=B(JS(K),J)

B(JS(K),J)=D

110 CONTINUE

RETURN

END

主程序：

DIMENSION A(5,5),B(5,1),JS(5)

DOUBLE PRECISION A,B

DA TA

A/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/

DA TA B/24.0,34.0,35.0,36.0,15.0/

N=5

M=1

CALL AGJDN(A,B,N,M,L,JS)

IF (L.NE.0) THEN

WRITE(*,10) ((B(I,J),I=1,5),J=1,1)

END IF

10 FORMAT(1X,4D15.6)

END

运行结果：

问题二：

二、 利用追赶法求解如下方程组，并给出详细的程序注解和说明： ??????

????????????

????????-----=?????????????????????????????????????????????????????---------------72922206431613715211492316521131655232131165214387654321x x x x x x x x 问题描述：

这是一个五条对角线的带型的方程组，问题要求需要用追赶法来解。

方法一：追赶法

程序说明：

pendag(a,b,c,d,e,r,u,n)

n ——整型变量，输入参数，方程阶数

a ——n 个元素的一维实型数组，输入参数，存放系数矩阵的下次对角元素下面的元素（,,,321a a a ...n a ),其中21,a a 任意

b ——n 个元素的一维实型数组，输入参数，存放系数矩阵的下次对角元素（,,,321b b b ...n a ),其中1b 任意

c ——n 个元素的一维实型数组，输入参数，存放系数矩阵的对角元素

d ——n 个元素的一维实型数组，输入参数，存放系数矩阵的上次对角元素（,,,321d d d ...n d ),其中n d 任意

e ——n 个元素的一维实型数组，输入参数，存放系数矩阵的上次对角元素上面的元素（,,21e e ...,n n n e e e ,,12--)，其中n n e e ,1-任意

r ——n 个元素的一维实型数组，输入参数，存放方程的右端向量

u ——n 个元素的一维实型数组，输出参数，输出方程的解向量

追赶法带型的计算程序及注解：

SUBROUTINE pendag(a,b,c,d,e,r,u,n)

PARAMETER(nmax=100)

REAL a(n),b(n),c(n),d(n),e(n),r(n),u(n),w(nmax),beta(nmax),alpha(nmax),cg(nmax),h(nmax) INTEGER k,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

Do k=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.0 in pendag'

Beta(k+1)=(d(k)-cg(k)*alpha(k))/w(k)

alpha(k+1)=e(k)/w(k)

End do

H(1)=0.0

H(2)=r(1)/w(1)

Do k=2,n

H(k+1)=(r(k)-a(k)*h(k-1)-cg(k)*h(k))/w(k)

End do

U(n)=h(n+1)

U(n-1)=h(n)-beta(n)*u(n)

Do k=n-2,1,-1

U(k)=h(k+1)-beta(k+1)*u(k+1)-alpha(k+1)*u(k+2)

End do

End SUBROUTINE pendag

PROGRAM D1R4

! Driver program for routine PENDAG

PARAMETER(N=8)

DIMENSION A1(N,N),A(N),B(N),C(N),D(N),E(N),R(N),U(N),X(N) DA TA A1/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./

DA TA R/13.,-6.,-31.,64.,-20.,-22.,-29.,7./

Print*,'已知的方程的右端向量'

Do I=1,N

WRITE(*,'(1X,3F12.6)') R(I)

END DO

DO I=3,N

A(I)=A1(I,I-2)

END DO

DO I=2,N

B(I)=A1(I,I-1)

END DO

DO I=1,N-1

D(I)=A1(I,I+1)

END DO

DO I=1,N-2

E(I)=A1(I,I+2)

END DO

DO I=1,N

C(I)=A1(I,I)

END DO

Call PENDAG(A,B,C,D,E,R,U,N)

WRITE(*,*)

Print*,'计算出的方程的解'

DO I=1,N

WRITE(*,'(1X,3F12.6)')U(I)

END DO

!将计算的解B乘以系数矩阵，以检验计算结果的正确

DO L=1,N

X(L)=0.

DO J=1,N

X(L)=X(L)+A1(L,J)*U(J)

END DO

END DO

WRITE(*,*)

Print*,'计算出的解乘以系数矩阵的结果'

DO I=1,N

WRITE(*,'(1X,3F12.6)') X(I)

END DO

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，表示正常返回。

列选主元高斯消去法算带型计算程序及注解：

子程序：

SUBROUTINE ABAND(B,D,N,L,IL,M,IT)

DIMENSION B(N,IL),D(N,M)

DOUBLE PRECISION B,D,T

IT=1

IF (IL.NE.2*L+1) THEN

IT=-1

WRITE(*,20)

RETURN

END IF

LS=L+1

DO 100 K=1,N-1

P=0.0

DO 10 I=K,LS

IF (ABS(B(I,1)).GT.P) THEN

P=ABS(B(I,1))

IS=I

END IF

10 CONTINUE

IF (P+1.0.EQ.1.0) THEN

IT=0

WRITE(*,20)

RETURN

END IF

20 FORMA T(1X,'***FAIL***')

DO 30 J=1,M

T=D(K,J)

D(K,J)=D(IS,J)

D(IS,J)=T

30 CONTINUE

DO 40 J=1,IL

T=B(K,J)

B(K,J)=B(IS,J)

B(IS,J)=T

40 CONTINUE

DO 50 J=1,M

50 D(K,J)=D(K,J)/B(K,1)

DO 60 J=2,IL

60 B(K,J)=B(K,J)/B(K,1)

DO 90 I=K+1,LS

T=B(I,1)

DO 70 J=1,M

70 D(I,J)=D(I,J)-T*D(K,J)

DO 80 J=2,IL

80 B(I,J-1)=B(I,J)-T*B(K,J)

B(I,IL)=0.0

90 CONTINUE

IF (LS.NE.N) LS=LS+1

100 CONTINUE

IF (ABS(B(N,1))+1.0.EQ.1.0) THEN

IT=0

WRITE(*,20)

RETURN

END IF

DO 110 J=1,M

110 D(N,J)=D(N,J)/B(N,1)

JS=2

DO 150 I=N-1,1,-1

DO 120 K=1,M

DO 120 J=2,JS

120 D(I,K)=D(I,K)-B(I,J)*D(I+J-1,K) IF (JS.NE.IL) JS=JS+1

150 CONTINUE

RETURN

END

主程序：

DIMENSION B(8,5),D(8,1)

DOUBLE PRECISION B,D

DA TA

B/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/

DA TA D/13.0,-6.0,-31.0,64.0,-20.0,-22.0,-29.0,7.0/

CALL ABAND(B,D,8,2,5,1,IT)

IF (IT.GT.0) THEN

WRITE(*,10) ((D(I,J),J=1,1),I=1,8)

END IF

10 FORMAT(1X,3D15.6)

END

运行结果：

反思与对比（收获）：

对比发现二者的计算结果有点出入，咨询过老师后，才知道这是由于双精度的原因，计算机计算存在点误差，属正常现象。

问题三：

三、编写Fortran 完整程序，完成对图示刚架结构的节点力、约束力的求解，已知各杆E = 30Mpa, A = 0.18m2, I = 5.4 ×10?3m4 ：

问题描述：

这是一个梁单元问题，需要用Fortran程序来计算总刚矩阵，最终求出想要的结果。

程序框图：

程序特点：

问题类型：可用于计算结构力学的平面刚架问题

单元类型：直接利用杆梁单元

载荷类型：节点载荷及非节点载荷，其中非节点载荷包括均布荷载和垂直于杆件的集中荷载 材料性质：所有杆单元几何性质相同，且由相同的均匀材料组成

方程求解：结构刚度矩阵采用满阵存放，Gauss 消元过程采用《数值分析》中的列主元素消

去法

输入文件：按先处理法的要求，由手工生成输入数据文件 Sub--lsc 单元长及角度

Sub--elv 单元定位向量

Sub--esm 单刚矩阵

Sub--lsc 单元长及角度

Sub--elv 单元定位向量

Sub--eff 单元固端力 Sub--lsc 单元长及角度

Sub--elv 单元定位向量

Sub--esm 单刚矩阵

Sub--eff 单元固端力

开始

调用input 输入数据

tsm 形成结构刚度矩阵 jlp 形成结点荷载列阵

gauss 解方程组

mvn 计算和输出杆端力 结束

输出output 数据文件

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. subroutine lsc(Length & Sin & Cos):

输入哑元：m(单元号),nj,ne,x,y,ij

输出哑元：bl(杆件长度),si(正弦值),co(余弦值)

II. subroutine elv(Element Location Vector):

输入哑元：m,ne,nj,ij,jn

输出哑元：lv(单元定位数组)

III. subroutine esm(Element Stiffness Matrix):

输入哑元：e,a,zi,bl,si,co

输出哑元：ek(整体坐标系下的单刚矩阵)

IV. subroutine eff(Element Fixed-end Forces)

输入哑元：i,pf,nf,bl

输出哑元：fo(局部坐标系下单元固端力)

第二部分：主程序直接调用部分

I. subroutine tsm(Total Stiffness Matrix 计算总刚矩阵)

输入哑元：ne,nj,n,e,x,y,ij,a,zi,jn

输出哑元：tk

II. subroutine jlp(Joint Load Vector 计算结点荷载)

输入哑元：ne,nj,n,np,nf,x,y,ij,jn,pj,pf

输出哑元：p(结点荷载列矩阵)

III. subroutine gauss(带列主元素消去的高斯法)

输入(输出)哑元：tk,p,n ;(注意，算出位移后，直接存储到结点荷载列矩阵) IV. subroutine mvn(Member-end forces of elements 计算各单元的杆端力) 输入哑元：ne,nj,n,nf,e,x,y,ij,a,zi,jn,pf,p

3.文件管理：

源程序文件：pff.for

程序需读入的数据文件：input.txt

程序输出的数据文件：output

4.数据文件格式： 【输入文件格式】：

【输出文件格式】：

1. 第1部分：

每行数据依次为：结点号，结点x 方向位移，结点y 方向位移，结点转角位移

2. 第2部分：

每行数据依次为：单元号，xi F ，yi F ，i M ，xj F ，yj F ，j M

栏目

格式说明 实际需输入的数据 基本模型数据 第1行，每两个数之间用“，”

号隔开

单元个数，结点个数，总自由度，弹性模量，杆截面积，惯性矩，结点荷载个数，非结点荷载个数 结点位置信息 第2行，每两个

数之间用“，”

号隔开

依次输入各结点的坐标（x,y ） 单元结点信息 每输入一个单元

换行（回车），

两个数之间用

“，”号隔开

依次输入各单元的起点结点号和终点结点号 结点约束信息 每输入一个结点

换行（回车），

两个数之间用

“，”号隔开

按先处理法要求，输入各结点编号 结点荷载信息 每个结点荷载成

一行，每两个数

之间用“，”号

隔开

每行依次输入荷载作用的结点号，荷载方向代码，荷载大小（参考“主要变量”的叙述） 非结点荷载信息 每个非结点荷载

成一行，每两个

数之间用“，”

号隔开 每行依次输入荷载作用单元，荷载代码，荷载大小，荷载作用”长度”（参考“主要向量”的叙述）

源程序：

program PFF

implicit none

real tk(100,100),x(50),y(50),p(100),pj(50,3),pf(50,4) integer ij(50,2),jn(50,3)

integer ne,nj,n,np,nf

real e,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

call input(ne,nj,x,y,ij,jn,np,nf,pj,pf)

call tsm(ne,nj,n,e,x,y,ij,a,zi,jn,tk)

call jlp(ne,nj,n,np,nf,x,y,ij,jn,pj,pf,p)

call gauss(tk,p,n)

call mvn(ne,nj,n,nf,e,x,y,ij,a,zi,jn,pf,p)

end

subroutine input(ne,nj,x,y,ij,jn,np,nf,pj,pf)

dimension x(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

subroutine tsm(ne,nj,n,e,x,y,ij,a,zi,jn,tk)

dimension x(nj),y(nj),ij(ne,2),jn(nj,3),tk(n,n),ek(6,6),lv(6) do i=1,n

do j=1,n

tk(i,j)=0

enddo

enddo

do m=1,ne

call lsc(m,ne,nj,x,y,ij,bl,si,co)

call esm(e,a,zi,bl,si,co,ek)

call elv(m,ne,nj,ij,jn,lv)

do l=1,6

i=lv(l)

if (i/=0) then

do k=1,6

j=lv(k)

if (j/=0) tk(i,j)=tk(i,j)+ek(l,k)

enddo

endif

enddo

enddo

end

subroutine lsc(m,ne,nj,x,y,ij,bl,si,co) dimension x(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)

si=dy/bl

co=dx/bl

end

subroutine esm(e,a,zi,bl,si,co,ek) dimension ek(6,6)

c1=e*a/bl

c2=2.0*e*zi/bl

c3=3.0*c2/bl

c4=2.0*c3/bl

s1=c1*co*co+c4*si*si

s2=(c1-c4)*si*co

s3=c3*si

s4=c1*si*si+c4*co*co

s5=c3*co

s6=c2

ek(1,1)=s1

ek(1,2)=s2

ek(1,3)=-s3

ek(1,4)=-s1

ek(1,5)=-s2

ek(1,6)=-s3

ek(2,2)=s4

ek(2,3)=s5

ek(2,4)=-s2

ek(2,5)=-s4

ek(2,6)=s5

ek(3,3)=2*s6

ek(3,4)=s3

ek(3,5)=-s5

ek(3,6)=s6

ek(4,4)=s1

ek(4,5)=s2

ek(4,6)=s3

ek(5,5)=s4

ek(5,6)=-s5

ek(6,6)=2.0*s6

do i=1,5

do j=i+1,6

ek(j,i)=ek(i,j)

enddo

enddo

end

subroutine elv(m,ne,nj,ij,jn,lv)

dimension ij(ne,2),jn(nj,3),lv(6)

i=ij(m,1)

j=ij(m,2)

do k=1,3

lv(k)=jn(i,k)

lv(k+3)=jn(j,k)

enddo

end

subroutine jlp(ne,nj,n,np,nf,x,y,ij,jn,pj,pf,p)

dimension x(nj),y(nj),ij(ne,2),jn(nj,3),pj(np,3),pf(nf,4),p(n),fo(6),pe(6),lv(6) do i=1,n

p(i)=0.0

enddo

if (np>0) then

do i=1,np

j=int(pj(i,1))

k=int(pj(i,2))

l=jn(j,k)

if (l/=0) p(l)=pj(i,3)

enddo

endif

if(nf>0) then

do i=1,nf

m=int(pf(i,1))

call lsc(m,ne,nj,x,y,ij,bl,si,co)

call eff(i,pf,nf,bl,fo)

call elv(m,ne,nj,ij,jn,lv)

pe(1)=-fo(1)*co+fo(2)*si

pe(2)=-fo(1)*si-fo(2)*co