program f90_FEM_GFRMEL !*************************************************************************** !2D Frame Geometrical Nonlinear Analysis using arc-length increment method !(Water pressure problem: loads act to perpendicular to structural member) !*************************************************************************** implicit none integer::i,j,k,l,m,kk,ne,n,ia,ja character::strcom*50 !Comment integer::NODT !Number of nodes integer::NELT !Number of elements integer::MATEL !Number of material sets integer::KOX !Number of fixed nodes in x-direction integer::KOY !Number of fixed nodes in y-direction integer::KOZ !Number of fixed nodes in rotation integer::NFE !Number of loaded elements (not 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::AA(:) !Section area of element real(8),allocatable::AI(:) !Moment of inertia of element real(8),allocatable::x(:) !x-coordinate real(8),allocatable::y(:) !y-coordinate integer,allocatable::nokx(:) !Fixed node number in x-direction integer,allocatable::noky(:) !Fixed node number in y-direction integer,allocatable::nokz(:) !Fixed node number in rotation real(8)::Dpress !Water pressure increment real(8)::Tpress !Total water pressure real(8),allocatable::fel(:,:) !External force vector in local coordinate system real(8),allocatable::df(:) !Increment of external force vector real(8),allocatable::ff(:) !Accumulated external force vector integer,allocatable::mindex(:) !Index for treatment of fixed nodes real(8),allocatable::ftvec(:) !Unbalanced force vector real(8),allocatable::tsm(:,:) !Total stiffness matrix real(8),allocatable::dis(:) !Displacement increment real(8),allocatable::dist(:) !Accumulated displacement real(8),allocatable::reac(:) !work array real(8),allocatable::sresul(:,:) !Stress resultants real(8),allocatable::sm(:,:) !Element stiffness matrix real(8),allocatable::hen(:,:) !Transform matrix real(8),allocatable::fevec(:) !Element load vector integer::nt !Number of degree of freedom of all nodes (NODT x nhen) integer::nod=2 !Number of nodes per element integer::nhen=3 !Number of degree of freedom per node integer::mm !Number of degree of freedom of reduced FE equation integer::nload,iteration !Counter for loading and iteration step integer::nlmax !Limit of loading step integer::itmax !limit of iteration step real(8)::sumd1,sumd2,sumf1,sumf2 !work for judgement of convergence !variables for arc-length increment method real(8),allocatable::dis0(:) !Displacement increment due to base load increment real(8),allocatable::disR(:) !Displacement increment due to Unbalanced force real(8),allocatable::ddis(:) !Displacement increment in iterative step real(8)::deltaS0,deltaS !Arc-length real(8)::coefA,Spara !Scaling parameter real(8)::lam,dlam,dlam0 !Load control parameter real(8)::sum1,sum2,sum3,sum4 !work for calculation character::fnameR*50,fnameW*50 !Input & output file name character :: linebuf*1000,dummy*100 character :: fmt1*200,fmt2*200,fmt3*200,fmt4*200 real(4) :: t1,t2 fmt1="(i5,2(',',e15.7),3(',',e15.7),3(',',i2))" !Node fmt2="(i5,2(',',i5),3(',',e15.7),',',i5)" !Element fmt3="(i5,2(',',e15.7),3(',',e15.7),3(',',e15.7),3(',',e15.7))" !Displacement & reaction fmt4="(i5,3(',',e15.7),3(',',e15.7))" !Stress resultant call getarg(1,fnameR) !Input file name call getarg(2,fnameW) !Output file name call getarg(3,dummy) !nlmax read(dummy,*) nlmax !To convert character 'dummy' to integer 'nlmax' call getarg(4,dummy) !itmax read(dummy,*) itmax !To convert character 'dummy' to integer 'itmax' call getarg(5,dummy) !coefA read(dummy,*) coefA !To convert character 'dummy' to real number of double presision 'coefA' !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 analisys !------------------------------------------------------------------------------------------ read(11,*) NODT,NELT,MATEL,KOX,KOY,KOZ,NFE !------------------------------------------------------------------------------------------ !To declare array size !------------------------------------------------------------------------------------------ nt=NODT*nhen allocate(kakom(1:NELT,1:nod)) allocate(matno(1:NELT)) allocate(wm(1:NELT,1:3)) allocate(Em(1:NELT)) allocate(AA(1:NELT)) allocate(AI(1:NELT)) allocate(x(1:NODT)) allocate(y(1:NODT)) if(1<=KOX)allocate(nokx(1:KOX)) if(1<=KOY)allocate(noky(1:KOY)) if(1<=KOZ)allocate(nokz(1:KOZ)) allocate(fel(1:NELT,1:6)) allocate(df(1:nt)) allocate(ff(1:nt)) allocate(mindex(1:nt)) allocate(ftvec(1:nt)) allocate(tsm(1:nt,1:nt)) allocate(dis(1:nt)) allocate(dist(1:nt)) allocate(reac(1:nt)) allocate(sresul(1:NELT,1:nod*nhen)) allocate(sm(1:nod*nhen,1:nod*nhen)) allocate(hen(1:nod*nhen,1:nod*nhen)) allocate(fevec(1:nod*nhen)) allocate(dis0(1:nt)) allocate(disR(1:nt)) allocate(ddis(1:nt)) !------------------------------------------------------------------------------------------ !Material properties (Em,AA,AI) !------------------------------------------------------------------------------------------ do i=1,MATEL read(11,*) wm(i,1),wm(i,2),wm(i,3) end do !------------------------------------------------------------------------------------------ !Element-nodes relationship, material set number (node-1,node-2,matno) !------------------------------------------------------------------------------------------ do ne=1,NELT read(11,*) kakom(ne,1),kakom(ne,2),matno(ne) do i=1,MATEL if(i==matno(ne))then Em(ne) =wm(i,1) AA(ne) =wm(i,2) AI(ne) =wm(i,3) end if end do end do deallocate(wm) !------------------------------------------------------------------------------------------ ! (x,y) coordinate of node !------------------------------------------------------------------------------------------ do i=1,NODT read(11,*) x(i),y(i) end do !------------------------------------------------------------------------------------------ !Fixed node number !------------------------------------------------------------------------------------------ if(00)then mm=mm+1 mindex(mm)=i end if end do !===================================================== do nload=1,nlmax !Loading !===================================================== !----------------------------------------------------- !Initialization of matrix and vectors !----------------------------------------------------- lam=0.0D0 do i=1,nt reac(i)=0.0D0 df(i)=0.0D0 do j=1,nt tsm(i,j)=0.0D0 end do end do !----------------------------------------------------- !To assemble matrix and vector !----------------------------------------------------- do ne=1,NELT !Stiffness matrix call FRAMESM(ne,kakom,x,y,dist,Em,AA,AI,sresul,sm,NELT,NODT,nt) call ASMAT(ne,nod,nhen,kakom,sm,tsm,NELT,NODT) !To transform local loads into global loads call FHENVEC(ne,kakom,x,y,dist,Dpress,fel,fevec,NELT,NODT,nt) call ASVEC(ne,nod,nhen,kakom,fevec,df,NELT,NODT) end do !----------------------------------------------------- !Treatment of fixed nodes !----------------------------------------------------- do i=1,mm ia=mindex(i) do j=1,mm ja=mindex(j) tsm(i,j)=tsm(ia,ja) end do end do !----------------------------------------------------- !Inverse matrix !----------------------------------------------------- call MATINV(mm,tsm,nt) !dis0 do i=1,nt reac(i)=df(i) end do call CALDIS(mm,mindex,reac,tsm,nt) do i=1,nt dis0(i)=reac(i) end do !-------------------------------------------------------------------------------------- ! Treatment of 'dlam' and 'deltaS' ! + The sign of 'dlam' is defined using inner product of 'ddis' and 'dis0.' ! If inner product is positive, the angle is less than or equal to 90 degree. ! In this case, the sign of 'dlam' is defined as positive. ! + If load-displacement curve becomes horizontal, 'dlam' becomes nearly equal zero. ! It means displacement and load increments are zero. ! Therefore, when 'dlam' becomes less than or equal to 1/10 times initial value, ! arc-length 'deltaS' shall be increased. !-------------------------------------------------------------------------------------- sum1=0.0D0 sum2=0.0D0 sum4=0.0D0 do i=1,nt sum1=sum1+dis0(i)*dis0(i) sum2=sum2+df(i)*df(i) sum4=sum4+ddis(i)*dis0(i) end do if(nload==1)then deltaS0=sqrt(sum1+Spara*Spara*sum2) deltaS=deltaS0 dlam0=sqrt(deltaS*deltaS/(sum1+Spara*Spara*sum2)) end if dlam=sqrt(deltaS*deltaS/(sum1+Spara*Spara*sum2)) if(abs(dlam)<0.1D0*abs(dlam0))deltaS=deltaS*1.2D0 if(1.0D0*abs(dlam0)