module defpi implicit none real(8),parameter::pi=3.14159265358979323846D0 end module defpi program f90_FEM_PLNT !***************************************************************** !2D Stress Analysis with the function of no-tension analysis !***************************************************************** use defpi implicit none integer::i,j,k,l,kk,ne,n,ia,ja character::strcom*50 !Comment integer::NSTRES !Stress state (0: plane strain, 1: plane stress) integer::NODT !Number of nodes integer::NELT !Number of elements integer::MATEL !Number of material sets integer::KOX !Number of restricted nodes in x-direction integer::KOY !Number of restricted nodes in y-direction integer::NF !Number of loaded nodes integer,allocatable::kakom(:,:) !Element connectivity integer,allocatable::matno(:) !Material set number real(8),allocatable::wm(:,:) !work for material properties real(8),allocatable::Em(:) !Elastic modulus of element real(8),allocatable::po(:) !Poisson's ratio of element real(8),allocatable::t(:) !Thickness of element real(8),allocatable::gamma(:) !Unit weight of element real(8),allocatable::gkh(:) !Horizontal acceleration of element in the ratio to 'g' real(8),allocatable::gkv(:) !Vertical acceleration of element in the ratio of 'g' real(8),allocatable::alpha(:) !Coefficient of thermal expantion of element real(8),allocatable::ts(:) !Tensile strength of element real(8),allocatable::x(:) !x-coordinate of node real(8),allocatable::y(:) !y-coordinate of node real(8),allocatable::deltaT(:) !Temperature change of node (+: temperature rize) integer,allocatable::nokx(:) !Restricted node number in x-direction integer,allocatable::noky(:) !Restricted node number in y-direction real(8),allocatable::f(:) !External load vector integer,allocatable::mindex(:) !Index for treatment of restricted nodes real(8),allocatable::ftvec(:) !Total load vector real(8),allocatable::fmass(:) !Body force vector real(8),allocatable::rdis(:) !Displacement of restricted nodes real(8),allocatable::tsm(:,:) !Stiffness matrix real(8),allocatable::dis(:) !Displacement increment of nodes real(8),allocatable::wdis(:) !work for displacement real(8),allocatable::dist(:) !Cumulated displacement of nodes real(8),allocatable::reac(:) !Internal force vector & work real(8),allocatable::strsave(:,:,:) !Memory for stresses integer,allocatable::noten(:,:) !Memory for stress state (no-tension or not) real(8),allocatable::sm(:,:) !Element stiffness matrix real(8),allocatable::fevec(:) !Element load vector integer::nt !Number of degrees of freedom of all nodes (NODT x nhen) integer::nod !Number of nodes per element integer::nhen=2 !Number of degree of freedom per node integer::mm !Number of degrees of freedom of reduced FE equation integer::ib !Band width of reduced matrix integer::kkmax !Number of integral points real(8)::elarea !Area of element real(8)::fact0=1.0D-30 !0 judgement real(8):: sigx,sigy,tauxy,ps1,ps2,ang !Stresses real(8)::xg,yg !Typical coordinate for stresses integer:: nnn !Increment for iterative calculation integer:: maxita=2000 !Limit value of iterative calculation integer:: ntsum !Number of no-tension state Gauss points integer:: icount !Number of converged degrees of freedom real(8):: dtest,ftest !Sum total of increment of displacement, sum total of unbalanced force integer:: IPR !Flag for print of stresses (0: all Gauss points, 1: element mean stresses) character::fnameR*50,fnameW*50 !Input file name & output file name character :: linebuf*1000 character :: fmt1*200,fmt2_3*200,fmt2_4*200,fmt3*200,fmt4*200 real(4) :: t1,t2 fmt1="(i5,2(',',e15.7),2(',',e15.7),2(',',i2),2(',',e15.7),',',e15.7)" !Node fmt2_3="(i5,3(',',i5),8(',',e15.7),',',i5)" !Element(3 nodes) fmt2_4="(i5,4(',',i5),8(',',e15.7),',',i5)" !Element(4 nodes) fmt3="(i5,2(',',e15.7),2(',',e15.7),2(',',e15.7),2(',',e15.7))" !Displacement & reaction fmt4="(i5,',',i5,2(',',e15.7),3(',',e15.7),3(',',e15.7),2(',',i5))" !Stresses call getarg(1,fnameR) !Input file name call getarg(2,fnameW) !Output file name !To start the measurement of calculation time t1=DIFT() !*************************** !data input !*************************** open(11,file=fnameR,status='old') !!-------------------------------------------------------------------------------------- !!Comment !!-------------------------------------------------------------------------------------- read(11,'(a)') strcom !------------------------------------------------------------------------------------------ !Basic values for analysis !------------------------------------------------------------------------------------------ read(11,*) nod,NODT,NELT,MATEL,KOX,KOY,NF,NSTRES,IPR if(nod==3)IPR=0 !------------------------------------------------------------------------------------------ !To declare array size !------------------------------------------------------------------------------------------ nt=NODT*nhen if(nod==3)kkmax=1 if(nod==4)kkmax=4 allocate(kakom(1:NELT,1:nod)) allocate(matno(1:NELT)) allocate(wm(1:NELT,1:8)) allocate(Em(1:NELT)) allocate(po(1:NELT)) allocate(t(1:NELT)) allocate(gamma(1:NELT)) allocate(gkh(1:NELT)) allocate(gkv(1:NELT)) allocate(alpha(1:NELT)) allocate(ts(1:NELT)) allocate(x(1:NODT)) allocate(y(1:NODT)) allocate(deltaT(1:NODT)) allocate(nokx(1:KOX)) allocate(noky(1:KOY)) allocate(f(1:nt)) allocate(mindex(1:nt)) allocate(ftvec(1:nt)) allocate(fmass(1:nt)) allocate(rdis(1:nt)) allocate(tsm(1:nt,1:nt)) allocate(dis(1:nt)) allocate(wdis(1:nt)) allocate(dist(1:nt)) allocate(reac(1:nt)) allocate(strsave(1:NELT,1:kkmax,1:3)) allocate(noten(1:NELT,1:kkmax)) allocate(sm(1:nod*nhen,1:nod*nhen)) allocate(fevec(1:nod*nhen)) !------------------------------------------------------------------------------------------ !Material properties (t,Em,po,gamma,gkh,gkv,alpha,ts) !------------------------------------------------------------------------------------------ do i=1,MATEL read(11,*) (wm(i,j),j=1,8) end do !------------------------------------------------------------------------------------------ !Element-nodes relationship, material set number !------------------------------------------------------------------------------------------ do ne=1,NELT read(11,*) (kakom(ne,j),j=1,nod),matno(ne) do i=1,MATEL if(i==matno(ne))then t(ne) =wm(i,1) Em(ne) =wm(i,2) po(ne) =wm(i,3) gamma(ne)=wm(i,4) gkh(ne) =wm(i,5) gkv(ne) =wm(i,6) alpha(ne)=wm(i,7) ts(ne) =wm(i,8) if(NSTRES==0)t(ne)=1.0D0 end if end do end do deallocate(wm) !------------------------------------------------------------------------------------------ !(x,y) coordinates, temperature change !------------------------------------------------------------------------------------------ do i=1,NODT read(11,*) x(i),y(i),deltaT(i) end do !------------------------------------------------------------------------------------------ !Restricted node number, displacement !------------------------------------------------------------------------------------------ do i=1,nt rdis(i)=0.0D0 end do if(00)then mm=mm+1 mindex(mm)=i end if end do do i=1,mm ia=mindex(i) ftvec(i)=ftvec(ia) do j=1,mm ja=mindex(j) tsm(i,j)=tsm(ia,ja) end do end do !----------------------------------------------------- !Calculation of band width & memory of band matrix !----------------------------------------------------- ib=IBAND(mm,tsm,fact0,nt) call BAND(mm,ib,tsm,nt) !----------------------------------------------------- !To verify diagonal element & to stop if zero value !----------------------------------------------------- do i=1,mm if(abs(tsm(i,1))0.0D0)ang=45.0D0 if(tauxy<0.0D0)ang=135.0D0 if(tauxy==0.0D0)ang=0.0D0 else ang=0.5D0*atan(2.0D0*tauxy/(sigx-sigy)) ang=180.0D0/pi*ang if(sigx>sigy.and.tauxy>=0.0D0)ang=ang if(sigx>sigy.and.tauxy<0.0D0)ang=ang+180.0D0 if(sigxm1)then mm=m1 else mm=n1-i+1 end if do j=1,mm z=array(i,j) if(j/=m1)then if(m1-j+11)then array(i,j)=z/array(i,1) else array(i,1)=sqrt(z) end if else array(i,j)=z/array(i,1) end if end do end do end subroutine CHBAND10 subroutine CHBAND11(n1,m1,array,vector,nt) integer,intent(in)::n1,m1,nt real(8),intent(in)::array(1:nt,1:nt) real(8),intent(out)::vector(nt) integer :: i,j,mm real(8) :: z vector(1)=vector(1)/array(1,1) do i=2,n1 z=vector(i) if(i>m1)then mm=m1 else mm=i end if do j=2,mm z=z-array(i-j+1,j)*vector(i-j+1) end do vector(i)=z/array(i,1) end do vector(n1)=vector(n1)/array(n1,1) do i=n1-1,1,-1 z=vector(i) do j=2,m1 if(i+j-1>n1)exit z=z-array(i,j)*vector(i+j-1) end do vector(i)=z/array(i,1) end do end subroutine CHBAND11 integer function IBAND(mm,array,fact0,nt) integer,intent(in)::mm,nt real(8),intent(in)::array(1:nt,1:nt),fact0 integer :: i,j !----------------------------------------------------- !Calculation of band width !----------------------------------------------------- IBAND=1 do i=1,mm-1 do j=mm,i+1,-1 if(fact0<=abs(array(i,j)))exit end do if(IBAND<=j-i+1)IBAND=j-i+1 end do if(mm