您所在位置：主题：平面八节点等参单元中怎么输出高斯点的应力和应变
花败花开
[专家分：0] 发布于 2011-08-18 19:36:00
以下是主程序的一部分，单元节点应力输出正确，但是输出的高斯点的应力结果太小，不对。希望请高手帮帮忙，谢谢！
!八节点等参单元计算程序:
!******************************************************
!主程序
!******************************************************
program stress
character*64 fname1,fname2
dimension ake(16,16)
dimension ld(1000),is(16),a(30000)
dimension xjacm(2,2),xjaci(2,2)
common/sol1/nn,ne,nm,nc
common/sol2/nf,nd,n,nfd
common/shu1/mea(200,9) !单元信息
common/shu2/x(1000),y(1000) !节点坐标
common/shu3/aeu(10,3) !单元材料的弹性模量，泊松比，板厚
common/shu4/p(1000) !节点荷载
common/shu5/jz(1000,2) !节点约束信息
common/bmdmp/bmatx(3,16),dmatx(3,3)
common/sd/shape(8),deriv(2,8)
write(*,'(a22)')'input data file'
read(*,'(a)')fname1
write(*,'(a22)')'output data file'
read(*,'(a)')fname2
open(1,file=fname1)
open(2,file=fname2,status='new')
read(1,400)nn,ne,nm,nc
nf=2 !节点自由度
nd=8 !单元的节点数
nfd=nf*nd !单元自由度
n=nn*nf !单元位移的总个数
do 50 i=1,nn
50 read(1,405)k,x(i),y(i),p(i*2-1),p(i*2)
read(1,410)((jz(i,j),j=1,2),i=1,nc)
do 100 i=1,ne
read(1,415)ie,(mea(i,j),j=1,9)
100 continue
read(1,420)((aeu(i,j),j=1,3),i=1,nm)
400 format(4i5)
405 format(i5,4f15.6)
410 format(5x,i5,5x,i5)
415 format(5x,i5,5x,9i5)
420 format(3f15.6)
write(2,460)nn,ne,nm,nc
write(2,465)(i,x(i),y(i),p(i*2-1),p(i*2),i=1,nn)
write(2,470)((jz(i,j),j=1,2),i=1,nc)
write(2,475)(i,(mea(i,j),j=1,9),i=1,ne)
write(2,480)(i,(aeu(i,j),j=1,3),i=1,nm)
460 format(/5x,'总的节点数,单元数,材料数,约束数'//(10x,4i5))
465 format(/5x'节点数',5x,'横坐标',10x,'纵坐标'15x,'节点载荷'//(5x,i5,4f15.6))
470 format(/5x,'约束的节点号及其类型'//(5x,i5,5x,i5))
475 format(/'节点编号矩阵：单元号',15x,'八个节点的整体编号',12x,'材料号'//(10x,i5,5x,8i5,5x,i5))
480 format(/3x,'材料的材料号',5x,'杨氏模量',10x,'泊松比',10x,'厚度'//(5x,i5,5x,3f15.6))
call fld(nt,ld) !形成变带宽存储指示数组ld的子程序
do 500 i=1,nt !一维数组A的容量
500 a(i)=0.0 !对一维数组A赋零值
do 600 ie=1,ne !对单元循环，叠加各单元刚度得总刚度
call ke(ie,ake) !形成单元刚度矩阵的子程序
call fis(ie,is) !形成局部编号与整体编号对应关系IS数组的子程序。计算IS数组，确定单元刚度在总刚中的位置
do 560 i=1,nfd
do 560 j=1,nfd
if(is(i)-is(j))560,520,520 !判断是否属于下三角部分，是，叠加单元刚度元素；否，不做
520 ni=is(i) !总刚度元素一维存储
ij=ld(ni)-ni+is(j)
a(ij)=a(ij)+ake(i,j)
560 continue
600 continue
call fcc(nt,ld,a) !处理约束的子程序，使得方程可解
call decom(nt,a,ld) !解线性方程组的子程序
write(2,850)(i,p(2*i-1),p(2*i),i=1,nn)
call str !求应力，输出应力的子程序
call stre
850 format(//25x,'节点位移'//5x,'node no.',12x,'u',12x,'v'/(7x,i5,2f15.5))
end
!************************************************************
! 这是一个求应力输出应力的子程序
!************************************************************
subroutine str
common/sol1/nn,ne,nm,nc
common/shu3/aeu(10,3)
common/shu4/p(1000)
common/bmdmp/bmatx(3,16),dmatx(3,3) !存放应变矩阵B的数组bmatx(3,16)，存放材料矩阵D的数组dmatx(3,3)
common/shu1/mea(200,9)
dimension xjacm(2,2),xjaci(2,2)
dimension s(3,16),st(3),ve(16),is(16)
dimension ss(8),tt(8)
data ss/-1.,1.,1.,-1.,0.,1.,0.,-1./
data tt/-1.,-1.,1.,1.,-1.,0.,1.,0./
call zero(3,3,dmatx)
call zero(3,16,bmatx)
do 800 ie=1,ne !对单元循环
write(2,850)ie
nm1=mea(ie,9) !材料类别号
young=aeu(nm1,1) !弹性模量
poiss=aeu(nm1,2) !泊松比
thick=aeu(nm1,3) !板厚
call modps(young,poiss) !计算弹性矩阵D
call fis(ie,is) !计算IS数组
do 750 i=1,16 !对单元自由度循环
ni=is(i) !单元位移与整体位移的对应关系
ve(i)=p(ni) !从整体位移向量中分离出单元位移
750 continue
do 300 i=1,8
exisp=ss(i)
etasp=tt(i)
k=mea(ie,i)
call sfr2(exisp,etasp) !计算形函数及其导数
call jacob2(ie,djacb) !计算雅可比矩阵
call bmatps !计算应变矩阵B
call dot(3,3,16,dmatx,bmatx,s) !计算应力矩阵S
call dot(3,16,1,s,ve,st) !计算应力向量σ
write(2,880)i,k,st(1),st(2),st(3) !输出节点编号及节点应力
300 continue
800 continue
880 format(i5,5x,i5,5x,f15.5,5x,f15.5,5x,f15.5)
850 format(//10x,'第',i3,'个单元的节点应力'/2x,'局部编号',2x,'整体编号',5x,'U方向应力',12x,'V方向应力',12x,'剪应力')
return
end
!************************************************************
! 这是一个求高斯积分点应力的子程序
!************************************************************
subroutine stre
common/sol1/nn,ne,nm,nc
common/shu3/aeu(10,3)
common/shu4/p(1000)
common/bmdmp/bmatx(3,8),dmatx(3,3) !存放应变矩阵B的数组bmatx(3,8)，存放材料矩阵D的数组dmatx(3,3)
common/shu1/mea(200,9)
dimension xjacm(2,2),xjaci(2,2)
dimension s(3,8),st(3),ve(8),is(8)
dimension pp(4),qq(4)
data pp/-0.5773502692,0.5773502692,0.5773502692,-0.5773502692/
data qq/-0.5773502692,-0.5773502692,0.5773502692,0.5773502692/
call zero(3,3,dmatx)
call zero(3,8,bmatx)
do 990 ie=1,ne !对单元循环
write(2,995)ie
nm1=mea(ie,9) !材料类别号
young=aeu(nm1,1) !弹性模量
poiss=aeu(nm1,2) !泊松比
thick=aeu(nm1,3) !板厚
call modps(young,poiss) !计算弹性矩阵D
call fis(ie,is) !计算IS数组
do 986 i=1,8 !对单元自由度循环
ni=is(i) !单元位移与整体位移的对应关系
ve(i)=p(ni) !从整体位移向量中分离出单元位移
986 continue
do 989 i=1,4
exisp=pp(i)
etasp=qq(i)
call sfr2(exisp,etasp) !计算形函数及其导数
call jacob2(ie,djacb) !计算雅可比矩阵
call bmatps !计算应变矩阵B
call dot(3,3,8,dmatx,bmatx,s) !计算应力矩阵S
call dot(3,8,1,s,ve,st) !计算应力向量σ
write(2,994)i,st(1),st(2),st(3) !输出高斯点编号及节点应力
989 continue
990 continue
994 format(i5,5x,f15.5,5x,f15.5,5x,f15.5)
995 format(//10x,'第',i3,'个单元的高斯点应力'/2x,'U方向应力',12x,'V方向应力',12x,'剪应力')
return
end
!*************************************************************************
! 这是一个形成局部编号与整体编号对应关系的IS数组的子程序
!*************************************************************************
subroutine fis(ie,is)
common/shu1/mea(200,9)
common/sol2/nf,nd,n,nfd
dimension is(nfd)
do 100 id=1,nd !对单元节点循环
do 100 jf=1,nf !对节点自由度循环
i1=(id-1)*nf+jf
is(i1)=(mea(ie,id)-1)*nf+jf !计算IS数组的值
100 continue
return
end
!*********************************************************
!这是一个形成变带宽存储指示数组LD的子程序
!*********************************************************
subroutine fld(nt,ld)
common/shu1/mea(200,9)
common/sol1/nn,ne,nm,nc
common/sol2/nf,nd,n,nfd
dimension ld(n)
ld(1)=1 !AK的第一个元素在A中为第一个数
do 300 k=1,nn !按节点K逐一搜索Nmin并求LD
nmin=3000 !Nmin赋初值
do 200 i=1,ne !对单元循环，判断是否有K点
do 150 j=1,nd !对单元中的节点循环
if(mea(i,j).ne.k)goto 150
do 100 l=1,nd
if(mea(i,j).lt.nmin)nmin=mea(i,l) !找与K有关的最小节点号Nmin
100 continue
150 continue
200 continue
do 250 i=1,nf
j=(k-1)*nf+i !K点对应的行号
if(j.ne.1)ld(j)=ld(j-1)+(k-nmin)*nf+i !计算K点对应各行的LD
250 continue
300 continue
nt=ld(n) !计算一维数组A的容量
return
end
!*************************************************************
! 这是一个处理约束的子程序
!*************************************************************
subroutine fcc(nt,ld,a)
common/sol1/nn,ne,nm,nc
common/sol2/nf,nd,n,nfd
common/shu5/jz(1000,2)
dimension ld(n),a(nt)
do 350 k=1,nc !对每个约束循环
i=jz(k,1) !约束的节点编号
j=jz(k,2) !约束的自由度号
ni=(i-1)*nf+j !被约束位移的总体位置（行号）
nj=ld(ni) !NI行主对角线元素在A中的位置
a(nj)=1.0e25 !主对角元素置大数
350 continue
return
end
!***************************************************************
! 这是一个解线性方程组的子程序
!***************************************************************
subroutine decom(nt,a,ld)
common/shu4/p(1000)
common/sol2/nf,nd,n,nfd
dimension a(nt),ld(n)
do 20 i=1,n !这一循环完成对系数矩阵A和右端项B的分解，此为I行循环
ldi=ld(i) !第I行对角线元素在一维数组中的位置
if(i.eq.1)goto 10
io=ldi-i
im1=i-1
mi=1-io+ld(im1)
if(mi.gt.im1)goto 10
do 8 j=mi,im1
jo=ld(j)-j
ij=io+j
if(j.eq.1)goto 6
jm1=j-1
mj=1-jo+ld(jm1)
mij=max0(mi,mj)
if(mij.gt.jm1)goto 6
s=0.0
do 2 k=mij,jm1
ik=io+k
jk=jo+k
s=s+a(ik)*a(jk)
2 continue
a(ij)=a(ij)-s
6 p(i)=p(i)-a(ij)*p(j)
8 continue
10 aldi=a(ldi)
if(i.eq.1.or.mi.gt.im1)goto 13
do 12 j=mi,im1 !按列循环，完成系数矩阵的分解
ij=io+j
ldj=ld(j)
s=a(ij) !S中间变量
a(ij)=s/a(ldj)
aldi=aldi-s*a(ij)
12 continue
13 a(ldi)=aldi
p(i)=p(i)/aldi
20 continue
do 40 ldi=2,n !回代解方程
i=n-ldi+2
io=ld(i)-i
mi=1-io+ld(i-1)
im1=i-1
if(mi.gt.im1)goto 40
do 30 j=mi,im1
ij=io+j
p(j)=p(j)-a(ij)*p(i)
30 continue
40 continue
return
end
!**************************************************************
!这是一个形成单元刚度矩阵的子程序
!**************************************************************
subroutine ke(ie,ake)
dimension posgp(2) !存放高斯积分点坐标的数组
common/sol1/nn,ne,nm,nc
common/bmdmp/bmatx(3,16),dmatx(3,3)
common/shu3/aeu(10,3)
common/shu1/mea(200,9)
common/shu2/x(1000),y(1000)
dimension bt(16,3),bd(16,3)
dimension ake(16,16),ake1(16,16)
posgp(1)=-0.5773502692 !给出高斯积分点的坐标值
posgp(2)=0.5773502692
nm1=mea(ie,9)
young=aeu(nm1,1)
poiss=aeu(nm1,2)
thick=aeu(nm1,3)
call modps(young,poiss) !计算弹性矩阵D
call zero(16,16,ake)
do 80 igaus=1,2 !按Gauss积分计算单元刚度
do 80 jgaus=1,2
exisp=posgp(igaus)
etasp=posgp(jgaus)
call sfr2(exisp,etasp) !计算局部坐标系下的形状函数及其导数
call jacob2(ie,djacb) !计算Jacobi矩阵及形函数对整体坐标系的导数
call bmatps !计算应变矩阵B
call tran(3,16,bmatx,bt)
call dot(16,3,3,bt,dmatx,bd)
call dot(16,3,16,bd,bmatx,ake1)
do 70 i=1,16 !完成公式中被积函数的计算，并累加得单元刚度矩阵
do 70 j=1,16
70 ake(i,j)=ake(i,j)+djacb*thick*ake1(i,j)
80 continue
return
end
!**********************************************************
! 这是一个形成弹性矩阵的子程序
!**********************************************************
subroutine modps(y,p)
common/bmdmp/bmatx(3,16),dmatx(3,3)
call zero(3,3,dmatx)
const=y/(1.0-p*p)
dmatx(1,1)=const
dmatx(2,2)=const
dmatx(1,2)=const*p
dmatx(2,1)=const*p
dmatx(3,3)=const*(1.0-p)/2.0
return
end
!八节点等参单元计算程序:
!******************************************************
!主程序
!******************************************************
program stress
character*64 fname1,fname2
dimension ake(16,16)
dimension ld(1000),is(16),a(30000)
dimension xjacm(2,2),xjaci(2,2)
common/sol1/nn,ne,nm,nc
common/sol2/nf,nd,n,nfd
common/shu1/mea(200,9) !单元信息
common/shu2/x(1000),y(1000) !节点坐标
common/shu3/aeu(10,3) !单元材料的弹性模量，泊松比，板厚
common/shu4/p(1000) !节点荷载
common/shu5/jz(1000,2) !节点约束信息
common/bmdmp/bmatx(3,16),dmatx(3,3)
common/sd/shape(8),deriv(2,8)
write(*,'(a22)')'input data file'
read(*,'(a)')fname1
write(*,'(a22)')'output data file'
read(*,'(a)')fname2
open(1,file=fname1)
open(2,file=fname2,status='new')
read(1,400)nn,ne,nm,nc
nf=2 !节点自由度
nd=8 !单元的节点数
nfd=nf*nd !单元自由度
n=nn*nf !单元位移的总个数
do 50 i=1,nn
50 read(1,405)k,x(i),y(i),p(i*2-1),p(i*2)
read(1,410)((jz(i,j),j=1,2),i=1,nc)
do 100 i=1,ne
read(1,415)ie,(mea(i,j),j=1,9)
100 continue
read(1,420)((aeu(i,j),j=1,3),i=1,nm)
400 format(4i5)
405 format(i5,4f15.6)
410 format(5x,i5,5x,i5)
415 format(5x,i5,5x,9i5)
420 format(3f15.6)
write(2,460)nn,ne,nm,nc
write(2,465)(i,x(i),y(i),p(i*2-1),p(i*2),i=1,nn)
write(2,470)((jz(i,j),j=1,2),i=1,nc)
write(2,475)(i,(mea(i,j),j=1,9),i=1,ne)
write(2,480)(i,(aeu(i,j),j=1,3),i=1,nm)
460 format(/5x,'总的节点数,单元数,材料数,约束数'//(10x,4i5))
465 format(/5x'节点数',5x,'横坐标',10x,'纵坐标'15x,'节点载荷'//(5x,i5,4f15.6))
470 format(/5x,'约束的节点号及其类型'//(5x,i5,5x,i5))
475 format(/'节点编号矩阵：单元号',15x,'八个节点的整体编号',12x,'材料号'//(10x,i5,5x,8i5,5x,i5))
480 format(/3x,'材料的材料号',5x,'杨氏模量',10x,'泊松比',10x,'厚度'//(5x,i5,5x,3f15.6))
call fld(nt,ld) !形成变带宽存储指示数组ld的子程序
do 500 i=1,nt !一维数组A的容量
500 a(i)=0.0 !对一维数组A赋零值
do 600 ie=1,ne !对单元循环，叠加各单元刚度得总刚度
call ke(ie,ake) !形成单元刚度矩阵的子程序
call fis(ie,is) !形成局部编号与整体编号对应关系IS数组的子程序。计算IS数组，确定单元刚度在总刚中的位置
do 560 i=1,nfd
do 560 j=1,nfd
if(is(i)-is(j))560,520,520 !判断是否属于下三角部分，是，叠加单元刚度元素；否，不做
520 ni=is(i) !总刚度元素一维存储
ij=ld(ni)-ni+is(j)
a(ij)=a(ij)+ake(i,j)
560 continue
600 continue
call fcc(nt,ld,a) !处理约束的子程序，使得方程可解
call decom(nt,a,ld) !解线性方程组的子程序
write(2,850)(i,p(2*i-1),p(2*i),i=1,nn)
call str !求应力，输出应力的子程序
call stre
850 format(//25x,'节点位移'//5x,'node no.',12x,'u',12x,'v'/(7x,i5,2f15.5))
end
!************************************************************
! 这是一个求应力输出应力的子程序
!************************************************************
subroutine str
common/sol1/nn,ne,nm,nc
common/shu3/aeu(10,3)
common/shu4/p(1000)
common/bmdmp/bmatx(3,16),dmatx(3,3) !存放应变矩阵B的数组bmatx(3,16)，存放材料矩阵D的数组dmatx(3,3)
common/shu1/mea(200,9)
dimension xjacm(2,2),xjaci(2,2)
dimension s(3,16),st(3),ve(16),is(16)
dimension ss(8),tt(8)
data ss/-1.,1.,1.,-1.,0.,1.,0.,-1./
data tt/-1.,-1.,1.,1.,-1.,0.,1.,0./
call zero(3,3,dmatx)
call zero(3,16,bmatx)
do 800 ie=1,ne !对单元循环
write(2,850)ie
nm1=mea(ie,9) !材料类别号
young=aeu(nm1,1) !弹性模量
poiss=aeu(nm1,2) !泊松比
thick=aeu(nm1,3) !板厚
call modps(young,poiss) !计算弹性矩阵D
call fis(ie,is) !计算IS数组
do 750 i=1,16 !对单元自由度循环
ni=is(i) !单元位移与整体位移的对应关系
ve(i)=p(ni) !从整体位移向量中分离出单元位移
750 continue
do 300 i=1,8
exisp=ss(i)
etasp=tt(i)
k=mea(ie,i)
call sfr2(exisp,etasp) !计算形函数及其导数
call jacob2(ie,djacb) !计算雅可比矩阵
call bmatps !计算应变矩阵B
call dot(3,3,16,dmatx,bmatx,s) !计算应力矩阵S
call dot(3,16,1,s,ve,st) !计算应力向量σ
write(2,880)i,k,st(1),st(2),st(3) !输出节点编号及节点应力
300 continue
800 continue
880 format(i5,5x,i5,5x,f15.5,5x,f15.5,5x,f15.5)
850 format(//10x,'第',i3,'个单元的节点应力'/2x,'局部编号',2x,'整体编号',5x,'U方向应力',12x,'V方向应力',12x,'剪应力')
return
end
!************************************************************
! 这是一个求高斯积分点应力的子程序
!************************************************************
subroutine stre
common/sol1/nn,ne,nm,nc
common/shu3/aeu(10,3)
common/shu4/p(1000)
common/bmdmp/bmatx(3,8),dmatx(3,3) !存放应变矩阵B的数组bmatx(3,8)，存放材料矩阵D的数组dmatx(3,3)
common/shu1/mea(200,9)
dimension xjacm(2,2),xjaci(2,2)
dimension s(3,8),st(3),ve(8),is(8)
dimension pp(4),qq(4)
data pp/-0.5773502692,0.5773502692,0.5773502692,-0.5773502692/
data qq/-0.5773502692,-0.5773502692,0.5773502692,0.5773502692/
call zero(3,3,dmatx)
call zero(3,8,bmatx)
do 990 ie=1,ne !对单元循环
write(2,995)ie
nm1=mea(ie,9) !材料类别号
young=aeu(nm1,1) !弹性模量
poiss=aeu(nm1,2) !泊松比
thick=aeu(nm1,3) !板厚
call modps(young,poiss) !计算弹性矩阵D
call fis(ie,is) !计算IS数组
do 986 i=1,8 !对单元自由度循环
ni=is(i) !单元位移与整体位移的对应关系
ve(i)=p(ni) !从整体位移向量中分离出单元位移
986 continue
do 989 i=1,4
exisp=pp(i)
etasp=qq(i)
call sfr2(exisp,etasp) !计算形函数及其导数
call jacob2(ie,djacb) !计算雅可比矩阵
call bmatps !计算应变矩阵B
call dot(3,3,8,dmatx,bmatx,s) !计算应力矩阵S
call dot(3,8,1,s,ve,st) !计算应力向量σ
write(2,994)i,st(1),st(2),st(3) !输出高斯点编号及节点应力
989 continue
990 continue
994 format(i5,5x,f15.5,5x,f15.5,5x,f15.5)
995 format(//10x,'第',i3,'个单元的高斯点应力'/2x,'U方向应力',12x,'V方向应力',12x,'剪应力')
return
end
!*************************************************************************
! 这是一个形成局部编号与整体编号对应关系的IS数组的子程序
!*************************************************************************
subroutine fis(ie,is)
common/shu1/mea(200,9)
common/sol2/nf,nd,n,nfd
dimension is(nfd)
do 100 id=1,nd !对单元节点循环
do 100 jf=1,nf !对节点自由度循环
i1=(id-1)*nf+jf
is(i1)=(mea(ie,id)-1)*nf+jf !计算IS数组的值
100 continue
return
end
!*********************************************************
!这是一个形成变带宽存储指示数组LD的子程序
!*********************************************************
subroutine fld(nt,ld)
common/shu1/mea(200,9)
common/sol1/nn,ne,nm,nc
common/sol2/nf,nd,n,nfd
dimension ld(n)
ld(1)=1 !AK的第一个元素在A中为第一个数
do 300 k=1,nn !按节点K逐一搜索Nmin并求LD
nmin=3000 !Nmin赋初值
do 200 i=1,ne !对单元循环，判断是否有K点
do 150 j=1,nd !对单元中的节点循环
if(mea(i,j).ne.k)goto 150
do 100 l=1,nd
if(mea(i,j).lt.nmin)nmin=mea(i,l) !找与K有关的最小节点号Nmin
100 continue
150 continue
200 continue
do 250 i=1,nf
j=(k-1)*nf+i !K点对应的行号
if(j.ne.1)ld(j)=ld(j-1)+(k-nmin)*nf+i !计算K点对应各行的LD
250 continue
300 continue
nt=ld(n) !计算一维数组A的容量
return
end
!*************************************************************
! 这是一个处理约束的子程序
!*************************************************************
subroutine fcc(nt,ld,a)
common/sol1/nn,ne,nm,nc
common/sol2/nf,nd,n,nfd
common/shu5/jz(1000,2)
dimension ld(n),a(nt)
do 350 k=1,nc !对每个约束循环
i=jz(k,1) !约束的节点编号
j=jz(k,2) !约束的自由度号
ni=(i-1)*nf+j !被约束位移的总体位置（行号）
nj=ld(ni) !NI行主对角线元素在A中的位置
a(nj)=1.0e25 !主对角元素置大数
350 continue
return
end
!***************************************************************
! 这是一个解线性方程组的子程序
!***************************************************************
subroutine decom(nt,a,ld)
common/shu4/p(1000)
common/sol2/nf,nd,n,nfd
dimension a(nt),ld(n)
do 20 i=1,n !这一循环完成对系数矩阵A和右端项B的分解，此为I行循环
ldi=ld(i) !第I行对角线元素在一维数组中的位置
if(i.eq.1)goto 10
io=ldi-i
im1=i-1
mi=1-io+ld(im1)
if(mi.gt.im1)goto 10
do 8 j=mi,im1
jo=ld(j)-j
ij=io+j
if(j.eq.1)goto 6
jm1=j-1
mj=1-jo+ld(jm1)
mij=max0(mi,mj)
if(mij.gt.jm1)goto 6
s=0.0
do 2 k=mij,jm1
ik=io+k
jk=jo+k
s=s+a(ik)*a(jk)
2 continue
a(ij)=a(ij)-s
6 p(i)=p(i)-a(ij)*p(j)
8 continue
10 aldi=a(ldi)
if(i.eq.1.or.mi.gt.im1)goto 13
do 12 j=mi,im1 !按列循环，完成系数矩阵的分解
ij=io+j
ldj=ld(j)
s=a(ij) !S中间变量
a(ij)=s/a(ldj)
aldi=aldi-s*a(ij)
12 continue
13 a(ldi)=aldi
p(i)=p(i)/aldi
20 continue
do 40 ldi=2,n !回代解方程
i=n-ldi+2
io=ld(i)-i
mi=1-io+ld(i-1)
im1=i-1
if(mi.gt.im1)goto 40
do 30 j=mi,im1
ij=io+j
p(j)=p(j)-a(ij)*p(i)
30 continue
40 continue
return
end
!**************************************************************
!这是一个形成单元刚度矩阵的子程序
!**************************************************************
subroutine ke(ie,ake)
dimension posgp(2) !存放高斯积分点坐标的数组
common/sol1/nn,ne,nm,nc
common/bmdmp/bmatx(3,16),dmatx(3,3)
common/shu3/aeu(10,3)
common/shu1/mea(200,9)
common/shu2/x(1000),y(1000)
dimension bt(16,3),bd(16,3)
dimension ake(16,16),ake1(16,16)
posgp(1)=-0.5773502692 !给出高斯积分点的坐标值
posgp(2)=0.5773502692
nm1=mea(ie,9)
young=aeu(nm1,1)
poiss=aeu(nm1,2)
thick=aeu(nm1,3)
call modps(young,poiss) !计算弹性矩阵D
call zero(16,16,ake)
do 80 igaus=1,2 !按Gauss积分计算单元刚度
do 80 jgaus=1,2
exisp=posgp(igaus)
etasp=posgp(jgaus)
call sfr2(exisp,etasp) !计算局部坐标系下的形状函数及其导数
call jacob2(ie,djacb) !计算Jacobi矩阵及形函数对整体坐标系的导数
call bmatps !计算应变矩阵B
call tran(3,16,bmatx,bt)
call dot(16,3,3,bt,dmatx,bd)
call dot(16,3,16,bd,bmatx,ake1)
do 70 i=1,16 !完成公式中被积函数的计算，并累加得单元刚度矩阵
do 70 j=1,16
70 ake(i,j)=ake(i,j)+djacb*thick*ake1(i,j)
80 continue
return
end
!**********************************************************
! 这是一个形成弹性矩阵的子程序
!**********************************************************
subroutine modps(y,p)
common/bmdmp/bmatx(3,16),dmatx(3,3)
call zero(3,3,dmatx)
const=y/(1.0-p*p)
dmatx(1,1)=const
dmatx(2,2)=const
dmatx(1,2)=const*p
dmatx(2,1)=const*p
dmatx(3,3)=const*(1.0-p)/2.0
return
end