!************************************************* !************************************************* !************************************************* ! ELASTO-PLASTIC MINDLIN PLATE BENDING ANALYSIS ! Overall structure of MINDLAY(Continued). ! Program FEAM(Input,Output,Tape) !************************************************* ! 明德林板分层弹塑性材料力学弯曲分析 !******************************************************************** !******************************************************************** ! SINOMACH ! ! FORTRAN95,2013&2018 version ! ! !PotsyYZ,Mike,周勇，Tom, Speagaul,Fiddler,Goldberg, Johnny,SwanseaUK ! 20230427 !********************************************************************* Program main !c******************************************************************** !c !c*** ELASTO-PLASTIC ANALYSIS OF LAYERED MINDLIN PLATES USING !C*** 4-,8- , 9-NODED OR HETEROSIS ISOPARAMETRIC QUADRILATERALS !C !C******************************************************************* Dimension asdis(240), coord(80, 2), effst(225), eload(25, 27),& epstn(225), estif(27, 27), & eqrhs(10), equat(40, 10), fixed(240),& iffix(240), gload(40), gstif(860), lnods(25, 9), locel(27), & matno(25), nacva(40), namev(10), ncdis(4), ncres(4),& ndest(27), ndfro(25), nofix(40), noutp(2), npivo(10),& posgp(4), presc(40, 3), props(10, 8), refor(240),& rload(25, 27), strsg(5, 225), tofor(240),& tdisp(240), tload(25, 27), treac(40, 3), vecrv(40),& weigp(4) !c******************************************************** Open (5, File='data\input.dat', Status='unknown') Open (6, File='data\output.dat', Status='unknown') !c******************************************************** !C !C*** PRESET VARIABLES ASSOCIATED WITH DYNAMIC DIMENSIONS !C Call dimmp(mbufa, melem, mevab, mfron, mmats, mpoin,& mstif, mtotg, mtotv, mvfix, ndime, ndofn,& nprop, nstre) !C !C*** CALL THE SUBROUTINE WHICH READS MOST OF THE PROBLEM DATA !c Call input(coord, iffix, lnods, matno, melem, mevab, & mfron, mmats, mpoin, mtotv, mvfix, nalgo, & ncrit, ndfro, ndime, ndofn, nelem, nevab, & ngaus, nlaps, nincs, nmats, nnode, nofix, & npoin, nprop, nstre, nstrl, nswit, ntotg, & ntotv, ntype, nvfix, posgp, presc, props, & weigp) !c !C**** INITIALIZE ARRAYS TO ZERO !C Call zeromp(effst, eload, epstn, melem, mevab, mtotg, & mtotv, mvfix, ndofn, nelem, nevab, ngaus, & ntotg, ntotv, nvfix, strsg, tdisp, tfact, & tload, treac) !C !C**** !C Call mindpb(ifdis, iffix, ifres, lnods, melem, mtotv, & ncdis, ncres, nelem, ntype) !C !C !C !C*** COMPUTE LOAD AFTER READING RELEVANT EXTRA DATA !C Call loadpb(coord, lnods, matno, melem, mmats, mpoin, & nelem, nevab, ngaus, nnode, npoin, props, & rload) !C !C*** LOOP OVER EACH INCREMENT !C !c ELASTO-PLASTIC MINDLIN PLATE BENDINGING ANALYSIS !c Do iincs = 1, nincs !c !c*** READ DATA FOR CURRENT INCREMENT !c Call increm(eload, fixed, iincs, melem, mevab, miter, & mtotv, mvfix, ndofn, nelem, nevab, noutp, & nofix, ntotv, nvfix, presc, rload, tfact, & tload, toler) !C !C*** LOOP OVER EACH ITERATION !C Do iiter = 1, miter !c !C*** CALL ROUTINE WHICH SELECTS SOLUTION ALGORITHM VARIABLE KRESL !c Call algor(fixed, iincs, iiter, kresl, mtotv, nalgo, & ntotv) !C !C*** CHECK WHETHER A NEW EVALUATION OF THE STIFFNESS MATRICES IS NEEDED. !c If (kresl==1) Call stimpa(coord, epstn, iincs, lnods, matno, melem, & mevab, mmats, mpoin, mtotg, ncrit, nelem, & nevab, ngaus, nnode, nlaps, props, strsg) !C !C*** SOLVE EQUATIONS !C Call front(asdis, eload, eqrhs, equat, estif, fixed, & iffix, iincs, iiter, gload, gstif, kresl, & lnods, locel, mbufa, melem, mevab, mfron, & mstif, mtotv, mvfix, nacva, namev, ndest, & ndofn, nelem, nevab, nnode, nofix, npivo, & npoin, ntotv, tdisp, tload, treac, vecrv) !c !C*** CALCULATE RESIDUAL FORCES !C Call resmp(asdis, coord, effst, eload, epstn, lnods, & matno, melem, mmats, mpoin, mtotg, mtotv, & ncrit, nelem, nevab, ngaus, nnode, nlaps, & props, strsg) !C !C*** CHECK FOR CONVERGENCE !C Call convmp(asdis, eload, iiter, ifdis, ifres, lnods, & melem, mevab, mtotv, nchek, ncdis, ncres, & ndofn, nelem, nevab, nnode, npoin, ntotv, & refor, tofor, tdisp, tload, toler) !C !C*** OUTPUT RESULTS IF REQUIRED !C If (iiter==1 .And. noutp(1)>0) Call outmpa(epstn, iiter, mtotg, ktotv, mvfix, nelem, & ngaus, nlaps, nofix, noutp, npoin, nvfix, & strsg, tdisp, treac) !c !c*** IF SOLUTION HAS CONVERGED STOP ITERATING AND OUTPUT RESULTS !C If (nchek==0) Goto 100 End Do !c !c**** !c If (nalog==2) Goto 100 Stop 100 Call outmpa(epstn, iiter, mtotg, mtotv, mvfix, nelem, & ngaus, nlaps, nofix, noutp, npoin, nvfix, & strsg, tdisp, treac) End Do !c******************************************************** Close (5) Close (6) !c******************************************************** Stop End Program main !********************************************************* !********************************************************* !*********************************************************** !*********************************************************** ! DIMMP !*********************************************************** Subroutine dimmp(mbufa, melem, mevab, mfron, mmats, mpoin,& mstif, mtotg, mtotv, mvfix, ndime, ndofn,& nprop, nstre) !c************************************************************ !c !c****Sets up dynamic dimensions-must agree with dimensions !c****in femp !c !c************************************************************ mbufa=10 melem=25 mfron=40 mmats=10 mpoin=80 mstif=(mfron*mfron-mfron)/2.0+mfron mtotg=melem*9 ndofn=3 mototv=mpoin*ndofn mvfix=40 ndime=2 nprop=8 nstre=5 mevab=ndofn*9 Return End Subroutine dimmp !*************************************************************** !*************************************************************** ! input !*************************************************************** Subroutine input (coord,iffix,lnods,matno,melem,mevab,& mfron,mmats,mpoin,mtotv,mvfix,nalgo,& ncrit,ndfro,ndime,ndofn,nelem,nevab,& ngaus,nlaps,nincs,nmats,nnode,nofix,& npoin,nprop,nstre,nstrl,nswit,ntotg,& ntotv,ntype,nvfix,posgp,presc,props,& weigp) !c********************************************************************* !c !c****This subroutine accepts most of the input data !c !c********************************************************************* Dimension coord(mpoin,2), iffix(mtotv), lnods(melem,9), & matno(melem), ndfro(melem), & nofix(mvfix), posgp(4), presc(mvfix,ndofn), & props(mmats,nprop), title(12), weigp(4) ! Rewind 1 ! Rewind 2 ! Rewind 3 ! Rewind 4 ! Rewind 8 Read(5,920) title Write(6,920) title !c !c***Read the first data card, and echo it immediately. !c Read(5,900) npoin, nelem, nvfix, ntype, nnode, nmats, ngaus,& nalgo, ncrit, nincs, nstre nevab=ndofn*nnode nstr1=nstre+1 If(ntype==3) nstr1=nstre ntotv=npoin*ndofn ngau2=ngaus*ngaus ntotg=nelem*ngau2 Write(6,901) npoin, nelem, nvfix, ntype, nnode, nmats, ngaus, nevab,& nalgo, ncrit, nincs, nstre Call check1(ndofn,nelem,ngaus,nmats,nnode,npoin,& nstre,ntype,nvfix,ncrit,nalgo,nincs) !c !c***Read the element nodal connections, and the property numbers. !c Write(6,902) Do ielem=1, nelem Read(5,900) numel, matno(numel), (lnods(numel,inode),inode=1,nnode) Write(6,903) numel, matno(numel), (lnods(numel,inode),inode=1,nnode) End Do !c !c***Zero all the nodal coordinates, prior to reading some of them. !c Do ipoin=1, npoin Do idime=1, 2 coord(ipoin,idime)=0.0 End Do End Do !c !c***Read some nodal coordinates, finishing with the node of all. !c Write(6,904) 6 Read(5,905) ipoin, (coord(ipoin,idime),idime=1,2) If(ipoin/=npoin) Goto 6 !c !c***Interpolate coordinates of mid-side nodes !c Call nodexy(coord,lnods,melem,mpoin,nelem,nnode) Do ipoin=1, npoin Write(6,906) ipoin, (coord(ipoin,idime),idime=1,2) End Do !c !c***Read the fixed values. !c Write(6,907) Do ivfix=1, nvfix Read(5,908) nofix(ivfix), ifpre, (presc(ivfix,idofn),idofn=1,ndofn) Write(6,908) nofix(ivfix), ifpre, (presc(ivfix,idofn),idofn=1,ndofn) nloca=(nofix(ivfix)-1)*ndofn ifdof=10**(ndofn-1) Do idofn=1, ndofn ngash=nloca+idofn If(ifpre1) Goto 455 Do ipoin = 1, npoin klast = 0 Do ielem = 1, nelem Do inode = 1, nnode If (lnods(ielem,inode)/=ipoin) Goto 120 klast = ielem nlast = inode 120 End Do ! !c FINITE ELEMENTS IN PLASTICITY ! End Do If (klast/=0) lnods(klast, nlast) = -ipoin End Do 455 Continue !C !c*** START BY INITIALIZING EVERYTHING THAT MATTERS TO ZERO !c Do ibufa = 1, mbufa eqrhs(ibufa) = 0.0 End Do Do istif = 1, mstif gstif(istif) = 0.0 End Do Do ifron = 1, mfron gload(ifron) = 0.0 vecrv(ifron) = 0.0 nacva(ifron) = 0 Do ibufa = 1, mbufa equat(ifron, ibufa) = 0.0 End Do End Do !C !C*** AND PREPARE FOR DISC READING AND WRITING OPERATIONS !c nbufa = 0 If (kresl>1) nbufa = mbufa !? REWIND 1 !? REWIND 2 !? REWIND 3 !? REWIND 4 !? REWIND 8 !C !C*** ENTER MAIN ELEMENT ASSEMBLY-REDUCTION LOOP !C nfron = 0 kelva = 0 Do ielem = 1, nelem If (kresl>1) Goto 400 kevab = 0 Read (1) estif Do inode = 1, nnode Do idofn = 1, ndofn nposi = (inode-1)*ndofn + idofn locno = lnods(ielem, inode) If (locno>0) locel(nposi) = (locno-1)*ndofn + idofn If (locnonfron) nfron = ndest(kevab) 210 End Do Write (8) locel, ndest, nacva, nfron 400 If (kresl>1) Read (8) locel, ndest, nacva, nfron !c !C**** ASSEMBLE ELEMENT LOADS !c Do ievab = 1, nevab idest = ndest(ievab) gload(idest) = gload(idest) + eload(ielem, ievab) !C !C*** ASSEMBLE THE ELEMENT STIFFNESSES-BUT NOT IN RESOLUTION !C If (kresl>1) Goto 402 Do jevab = 1, ievab jdest = ndest(jevab) ngash = nfunc(idest, jdest) ngish = nfunc(jdest, idest) If (jdest>=idest) gstif(ngash) = gstif(ngash) + estif(ievab, jevab) If (jdest0.0) Goto 235 Write (6, 900) nikno, pivot Stop 235 Continue equat(ifron, nbufa) = 0.0 !c !C*** ENQUIRE WHETHER PRESENT VARIABLE IS FREE OR PRESCRIBED !C If (iffix(nikno)==0) Goto 250 !c !C*** DEAL WITH A PRESCRIBED DEFLECTION !C Do jfron = 1, nfron gload(jfron) = gload(jfron) - fixed(nikno)*equat(jfron, nbufa) End Do Goto 280 !c !c*** ELIMINATE A FREE VARIABLE - DEAL WITH THE RIGHT HAND SIDE FIRST !C 250 Do jfron = 1, nfron gload(jfron) = gload(jfron) - equat(jfron, nbufa)*eqrhs(nbufa)/pivot !c !c**** NOW DEAL WITH THE COEFFICIENTS IN CORE !C If (kresl>1) Goto 418 If (equat(jfron,nbufa)==0.0) Goto 270 nloca = nfunc(0, jfron) cureq = equat(jfron, nbufa) Do lfron = 1, jfron ngash = lfron + nloca gstif(ngash) = gstif(ngash) - cureq*equat(lfron, nbufa)& /pivot End Do 418 Continue 270 End Do 280 equat(ifron, nbufa) = pivot !C !C*** RECORD THE NEW VACANT SPACE, AND REDUCE FRONTWIDTH IF POSSIBLE !C nacva(ifron) = 0 Goto 290 !c !C*** COMPLETE THE ELEMENT LOOP IN THE FORWARD ELIMINATION !c 300 End Do 290 If (nacva(nfron)/=0) Goto 310 nfron = nfron - 1 If (nfron>0) Goto 290 310 End Do End Do If (kresl==1) Write (2) equat, eqrhs, npivo, namev ! !c?? BACKSPACE 2 ??? !c !c*** ENTER BACKSUBSTITUTION PHASE. LOOP BACKWARDS THROUGH VARIABLES !C Do ielva = 1, kelva !C !C***READ A NEW BLOCK OF EQUATIONS - IF NEEDED !c ! PRELIMINARY THEORY AND STANDARD SUBROUTINES ! If (nbufa/=0) Goto 412 !c??? BACKSPACE 2 Read (2) equat, eqrhs, npivo, namev !c??? BACKSPACE 2 nbufa = mbufa If (kresl==1) Goto 412 !c??? BACKSPACE 4 Read (4) eqrhs !c??? BACKSPACE 4 412 Continue !c !c PREPARE TO BACK-SUBSTITUTE FROM THE CURRENT EQUATION !c ifron = npivo(nbufa) nikno = namev(nbufa) pivot = equat(ifron, nbufa) If (iffix(nikno)/=0) vecrv(ifron) = fixed(nikno) If (iffix(nikno)==0) equat(ifron, nbufa) = 0.0 !c !c*** BACK-SUBSTITUTE IN THE CURRENT EQUATION !c Do jfron = 1, mfron eqrhs(nbufa) = eqrhs(nbufa) - vecrv(jfron)*equat(jfron, nbufa) End Do !C !C*** PUT THE FINAL VALUES WHERE THEY BELONG !c If (iffix(nikno)==0) vecrv(ifron) = eqrhs(nbufa)/pivot If (iffix(nikno)/=0) fixed(nikno) = -eqrhs(nbufa) nbufa = nbufa - 1 asdis(nikno) = vecrv(ifron) End Do !C !C*** ADD DISPLACEMENTS TO PREVIOUS TOTAL VALUES !C Do itotv = 1, ntotv tdisp(itotv)=tdisp(itotv) + asdis(itotv) End Do !C !C*** STORE REACTIONS FOR PRINTING LATER !c kboun = 1 Do ipoin = 1, npoin nloca = (ipoin-1)*ndofn Do idofn = 1, ndofn ngush = nloca + idofn If (iffix(ngush)>0) Goto 360 End Do Goto 370 360 Do idofn = 1, ndofn ngash = nloca + idofn treac(kboun, idofn) = treac(kboun, idofn) + fixed(ngash) End Do kboun = kboun + 1 370 End Do !c !c*** ADD REACTIONS INTO THE TOTAL LOAD ARRAY !C Do ipoin = 1, npoin Do ielem = 1, nelem Do inode = 1, nnode nloca = iabs(lnods(ielem,inode)) If(ipoin==nloca) Goto 720 End Do End Do 720 Do idofn=1, ndofn ngash=(inode-1)*ndofn+idofn mgash=(ipoin-1)*ndofn+idofn tload(ielem,ngash)=tload(ielem,ngash)+fixed(mgash) End Do ! !c FINITE ELEMENTS IN PLASTICITY ! End Do Return 900 Format('0',3X,'NEGATIVE OR ZERO PIVOT ENCOUNTERED FOR VARIABLE & NO.','I4','OF VALUE',E17.6) End Subroutine front !c****************************************** !******************************************* ! RESMP !******************************************* Subroutine resmp(asdis, coord, effst, eload, epstn, lnods,& matno, melem, mmats, mpoin, mtotg, mtotv,& ncrit, nelem, nevab, ngaus, nnode, props,& strsg) !c******************************************************************** !c !c*** EVALUATES EQUIVALENT NODAL FORCES FOR THE STRESS RESULTANTS !c*** IN MINDLIN PLATES DURING EP ANALYSIS !c !c******************************************************************** Dimension asdis(mtotv), avect(5), cartd(2, 9), & coord(mpoin, 2), deriv(2, 9), desig(5), devia(4),& dvect(5),& effst(mtotg), elcod(2, 9),& eldis(3, 9), eload(melem, 27), epstn(mtotg), gpcod(2, 9),& lnods(melem, 9), matno(melem), posgp(4), & props(mmats, 8), sgtot(5), shape(9), sigma(5),& stres(5), strsg(5, mtotg), weigp(4), & dflex(3, 3), dsher(2, 2), bflei(3, 3), bshei(2, 3), & dummy(3, 3), force(3), dgrad(8) ntime = 1 Do ielem = 1, nelem Do ievab = 1, nevab eload(ielem, ievab) = 0.0 End Do End Do kgaus = 0 lgaus = 0 Do ielem = 1, nelem lprop = matno(ielem) !C !C*** COMPUTE COORDINATE AND INCREMENTAL DISPLACEMENTS OF THE !C ELEMENT NODAL POINTS !C Do inode = 1, nnode lnode = iabs(lnods(ielem,inode)) nposn = (lnode-1)*3 Do idofn = 1, 3 nposn = nposn + 1 eldis(idofn, inode) = asdis(nposn) End Do Do idime = 1, 2 elcod(idime, inode) = coord(lnode, idime) End Do End Do kcasp = 0 Call modpb(dflex, dummy, dsher, lprop, mmats, props,& 0, 1, 1) Call gaussq(ngaus, posgp, weigp) Do igaus = 1, ngaus Do jgaus = 1, ngaus bring = 1.0 kgaus = kgaus + 1 exisp = posgp(igaus) etasp = posgp(jgaus) Call sfr2(deriv, etasp, exisp, nnode, shape) kgasp = kgasp + 1 Call jacob2(cartd, deriv, djacb, elcod, gpcod, ielem,& kgasp, nnode, shape) ! !c***** ELASTO-PLASTIC MlNDLlN PLATE BENDING ANALYSIS ! darea = djacb*weigp(igaus)*weigp(jgaus) Call gradmp(cartd, dgrad, eldis, 3, nnode) Call strmp (cartd, dflex, dgrad, dsher, eldis, nnode,& shape, stres, 1, 0) preys = props(lprop, 6) + epstn(kgaus)*props(lprop, 7) Do istre = 1, 3 desig(istre) = stres(istre) sigma(istre) = strsg(istre, kgaus) + stres(istre) End Do Call invmp(devia, ncrit, sint3, steff, sigma, theta,& varj2, yield) espre = effst(kgaus) - preys If (espre>=0.0) Goto 50 escur = yield - preys If (escurtdler) nchek = l If (ncdis(4)==0) tdito = tdito 50 Continue Write (6, 600) Write (6, 601)(tdidf(idofn), idofn=1, ndofn) Write (6, 602) Write (6, 603) tdito !C**** RESIDUAL CONVERGENCE CHECK 1000 If (ifres==0) goto 2000 !C**** ASSEMBLE TOTAL AND RESIDUAL FORCE VECTORS Do itotv = 1, ntotv refor(itotv) = 0.0 tofor(itotv) = 0.0 End Do Do ielem = 1, nelem kevab = 0 Do inode = 1, nnode locno = iabs(lnods(ielem,inode)) Do idofn = 1, ndofn kevab = kevab + 1 nposi = (locno-1)*ndofn + idofn tofor(nposi) = tofor(nposi) + tload(ielem, kevab) refor(nposi) = refor(nposi) + eload(ielem, kevab) End Do End Do End Do !C*** COMPUTE TOTAL AND DIRECTIONAL NORMS OF RESIDUAL AND TOTAL FORCE refto = 0.0 tofto = 0.0 Call vzero(ndofn, refdf) Call vzero(ndofn, tofdf) nposi = 0 Do ipoin = 1, npoin Do idofn = 1, ndofn nposi = nposi + 1 refdf(idofn) = refdf(idofn) + refor(nposi)*refor(nposi) tofdf(idofn) = tofdf(idofn) + tofor(nposi)*tofor(nposi) End Do End Do Do idofn = 1, ndofn refto = refto + refdf(idofn) tofto = tofto + tofdf(idofn) refdf(idofn) = sqrt(refdf(idofn)) tofdf(idofn) = sqrt(tofdf(idofn)) End Do refto = sqrt(refto) tofto = sqrt(tofto) !C*** CHECK FOR CONVERGENCE AND PRINT ERRORS PER CENT Do idofn = 1, ndofn If (tofdf(idofn)==0.0) goto 90 tofdf(idofn) = 100.*refdf(idofn)/tofdf(idofn) If (ncres(idofn)/=0 .And. tofdf(idofn)>toler) nchek = 1 If (ncres(idofn)==0) tofdf(idofn) = -tofdf(idofn) 90 End Do If (tofto == 0.0) goto 100 tofto = 100.*refto/tofto If (ncres(4) /= 0 .And. tofto > toler) nchek = 1 If (ncres(4) == 0) tofto = -tofto 100 Continue Write (6, 604) Write (6, 601)(tofdf(idofn), idofn = 1, ndofn) !c*****FINITE ELEMENTS IN PLASTICITY Write (6, 602) Write (6, 603) tofto !c !c***** PRINT CONVERGENCE CODE !c 2000 Write (6, 605) nchek Return 606 Format (///, 'In Conver',10X,'Iteration number',I3,/) 600 Format (1X, 'DISPLACEMENT CHANGE NORMT', //) 601 Format (1X, 5(E10.3,5X)) 602 Format (5X, 'TOTAL') 603 Format (3X, E10.3) 604 Format (1X, 'RESIDUAL NORM', //) 605 Format (1X, 'CONVERGENCE CODE', I4, //) End Subroutine convmp !************************************* !************************************* ! OUTMPA?????? !************************************* ! 20230427 !************************************* Subroutine outmpa(epstn, iiter, mtotg, mtotv, mvfix, nelem, ngaus, nlaps, nofix, noutp, npoin, nvfix, strsg, tdisp, treac) !C*********************************************************************** !C !C****** OUTPUT DISPLACEMENTS,REACTIONS AND GAUSS POINT STRESSES !C****** IN EACH LAYER FOR EP MINDLIN PLATE ANALYSIS !C !C************************************************************************ Dimension epstn(mt0tg), gpcod(2, 9), nofix(mvf1x), noutp(2), strsg(5, tvotg), tdisp(mtotv), treac(mvfix, 3) !c !C***** OUTPUT DISPLACEMENTS !C If (koutp8) Goto 30 s2 = s*2.0 t2 = t*2.0 ss = s*s tt = t*t st = s*t sst = s*s*t stt = s*t*t st2 = s*t*2.0 ! c ! c***Shape functions for 8 noded element. ! c shape(1) = (-1.0+st+ss+tt-sst-stt)/4.0 shape(2) = (1.0-t-ss+sst)/2.0 shape(3) = (-1.0-st+ss+tt-sst+stt)/4.0 shape(4) = (1.0+s-tt-stt)/2.0 shape(5) = (-1.0+st+ss+tt+sst+stt)/4.0 shape(6) = (1.0+t-ss-sst)/2.0 shape(7) = (-1.0-st+ss+tt+sst-stt)/4.0 shape(8) = (1.0-s-tt+stt)/2.0 ! c ! c***Shape function derivertives. ! c deriv(1, 1) = (t+s2-st2-tt)/4.0 deriv(1, 2) = -s + st deriv(1, 3) = (-t+s2-st2+tt)/4.0 deriv(1, 4) = (1.0-tt)/2.0 deriv(1, 5) = (t+s2+st2+tt)/4.0 deriv(1, 6) = -s - st deriv(1, 7) = (-t+s2+st2-tt)/4.0 deriv(1, 8) = (-1.0+tt)/2.0 deriv(2, 1) = (s+t2-ss-st2)/4.0 deriv(2, 2) = (-1.0+ss)/2.0 deriv(2, 3) = (-s+t2-ss+st2)/4.0 deriv(2, 4) = -t - st deriv(2, 5) = (s+t2+ss+st2)/4.0 deriv(2, 6) = (1.0-ss)/2.0 deriv(2, 7) = (-s+t2+ss-st2)/4.0 deriv(2, 8) = -t + st 30 If (nnode==8) Return !C*** BUBBLE FUNCTION FOR HIERARCHICAL AND HETEROSIS ELEMENTS shape(9) = (1.0-ss)*(1.0-tt) deriv(1, 9) = -s2*(1.0-tt) deriv(2, 9) = -t2*(1.0-ss) Return End Subroutine sfr2 !***************************************************************** !***************************************************************** !c jacob2 !***************************************************************** Subroutine jacob2(cartd, deriv, djacb, elcod, gpcod, ielem, kgasp, & nnode, shape) !C************************************************************************ !C !C**** THIS SUBROUTINE EVALUATES THE JACOBIAN MATRIX AND THE CARTESIAN !C SHAPE FUNCTION DERIVATIVES !C !C************************************************************************ Dimension cartd(2, 9), deriv(2, 9), elcod(2, 9), gpcod(2, 9), shape(9),& xjaci(2, 2), xjacm(2, 2) !c !c*****CALCULATE COORDINATES OF SAMPLING POINT !c Do idime = 1, 2 gpcod(idime, kgasp) = 0.0 Do inode = 1, nnode gpcod(idime, kgasp) = gpcod(idime, kgasp) + elcod(idime, inode)& *shape(inode) End Do End Do !C !C*** CREATE JACOBIAN MATRIX XJACM !C Do idime = 1, 2 Do jdime = 1, 2 xjacm(idime, jdime) = 0.0 Do inode = 1, nnode xjacm(idime, jdime) = xjacm(idime, jdime) + deriv(idime, inode)*& elcod(jdime, inode) End Do End Do End Do !C !C**** CALCULATE DETERMINANT AND INVERSE OF JACOBIAN MATRIX !C djacb = xjacm(1, 1)*xjacm(2, 2) - xjacm(1, 2)*xjacm(2, 1) If (djacb) 6, 6, 8 6 Write (6, 600) ielem Stop 8 Continue xjaci(1, 1) = xjacm(2, 2)/djacb xjaci(2, 2) = xjacm(1, 1)/djacb xjaci(1, 2) = -xjacm(1, 2)/djacb xjaci(2, 1) = -xjacm(2, 1)/djacb !c !c**** CALCULATE CARTESIAN DERIVATIVES !c Do idime = 1, 2 Do inode = 1, nnode cartd(idime, inode) = 0.0 Do jdime = 1, 2 cartd(idime, inode) = cartd(idime, inode) + xjaci(idime, jdime)*& deriv(jdihe, inode) !c !c PRELIMINARY THEORY AND STANDARD SUBROUTINES !c End Do End Do End Do Return 600 Format (//, ' PROGRAM HALTED IN SUBROUTINE JACOB2', /,11X, & 'ZERO OR NEGATIVE AREA',/,10X, 'ELEMENT NUMBER',I5) End Subroutine jacob2 !******************************************************************** !******************************************************************** ! modpb !******************************************************************** Subroutine modpb(dflex, dplan, dsher, lprop, mmats, props,& ifpla, iffle, ifshe) !c************************************************************* !c !c****CALCULATES MATRIX OF ELASTIC RIGIDITIES !c****FOR MINDLIN PLATE !c !c************************************************************* Dimension dflex(3, 3), dplan(3, 3), dsher(2, 2),& props(mmats, 8) young = props(lprop, 1) poiss = props(lprop, 2) thick = props(lprop, 3) !C*** FORM DPLAN If (ifpla==0) Goto 10 Do irows = 1, 3 Do jcols = 1, 3 dplan(irows, jcols) = 0.0 End Do End Do const = (young*thick)/(1.0-poiss*poiss) dplan(1, 1) = const dplan(2, 2) = const dplan(1, 2) = const*poiss dplan(2, 1) = const*poiss dplan(3, 3) = const*(1.0-poiss)/2.0 !c**** FORM DFLEX 10 If (ifflex==0) Goto 20 Do irows = 1, 3 Do jcols = 1, 3 dflex(irows, jcols) = 0.0 End Do End Do const = (young*thick**3)/(12.*(1.0-poiss*poiss)) dflex(1, 1) = const dflex(2, 2) = const dflex(1, 2) = const*poiss dflex(2, 1) = const*poiss dflex(3, 3) = const*(1.-poiss)/2 !C**** FORM DSHER 20 If (ifshe==0) Return Do irows = 1, 2 Do jcols = 1, 2 dsher(irows, jcols) = 0.0 End Do End Do dsher(1, 1) = (young*thick)/(2.4+2.4*poiss) dsher(2, 2) = (young*thick)/(2.4+2.4*poiss) Return End Subroutine modpb !********************************************************************** !********************************************************************** ! nodexy !*********************************************************************** Subroutine nodexy(coord, lnods, melem, mpoin, nelem, nnode) !c********************************************************************** !c !c****This subroutine interpolates the mid side nodes of straight !c sides of elements and the central node of 9 noded elements. !c !c********************************************************************** Dimension coord(mpoin, 2), lnods(melem, 9) If (nnode==4) Return !c !c***Loop over each element. !c Do ielem = 1, nelem !c !c***Loop over each element edge. !c nnod1 = 9 If (nnode==8) nnod1 = 7 Do inode = 1, nnod1, 2 If (inode==9) Goto 50 !c !c***Compute the node number of the first node. !c nodst = lnods(ielem, inode) igash = inode + 2 If (igash>8) igash = 1 !c !c***Compute the node number of the last node. !c nodfn = lnods(ielem, igash) midpt = inode + 1 !c !c***Compute the node number of the intermediate node. !c nodmd = lnods(ielem, midpt) total = abs(coord(nodmd,1)) + abs(coord(nodmd,2)) !c !c***If the coordinates of the intermediate node are both !c zero interpolate by a straight line. !c If (total>0.0) Goto 20 kount = 1 10 coord(nodmd, kount) =& (coord(nodst,kount)+coord(nodfn,kount))/2.0 kount = kount + 1 If (kount==2) Goto 10 20 End Do Goto 30 50 lnode = lnods(ielem, inode) total = abs(coord(lnode,1)) + abs(coord(lnode,2)) If (total>0.0) Goto 30 lnod1 = lnods(ielem, 1) lnod3 = lnods(ielem, 3) lnod5 = lnods(ielem, 5) lnod7 = lnods(ielem, 7) kount = 1 40 coord(lnode, kount) = (coord(lnod1,kount)+coord(lnod3,kount)+& coord(lnod5,kount)+coord(lnod7,kount))/4.0 kount = kount + 1 If (kount==2) Goto 40 30 End Do Return End Subroutine nodexy !************************************************************ !************************************************************ ! flowmp !*********************************************************** Subroutine flowmp(abeta, avect, devia, dmatx, dvect, hards, & ncrit, sint3, steff, theta, varj2) !c************************************************************ !c !C*** DETERMINES YIELD FUNCTION DERIVATIVES FOR MINDLIN PLATES !c*** 1 Von Mises !c*** 2 Tresca !c !c ELASTO-PLASTIC MINDLIN PLATE BENDING ANALYSIS !C !C************************************************************* Dimension avect(5), devia(4), dmatx(3, 3), dvect(5),& veca1(3), veca2(3), veca3(3) !c !C*** DETERMINE THE VECTOR DERIVATIVE OF F FOR VON-MISES !c sinth = sin(theta) costh = cos(theta) root3 = 1.73205080757 !c !C*** CALCULATE VECTOR A1 !C veca1(1) = 0.333333333333 veca1(2) = 0.333333333333 veca1(3) = 0.0 !C !C*** CALCULATE VECTOR A2 !C Do istre = 1, 3 veca2(istre) = devia(istre)/(2.0*steff) End Do veca2(3) = devia(3)/steff !c !C*** CALCULATE VECTOR A3 !c veca3(1) = devia(2)*devia(4+varj2/3.0) veca3(2) = devia(1)*devia(4+varj2/3.0) veca3(3) = -2.0*devia(3)*devia(4) Goto (1, 2) ncrit !C !C*** VON MISES !C 1 cons1 = 0.0 cons2 = root3 cons3 = 0.0 Goto 40 !C !C*** TRESCA !C 2 cons1 = 0.0 abthe = abs(theta*57.29577951308) If (abthe

!************************************************* !************************************************* !************************************************* ! ELASTO-PLASTIC MINDLIN PLATE BENDING ANALYSIS ! Overall structure of MINDLIN ! Program FEMP(Input,Output,Tape5=Input,Tape6=Output, ! Tape1,Tape2,Tape3,Tape4,Tape8,Tape9) !************************************************* ! 有限元法明德林层板非分层弯曲力学分析 !********************************************************************* !********************************************************************* !********************************************************************* ! SINOMACH ! ! FORTRAN95,2013&2018 version ! ! ! PotsyYZ,周勇,Speagaul,XiaolG,Fiddler,ChengnW,Johnny,SwanseaUK ! ! 20230424 !********************************************************************* Program main !c******************************************************************** !c !c*** ELASTO-PLASTIC ANALYSIS OF NON-LAYERED MINDLIN PLATES USING !C*** 4-,8- , 9-NODED OR HETEROSIS ISOPARAMETRIC QUADRILATERALS !C !C******************************************************************* Dimension asdis(240), coord(80, 2), effst(225), eload(25, 27),& epstn(225), estif(27, 27), & eqrhs(10), equat(40, 10), fixed(240),& iffix(240), gload(40), gstif(860), lnods(25, 9), locel(27),& matno(25), nacva(40), namev(10), ncdis(4), ncres(4),& ndest(27), ndfro(25), nofix(40), noutp(2), npivo(10),& posgp(4), presc(40,3), props(10,8), refor(240),& rload(25,27), strsg(5, 225), tofor(240),& tdisp(240), tload(25, 27), treac(40, 3), vecrv(40),& weigp(4) !C !C*** PRESET VARIABLES ASSOCIATED WITH DYNAMIC DIMENSIONS !C Call dimmp(mbufa, melem, mevab, mfron, mmats, mpoin,& mstif, mtotg, mtotv, mvfix, ndime, ndofn,& nprop, nstre) !C !C*** CALL THE SUBROUTINE WHICH READS MOST OF THE PROBLEM DATA !c Call input(coord,iffix,lnods,matno,melem,mevab,& mfron,mmats,mpoin,mtotv,mvfix,nalgo,& ncrit,ndfro,ndime,ndofn,nelem,nevab,& ngaus,nlaps,nincs,nmats,nnode,nofix,& npoin,nprop,nstre,nstrl,nswit,ntotg,& ntotv,ntype,nvfix,posgp,presc,props,& weigp) !c !C**** INITIALIZE ARRAYS TO ZERO !C Call zeromp(effst,eload,epstn,melem,mevab,mtotg,& mtotv,mvfix,ndofn,nelem,nevab,ngaus,& ntotg,ntotv,nvfix,strsg,tdisp,tfact,& tload,treac) !C !C**** !C Call mindpb(ifdis,iffix,ifres,lnods,melem,mtotv,& ncdis,ncres,nelem,ntype) !C !C !C !C*** CCMPUTE LOAD AFTER READING RELEVANT EXTRA DATA !C Call loadpb(coord,lnods,matno,melem,mmats,mpoin,& nelem,nevab,ngaus,nnode,npoin,props,& rload) !C !C*** LOOP OVER EACH INCREMENT !C !c ELASTO-PLASTIC MINDLIN PLATE BENDING ANALYSIS !c Do iincs=1, nincs !c !c*** READ DATA FOR CURRENT INCREMENT !c Call increm(eload, fixed, iincs, melem, mevab, miter, & mtotv, mvfix, ndofn, nelem, nevab, noutp, & nofix, ntotv, nvfix, presc, rload, tfact, & tload, toler) !C !C*** LOOP OVER EACH ITERATION !C Do iiter = 1, miter !c !C*** CALL ROUTINE WHICH SELECTS SOLUTION ALGORITHM VARIABLE KRESL !c Call algor(fixed, iincs, iiter, kresl, mtotv, nalgo,& ntotv) !C !C*** CHECK WHETHER A NEW EVALUATION OF THE STIFFNESS MATRICES IS NEEDED. !c If (kresl==1) & Call stifmp(coord, epstn, iincs, lnods, matno, melem,& mevab, mmats, mpoin, mtotg, ncrit, nelem,& nevab, ngaus, nnode, props, strsg) !C !C*** SOLVE EQUATIONS !C Call front(asdis, eload, eqrhs, equat, estif, fixed, & iffix, iincs, iiter, gload, gstif, kresl, & lnods, locel, mbufa, melem, mevab, mfron, & mstif, mtotv, mvfix, nacva, namev, ndest, & ndofn, nelem, nevab, nnode, nofix, npivo, & npoin, ntotv, tdisp, tload, treac, vecrv) !c !C*** CALCULATE RESIDUAL FORCES !C Call resmp(asdis, coord, effst, eload, epstn, lnods,& matno, melem, mmats, mpoin, mtotg, mtotv,& ncrit, nelem, nevab, ngaus, nnode, props,& strsg) !C !C*** CHECK FOR CONVERGENCE !C Call convmp(asdis, eload, iiter, ifdis, ifres, lnods,& melem, mevab, mtotv, nchek, ncdis, ncres,& ndofn, nelem, nevab, nnode, npoin, ntotv,& refor, tofor, tdisp, tload, toler) !C !C*** OUTPUT RESULTS IF REQUIRED !C If (iiter==1 .And. noutp(1)>0) & Call outmp(epstn, iiter, mtotg, mtotv, mvfix, nelem, & ngaus, nofix, noutp, npoin, nvfix, strsg, & tdisp, treac) !c !c*** IF SOLUTION HAS CONVERGED STOP ITERATING AND OUTPUT RESULTS !C If (nchek==0) Goto 100 End Do !c !c****FINITE ELEMENTS IN PLASTICITY !c If (nalog==2) Goto 100 Stop 100 Call outmp(epstn, iiter, mtotg, mtotv, mvfix,nelem, & ngaus, nofix, noutp, npoin, nvfix,strsg, & tdisp, treac) End Do Stop End Program main !*********************************************************** !*********************************************************** ! 20230423 !*********************************************************** ! DIMMP !*********************************************************** Subroutine dimmp(mbufa, melem, mevab, mfron, mmats, mpoin,& mstif, mtotg, mtotv, mvfix, ndime, ndofn,& nprop, nstre) !c************************************************************ !c !c****Sets up dynamic dimensions-must agree with dimensions !c****in femp !c !c************************************************************ mbufa=10 melem=25 mfron=40 mmats=10 mpoin=80 mstif=(mfron*mfron-mfron)/2.0+mfron mtotg=melem*9 ndofn=3 mototv=mpoin*ndofn mvfix=40 ndime=2 nprop=8 nstre=5 mevab=ndofn*9 Return End Subroutine dimmp !*************************************************************** !*************************************************************** ! input !*************************************************************** Subroutine input(coord,iffix,lnods,matno,melem,mevab,mfron,mmats,& mpoin,mtotv,mvfix,nalgo,& ncrit,ndfro,ndofn,nelem,& nevab,ngaus,ngau2,& nincs,nmats,nnode,nofix,npoin,nprop,nstre,nstr1,& ntotg,ntotv,ntype,nvfix,posgp,presc,props,weigp) !c********************************************************************* !c !c****This subroutine accepts most of the input data !c !c********************************************************************* Dimension coord(mpoin,2), iffix(mtotv), lnods(melem,9), & matno(melem), ndfro(melem), & nofix(mvfix), posgp(4), presc(mvfix,ndofn), & props(mmats,nprop), title(12), weigp(4) ! Rewind 1 ! Rewind 2 ! Rewind 3 ! Rewind 4 ! Rewind 8 Read(5,920) title Write(6,920) title !c !c***Read the first data card, and echo it immediately. !c Read(5,900) npoin, nelem, nvfix, ntype, nnode, nmats, ngaus,& nalgo, ncrit, nincs, nstre nevab=ndofn*nnode nstr1=nstre+1 If(ntype==3) nstr1=nstre ntotv=npoin*ndofn ngau2=ngaus*ngaus ntotg=nelem*ngau2 Write(6,901) npoin, nelem, nvfix, ntype, nnode, nmats, ngaus, nevab,& nalgo, ncrit, nincs, nstre Call check1(ndofn,nelem,ngaus,nmats,nnode,npoin,& nstre,ntype,nvfix,ncrit,nalgo,nincs) !c !c***Read the element nodal connections, and the property numbers. !c Write(6,902) Do ielem=1, nelem Read(5,900) numel, matno(numel), (lnods(numel,inode),inode=1,nnode) Write(6,903) numel, matno(numel), (lnods(numel,inode),inode=1,nnode) End Do !c !c***Zero all the nodal coordinates, prior to reading some of them. !c Do ipoin=1, npoin Do idime=1, 2 coord(ipoin,idime)=0.0 End Do End Do !c !c***Read some nodal coordinates, finishing with the node of all. !c Write(6,904) 6 Read(5,905) ipoin, (coord(ipoin,idime),idime=1,2) If(ipoin/=npoin) Goto 6 !c !c***Interpolate coordinates of mid-side nodes !c Call nodexy(coord,lnods,melem,mpoin,nelem,nnode) Do ipoin=1, npoin Write(6,906) ipoin, (coord(ipoin,idime),idime=1,2) End Do !c !c***Read the fixed values. !c Write(6,907) Do ivfix=1, nvfix Read(5,908) nofix(ivfix), ifpre, (presc(ivfix,idofn),idofn=1,ndofn) Write(6,908) nofix(ivfix), ifpre, (presc(ivfix,idofn),idofn=1,ndofn) nloca=(nofix(ivfix)-1)*ndofn ifdof=10**(ndofn-1) Do idofn=1, ndofn ngash=nloca+idofn If(ifpre1) Goto 455 Do ipoin = 1, npoin klast = 0 Do ielem = 1, nelem Do inode = 1, nnode If (lnods(ielem,inode)/=ipoin) Goto 120 klast = ielem nlast = inode 120 End Do ! !c FINITE ELEMENTS IN PLASTICITY ! End Do If (klast/=0) lnods(klast, nlast) = -ipoin End Do 455 Continue !C !c*** START BY INITIALIZING EVERYTHING THAT MATTERS TO ZERO !c Do ibufa = 1, mbufa eqrsh(ibufa) = 0.0 End Do Do istif = 1, mstif gstif(istif) = 0.0 End Do Do ifron = 1, mfron gload(ifron) = 0.0 vecrv(ifron) = 0.0 nacva(ifron) = 0 Do ibufa = 1, mbufa epuat(ifron, ibufa) = 0.0 End Do End Do !C !C*** AND PREPARE FOR DISC READING AND WRITING OPERATIONS !c nbufa = 0 If (kresl>1) nbufa = mbufa !? REWIND 1 !? REWIND 2 !? REWIND 3 !? REWIND 4 !? REWIND 8 !C !C*** ENTER MAIN ELEMENT ASSEMBLY-REDUCTION LOOP !C nfron = 0 kelva = 0 Do ielem = 1, nelem If (kresl>1) Goto 400 kevab = 0 Read (1) estif Do inode = 1, nnode Do idofn = 1, ndofn nposi = (inode-1)*ndofn + idofn locno = lnods(ielem, inode) If (locno>0) locel(nposi) = (locno-1)*ndofn + idofn If (locnonfron) nfron = ndest(kevab) 210 End Do Write (8) locel, ndest, nacva, nfron 400 If (kresl>1) Read (8) locel, ndest, nacva, nfron !c !C**** ASSEMBLE ELEMENT LOADS !c Do ievab = 1, nevab idest = ndest(ievab) gload(idest) = gload(idest) + eload(ielem, ievab) !C !C*** ASSEMBLE THE ELEMENT STIFFNESSES-BUT NOT IN RESOLUTION !C If (kresl>1) Goto 402 Do jevab = 1, ievab jdest = ndest(jevab) ngash = nfunc(idest, jdest) ngish = nfunc(jdest, idest) If (jdest>=idest) gstif(ngash) = gstif(ngash) + estif(ievab, jevab) If (jdest0.0) Goto 235 Write (6, 900) nikno, pivot Stop 235 Continue equat(ifron, nbufa) = 0.0 !c !C*** ENQUIRE WHETHER PRESENT VARIABLE IS FREE OR PRESCRIBED !C If (iffix(nikno)==0) Goto 250 !c !C*** DEAL WITH A PRESCRIBED DEFLECTION !C Do jfron = 1, nfron gload(jfron) = gload(jfron) - fixed(nikno)*equat(jfron, nbufa) End Do Goto 280 !c !c*** ELIMINATE A FREE VARIABLE - DEAL WITH THE RIGHT HAND SIDE FIRST !C 250 Do jfron = 1, nfron gload(jfron) = gload(jfron) - equat(jfron, nbufa)*eqrhs(nbufa)/pivot !c !c**** NOW DEAL WITH THE COEFFICIENTS IN CORE !C If (kresl>1) Goto 418 If (equat(jfron,nbufa)==0.0) Goto 270 nloca = nfunc(0, jfron) cureq = equat(jfron, nbufa) Do lfron = 1, jfron ngash = lfron + nloca gstif(ngash) = gstif(ngash) - cureq*equat(lfron, nbufa)& /pivot End Do 418 Continue 270 End Do 280 equat(ifron, nbufa) = pivot !C !C*** RECORD THE NEW VACANT SPACE, AND REDUCE FRONTWIDTH IF POSSIBLE !C nacva(ifron) = 0 Goto 290 !c !C*** COMPLETE THE ELEMENT LOOP IN THE FORWARD ELIMINATION !c 300 End Do 290 If (nacva(nfron)/=0) Goto 310 nfron = nfron - 1 If (nfron>0) Goto 290 310 End Do End Do If (kresl==1) Write (2) equat, eqrhs, npivo, namev ! !c?? BACKSPACE 2 ??? !c !c*** ENTER BACKSUBSTITUTION PHASE. LOOP BACKWARDS THROUGH VARIABLES !C Do ielva = 1, kelva !C !C***READ A NEW BLOCK OF EQUATIONS - IF NEEDED !c ! PRELIMINARY THEORY AND STANDARD SUBROUTINES ! If (nbufa/=0) Goto 412 !c??? BACKSPACE 2 Read (2) equat, eqrhs, npivo, namev !c??? BACKSPACE 2 nbufa = mbufa If (kresl==1) Goto 412 !c??? BACKSPACE 4 Read (4) eqrhs !c??? BACKSPACE 4 412 Continue !c !c PREPARE TO BACK-SUBSTITUTE FROM THE CURRENT EQUATION !c ifron = npivo(nbufa) nikno = namev(nbufa) pivot = equat(ifron, nbufa) If (iffix(nikno)/=0) vecrv(ifron) = fixed(nikno) If (iffix(nikno)==0) equat(ifron, nbufa) = 0.0 !c !c*** BACK-SUBSTITUTE IN THE CURRENT EQUATION !c Do jfron = 1, mfron eqrhs(nbufa) = eqrhs(nbufa) - vecrv(jfron)*equat(jfron, nbufa) End Do !C !C*** PUT THE FINAL VALUES WHERE THEY BELONG !c If (iffix(nikno)==0) vecrv(ifron) = eqrhs(nbufa)/pivot If (iffix(nikno)/=0) fixed(nikno) = -eqrhs(nbufa) nbufa = nbufa - 1 asdis(nikno) = vecrv(ifron) End Do !C !C*** ADD DISPLACEMENTS TO PREVIOUS TOTAL VALUES !C Do itotv = 1, ntotv tdisp(itotv)=tdisp(itotv) + asdis(itotv) End Do !C !C*** STORE REACTIONS FOR PRINTING LATER !c kboun = 1 Do ipoin = 1, npoin nloca = (ipoin-1)*ndofn Do idofn = 1, ndofn ngush = nloca + idofn If (iffix(ngush)>0) Goto 360 End Do Goto 370 360 Do idofn = 1, ndofn ngash = nloca + idofn treac(kboun, idofn) = treac(kboun, idofn) + fixed(ngash) End Do kboun = kboun + 1 370 End Do !c !c*** ADD REACTIONS INTO THE TOTAL LOAD ARRAY !C Do ipoin = 1, npoin Do ielem = 1, nelem Do inode = 1, nnode nloca = iabs(lnods(ielem,inode)) If(ipoin==nloca) Goto 720 End Do End Do 720 Do idofn=1, ndofn ngash=(inode-1)*ndofn+idofn mgash=(ipoin-1)*ndofn+idofn tload(ielem,ngash)=tload(ielem,ngash)+fixed(mgash) End Do ! !c FINITE ELEMENTS IN PLASTICITY ! End Do Return 900 Format('0',3X,'NEGATIVE OR ZERO PIVOT ENCOUNTERED FOR VARIABLE & NO.','I4','OF VALUE',E17.6) End Subroutine front !c****************************************** !******************************************* ! 20230423 !******************************************* ! RESMP !******************************************* Subroutine resmp(asdis, coord, effst, eload, epstn, lnods,& matno, melem, mmats, mpoin, mtotg, mtotv,& ncrit, nelem, nevab, ngaus, nnode, props,& strsg) !c******************************************************************** !c !c*** EVALUATES EQUIVALENT NODAL FORCES FOR THE STRESS RESULTANTS !c*** IN MINDLIN PLATES DURING EP ANALYSIS !c !c******************************************************************** Dimension asdis(mtotv), avect(5), cartd(2, 9), & coord(mpoin, 2), deriv(2, 9), desig(5), devia(4),& dvect(5),& effst(mtotg), elcod(2, 9),& eldis(3, 9), eload(melem, 27), epstn(mtotg), gpcod(2, 9),& lnods(melem, 9), matno(melem), posgp(4), & props(mmats, 8), sgtot(5), shape(9), sigma(5),& stres(5), strsg(5, mtotg), weigp(4), & dflex(3, 3), dsher(2, 2), bflei(3, 3), bshei(2, 3), & dummy(3, 3), force(3), dgrad(6) ntime = 1 Do ielem = 1, nelem Do ievab = 1, nevab eload(ielem, ievab) = 0.0 End Do End Do kgaus = 0 lgaus = 0 Do ielem = 1, nelem lprop = matno(ielem) !C !C*** COMPUTE COORDINATE AND INCREMENTAL DISPLACEMENTS OF THE !C ELEMENT NODAL POINTS !C Do inode = 1, nnode lnode = iabs(lnods(ielem,inode)) nposn = (lnode-1)*3 Do idofn = 1, 3 nposn = nposn + 1 eldis(idofn, inode) = asdis(nposn) End Do Do idime = 1, 2 elcod(idime, inode) = coord(lnode, idime) End Do End Do kcasp = 0 Call modpb(dflex, dummy, dsher, lprop, mmats, props,& 0, 1, 1) Call gaussq(ngaus, posgp, weigp) Do igaus = 1, ngaus Do jgaus = 1, ngaus bring = 1.0 kgaus = kgaus + 1 exisp = posgp(igaus) etasp = posgp(jgaus) Call sfr2(deriv, etasp, exisp, nnode, shape) kgasp = kgasp + 1 Call jacob2(cartd, deriv, djacb, elcod, gpcod, ielem,& kgasp, nnode, shape) ! !c***** ELASTO-PLASTIC MlNDLlN PLATE BENDING ANALYSIS ! darea = djacb*weigp(igaus)*weigp(jgaus) Call gradmp(cartd, dgrad, eldis, 3, nnode) Call strmp (cartd, dflex, dgrad, dsher, eldis, nnode,& shape, stres, 1, 0) preys = props(lprop, 6) + epstn(kgaus)*props(lprop, 7) Do istre = 1, 3 desig(istre) = stres(istre) sigma(istre) = strsg(istre, kgaus) + stres(istre) End Do Call invmp(devia, ncrit, sint3, steff, sigma, theta,& varj2, yield) espre = effst(kgaus) - preys If (espre>=0.0) Goto 50 escur = yield - preys If (escurtdler) nchek = l If (ncdis(4)==0) tdito = tdito 50 Continue Write (6, 600) Write (6, 601)(tdidf(idofn), idofn=1, ndofn) Write (6, 602) Write (6, 603) tdito !C**** RESIDUAL CONVERGENCE CHECK 1000 If (ifres==0) goto 2000 !C**** ASSEMBLE TOTAL AND RESIDUAL FORCE VECTORS Do itotv = 1, ntotv refor(itotv) = 0.0 tofor(itotv) = 0.0 End Do Do ielem = 1, nelem kevab = 0 Do inode = 1, nnode locno = iabs(lnods(ielem,inode)) Do idofn = 1, ndofn kevab = kevab + 1 nposi = (locno-1)*ndofn + idofn tofor(nposi) = tofor(nposi) + tload(ielem, kevab) refor(nposi) = refor(nposi) + eload(ielem, kevab) End Do End Do End Do !C*** COMPUTE TOTAL AND DIRECTIONAL NORMS OF RESIDUAL AND TOTAL FORCE refto = 0.0 tofto = 0.0 Call vzero(ndofn, refdf) Call vzero(ndofn, tofdf) nposi = 0 Do ipoin = 1, npoin Do idofn = 1, ndofn nposi = nposi + 1 refdf(idofn) = refdf(idofn) + refor(nposi)*refor(nposi) tofdf(idofn) = tofdf(idofn) + tofor(nposi)*tofor(nposi) End Do End Do Do idofn = 1, ndofn refto = refto + refdf(idofn) tofto = tofto + tofdf(idofn) refdf(idofn) = sqrt(refdf(idofn)) tofdf(idofn) = sqrt(tofdf(idofn)) End Do refto = sqrt(refto) tofto = sqrt(tofto) !C*** CHECK FOR CONVERGENCE AND PRINT ERRORS PER CENT Do idofn = 1, ndofn If (tofdf(idofn)==0.0) goto 90 tofdf(idofn) = 100.*refdf(idofn)/tofdf(idofn) If (ncres(idofn)/=0 .And. tofdf(idofn)>toler) nchek = 1 If (ncres(idofn)==0) tofdf(idofn) = -tofdf(idofn) 90 End Do If tofto == 0.0) goto 100 tofto = 100.*refto/tofto If (ncres(4) /= 0 .And. tofto > toler) nchek = 1 If (ncres(4) == 0) tofto = -tofto Write (6, 604) Write (6, 601)(tofdf(idofn), idofn = 1, ndofn) !c*****FINITE ELEMENTS IN PLASTICITY Write (6, 602) Write (6, 603) tofto !c !c***** PRINT CONVERGENCE CODE !c 2000 Write (6, 605) nchek Return 606 Format (///, 'In Conver',10X,'Iteration number',I3,/) 600 Format (1X, 'DISPLACEMENT CHANGE NORMT', //) 601 Format (1X, 5(E10.3,5X)) 602 Format (5X, 'TOTAL') 603 Format (3X, E10.3) 604 Format (1X, 'RESIDUAL NORM', //) 605 Format (1X, 'CONVERGENCE CODE', I4, //) End Subroutine convmp !************************************* !************************************* ! OUTMP !************************************* Subroutine outmp(epstn, iiter, mtotg, mtotv, mvfix, nelem,& ngaus, nofix, noutp, npoin, nvfix, strsg,& tdisp, treac) !c********************************************************** !C !C*** OUTPUT DISPLACEMENTS,REACTIONS AND GAUSS POINT STRESS !C*** RESULTANTS FOR EP MINDLIN PLATE ANALYSIS !c !C ELASTO-PLASTIC MINDLIN PLATE BENDING ANALYSIS !c !c******************************************************************* Dimension epstn(mtotg), gpcod(2, 9), nofix(mvfix), noutp(2), & strsg(5, mtotg), tdisp(mtotv), treac(mvfix, 3) koutp = noutp(1) If (iiter>1) koutp = noutp(2) !C !C*** OUTPUT DISPLACEMENTS !C If (koutp8) Goto 30 s2 = s*2.0 t2 = t*2.0 ss = s*s tt = t*t st = s*t sst = s*s*t stt = s*t*t st2 = s*t*2.0 ! c ! c***Shape functions for 8 noded element. ! c shape(1) = (-1.0+st+ss+tt-sst-stt)/4.0 shape(2) = (1.0-t-ss+sst)/2.0 shape(3) = (-1.0-st+ss+tt-sst+stt)/4.0 shape(4) = (1.0+s-tt-stt)/2.0 shape(5) = (-1.0+st+ss+tt+sst+stt)/4.0 shape(6) = (1.0+t-ss-sst)/2.0 shape(7) = (-1.0-st+ss+tt+sst-stt)/4.0 shape(8) = (1.0-s-tt+stt)/2.0 ! c ! c***Shape function derivertives. ! c deriv(1, 1) = (t+s2-st2-tt)/4.0 deriv(1, 2) = -s + st deriv(1, 3) = (-t+s2-st2+tt)/4.0 deriv(1, 4) = (1.0-tt)/2.0 deriv(1, 5) = (t+s2+st2+tt)/4.0 deriv(1, 6) = -s - st deriv(1, 7) = (-t+s2+st2-tt)/4.0 deriv(1, 8) = (-1.0+tt)/2.0 deriv(2, 1) = (s+t2-ss-st2)/4.0 deriv(2, 2) = (-1.0+ss)/2.0 deriv(2, 3) = (-s+t2-ss+st2)/4.0 deriv(2, 4) = -t - st deriv(2, 5) = (s+t2+ss+st2)/4.0 deriv(2, 6) = (1.0-ss)/2.0 deriv(2, 7) = (-s+t2+ss-st2)/4.0 deriv(2, 8) = -t + st If (nnode==8) Return !C*** BUBBLE FUNCTION FOR HIERARCHICAL AND HETEROSIS ELEMENTS shape(9) = (1.0-ss)*(1.0-tt) deriv(1, 9) = -s2*(1.0-tt) deriv(2, 9) = -t2*(1.0-ss) Return End Subroutine sfr2 !***************************************************************** !***************************************************************** !c jacob2 !***************************************************************** Subroutine jacob2(cartd, deriv, djacb, elcod, gpcod, ielem, kgasp, & nnode, shape) !C****************************************************************************** !C !C**** THIS SUBROUTINE EVALUATES THE JACOBIAN MATRIX AND THE CARTESIAN !C SHAPE FUNCTION DERIVATIVES !C !C****************************************************************************** Dimension cartd(2, 9), deriv(2, 9), elcod(2, 9), gpcod(2, 9), shape(9),& xjaci(2, 2), xjacm(2, 2) !c !c*****CALCULATE COORDINATES OF SAMPLING POINT !c Do idime = 1, 2 gpcod(idime, kgasp) = 0.0 Do inode = 1, nnode gpcod(idime, kgasp) = gpcod(idime, kgasp) + elcod(idime, inode)& *shape(inode) End Do End Do !C !C*** CREATE JACOBIAN MATRIX XJACM !C Do idime = 1, 2 Do jdime = 1, 2 xjacm(idime, jdime) = 0.0 Do inode = 1, nnode xjacm(idime, jdime) = xjacm(idime, jdime) + deriv(idime, inode)*& elcod(jdime, inode) End Do End Do End Do !C !C**** CALCULATE DETERMINANT AND INVERSE OF JACOBIAN MATRIX !C djacb = xjacm(1, 1)*xjacm(2, 2) - xjacm(1, 2)*xjacm(2, 1) If (djacb) 6, 6, 8 6 Write (6, 600) ielem Stop 8 Continue xjaci(1, 1) = xjacm(2, 2)/djacb xjaci(2, 2) = xjacm(1, 1)/djacb xjaci(1, 2) = -xjacm(1, 2)/djacb xjaci(2, 1) = -xjacm(2, 1)/djacb !c !c**** CALCULATE CARTESIAN DERIVATIVES !c Do idime = 1, 2 Do inode = 1, nnode cartd(idime, inode) = 0.0 Do jdime = 1, 2 cartd(idime, inode) = cartd(idime, inode) + xjaci(idime, jdime)*& deriv(jdihe, inode) !c !c PRELIMINARY THEORY AND STANDARD SUBROUTINES !c End Do End Do End Do Return 600 Format (//, ' PROGRAM HALTED IN SUBROUTINE JACOB2', /,11X, & 'ZERO OR NEGATIVE AREA',/,10X, 'ELEMENT NUMBER',I5) End Subroutine jacob2 !******************************************************************** !******************************************************************** ! modpb !******************************************************************** Subroutine modpb(dflex, dplan, dsher, lprop, mmats, props,& ifpla, iffle, ifshe) !c************************************************************* !c !c****CALCULATES MATRIX OF ELASTIC RIGIDITIES !c****FOR MINDLIN PLATE !c !c************************************************************* Dimension dflex(3, 3), dplan(3, 3), dsher(2, 2),& props(mmats, 8) young = props(lprop, 1) poiss = props(lprop, 2) thick = props(lprop, 3) !C*** FORM DPLAN If (ifpla==0) Goto 10 Do irows = 1, 3 Do jcols = 1, 3 dplan(irows, jcols) = 0.0 End Do End Do const = (young*thick)/(1.0-poiss*poiss) dplan(1, 1) = const dplan(2, 2) = const dplan(1, 2) = const*poiss dplan(2, 1) = const*poiss dplan(3, 3) = const*(1.0-poiss)/2.0 !c**** FORM DFLEX 10 If (ifflex==0) Goto 20 Do irows = 1, 3 Do jcols = 1, 3 dflex(irows, jcols) = 0.0 End Do End Do const = (young*thick**3)/(12.*(1.0-poiss*poiss)) dflex(1, 1) = const dflex(2, 2) = const dflex(1, 2) = const*poiss dflex(2, 1) = const*poiss dflex(3, 3) = const*(1.-poiss)/2 !C**** FORM DSHER 20 If (ifshe==0) Return Do irows = 1, 2 Do jcols = 1, 2 dsher(irows, jcols) = 0.0 End Do End Do dsher(1, 1) = (young*thick)/(2.4+2.4*poiss) dsher(2, 2) = (young*thick)/(2.4+2.4*poiss) Return End Subroutine modpb !********************************************************************** !********************************************************************** ! nodexy !*********************************************************************** Subroutine nodexy(coord, lnods, melem, mpoin, nelem, nnode) !c********************************************************************** !c !c****This subroutine interpolates the mid side nodes of straight !c sides of elements and the central node of 9 noded elements. !c !c********************************************************************** Dimension coord(mpoin, 2), lnods(melem, 9) If (nnode==4) Return !c !c***Loop over each element. !c Do ielem = 1, nelem !c !c***Loop over each element edge. ! c nnod1 = 9 If (nnode==8) nnod1 = 7 Do inode = 1, nnod1, 2 If (inode==9) Goto 50 ! c ! c***Compute the node number of the first node. ! c nodst = lnods(ielem, inode) igash = inode + 2 If (igash>8) igash = 1 ! c ! c***Compute the node number of the last node. ! c nodfn = lnods(ielem, igash) midpt = inode + 1 ! c ! c***Compute the node number of the intermediate node. ! c nodmd = lnods(ielem, midpt) total = abs(coord(nodmd,1)) + abs(coord(nodmd,2)) ! c ! c***If the coordinates of the intermediate node are both ! c zero interpolate by a straight line. ! c If (total>0.0) Goto 20 kount = 1 10 coord(nodmd, kount) = (coord(nodst,kount)+coord(nodfn,kount))/2.0 kount = kount + 1 If (kount==2) Goto 10 20 End Do Goto 30 50 lnode = lnods(ielem, inode) total = abs(coord(lnode,1)) + abs(coord(lnode,2)) If (total>0.0) Goto 30 lnod1 = lnods(ielem, 1) lnod3 = lnods(ielem, 3) lnod5 = lnods(ielem, 5) lnod7 = lnods(ielem, 7) kount = 1 40 coord(lnode, kount) = (coord(lnod1,kount)+coord(lnod3,kount)+& coord(lnod5,kount)+coord(lnod7,kount))/4.0 kount = kount + 1 If (kount==2) Goto 40 30 End Do Return End Subroutine nodexy !************************************************************ !************************************************************ ! 20230423 !************************************************************ ! flowmp !*********************************************************** Subroutine flowmp(abeta, avect, devia, dmatx, dvect, hards, ncrit, sint3, steff, theta, varj2) !c************************************************************ !c !C*** DETERMINES YIELD FUNCTION DERIVATIVES FOR MINDLIN PLATES !c*** 1 Von Mises !c*** 2 Tresca !c !c ELASTO-PLASTIC MINDLIN PLATE BENDING ANALYSIS !C !C************************************************************* Dimension avect(5), devia(4), dmatx(3, 3), dvect(5), veca1(3), veca2(3), veca3(3) !c !C*** DETERMINE THE VECTOR DERIVATIVE OF F FOR VON-MISES !c sinth = sin(theta) costh = cos(theta) root3 = 1.73205080757 !c !C*** CALCULATE VECTOR A1 !C veca1(1) = 0.333333333333 veca1(2) = 0.333333333333 vecal(3) = 0.0 !C !C*** CALCULATE VECTOR A2 !C Do istre = 1, 3 veca2(istre) = devia(istre)/(2.0*steff) End Do veca2(3) = devia(3)/steff !c !C*** CALCULATE VECTOR A3 !c veca3(1) = devia(2)*devia(4+varj2/3.0) veca3(2) = devia(1)*devia(4+varj2/3.0) veca3(3) = -2.0*devia(3)*devia(4) Goto (1, 2) ncrit !C !C*** VON MISES !C 1 cons1 = 0.0 cons2 = root3 cons3 = 0.0 Goto 40 !C !C*** TRESCA !C 2 cons1 = 0.0 abthe = abs(theta*57.29577951308) If (abthe

有限元隐含现象混合瞬态动力学分析 !********************************************************************** !********************************************************************** !********************************************************************** !C IMPLICIT-EXPLICIT TRANSIENT DYNAMIC ANALYSIS !C Overall structure of program MIXDYN !C Master routine MIXDYN !C PROGRAM MIXDYN (INPUT,TAPE5=INPUT,TAPE4,TAPElO,TAPEl2,TAPE3, !C OUTPUT,TAPE6=OUTPUT,TAPE7,TAPE11,TAPE13 !********************************************************************* !********************************************************************* ! ! MIXDYN IMPLICIT-EXPLICIT TRANSIENT DYNAMIC ANALYSIS ! 20230414 !********************************************************************* !********************************************************************* ! SINOMACH ! ! FORTRAN95,2013&2018 version ! ! ! Lean, Dobbsen,Beston, Tim, Andy, Young, Ted, Libberman SwanseaUK !********************************************************************* Program main !C******************************************************************** !C !C TIME INTEGRATION IMPLICIT-EXPLICIT ALGORITHM !C !C******************************************************************** Dimension coord(200,2),stiff(6000),dispi(400),posgp( 4),& ifpre(2,200),stifs(6000),veloi(400),weigp( 4),& lnods(50, 9),stifi(6000),accei(400),nprqd( 10),& rload(50,18),xmass(6000),displ(400),ngrqs( 10),& props(10,13),dampi(6000),velol(400),intgr( 50),& leqns(18,50),dampg(6000),accel(400),matno( 50),& strin(4,450),ymass( 400),accek(400),maxai(400),& strag(4,450),force( 400),accej(400),maxaj(400),& strsg(4,450),resid( 400),dispt(400),acceh(400),& epstn( 450),tempe( 400),dispq(400),accev(400),& niter( 2000),mhigh( 400),effst(450),velot(400) !C Common stiff, xmass, dampg, stifi, dampi !c !c******************************************************** !c**** call timestamp ( ) !c******************************************************** Open (5, File='data\input.dat', Status='unknown') Open (6, File='data\output.dat', Status='unknown') !c******************************************************** !C Call contol(ndofn, nelem, nmats, npoin) !C Call inputd(coord, ifpre, lnods, matno, nconm, ncrit, & ndime, ndofn, nelem, ngaum, ngaus, nlaps, & nmats, nnode, npoin, nprev, nstre, ntype, & posgp, props, weigp) !C Call intime(aalfa, acceh, accev, afact, azero, beeta, & bzero, delta, dtime, dtend, gaama, ifixd, & ifunc, intgr, kstep, miter, ndofn, nelem, & ngrqs, noutd, noutp, npoin, nprqd, nreqd, & nreqs, nstep, omega, dispi, toler, veloi, & ipred) ! IMPLICIT-EXPLICIT TRANSIENT DYNAMIC ANALYSIS Call prevos(force, ndofn, nelem, ngaus, npoin, nprev, & strin) Call loadpl(coord, force, lnods, matno, ndime, ndofn, & nelem, ngaus, nmats, nnode, npoin, nstre, & ntype, posgp, props, rload, strin, tempe, & weigp) !c Call lumass(coord, intgr, lnods, matno, nconm, ndime, & ndofn, nelem, ngaum, nmats, nnode, npoin, & ntype, props, ymass) Call linkin(force, ifpre, intgr, leqns, lnods, maxai, & maxaj, mhigh, ndofn, nelem, neqns, nnode, & npoin, nwktl, nwmtl, xmass, ymass) !c Do istep = 1, nstep !c Do iiter = 1, miter !c Call gstiff(coord, epstn, intgr, istep, kstep, leqns, & lnods, matno, maxai, maxaj, ncrit, ndime, & ndofn, nelem, ngaus, nlaps, nmats, nnode, & npoin, nstre, ntype, nwmtl, nwktl, posgp, & props, stiff, stifi, strsg, dispt, weigp) !c Call impexp(aalfa,acceh,accei,accej,accek,accel,& accev,afact,azero,beeta,bzero,consd,& consf,dampi,dampg,delta,dispi,displ,& dispt,dtend,dtime,gaama,ifixd,ifpre,& ifunc,iiter,istep,kstep,maxai,maxaj,& ndofn,neqns,npoin,nwktl,nwmtl,omega,& force,stiff,stifi,stifs,veloi,velol,& velot,xmass,ymass,ipred) !c Call resepl(coord,dispt,effst,rload,epstn,iiter,& intgr,leqns,lnods,matno,ncrit,ndime,& ndofn,nelem,ngaus,nlaps,nmats,nnode,& npoin,nstre,ntype,posgp,props,resid,& strag,strin,strsg,weigp,ipred,istep) !c Call itrate(accei,accel,consd,consf,xmass,dispi,& displ,dispt,maxai,nchek,neqns,nwmtl,& resid,stifs,toler,veloi,velol,velot,& iiter,miter) !c If(nchek==1) Goto 510 End Do !C 510 Call outdyn(dispq,dtime,epstn,ifpre,iiter,istep,& ndofn,nelem,ngaus,ngrqs,niter,noutd,& noutp,npoin,nprqd,nreqd,nreqs,ntype,& strsg,dispi) !c End Do !C !c******************************************************** Close (5) Close (6) !c******************************************************** !c**** call timestamp ( ) !c******************************************************** Stop End Program main !****************************************************************** !****************************************************************** ! blarge !****************************************************************** Subroutine blarge(bmatx, cartd, djacm, dlcod, gpcod, kgasp,& nlaps, nnode, ntype, shape) !c******************************************************* !c !c*** Large displacement B matrix !c !c******************************************************** Dimension bmatx(4, 18), cartd(2, 9), djacm(2, 2), dlcod(2, 9),& gpcod(2, 9), shape(9) ngash = 0 Do inode = 1, nnode mgash = ngash + 1 ngash = mgash + 1 bmatx(1, mgash) = cartd(1, inode)*djacm(1, 1) bmatx(1, ngash) = cartd(1, inode)*djacm(2, 1) bmatx(2, mgash) = cartd(2, inode)*djacm(1, 2) bmatx(2, ngash) = cartd(2, inode)*djacm(2, 2) bmatx(3, mgash) = cartd(2, inode)*djacm(1, 1) + cartd(1, inode)*& djacm(1, 2) bmatx(3, ngash) = cartd(1, inode)*djacm(2, 2) + cartd(2, inode)*& djacm(2, 1) End Do If (ntype/=3) Return fmult = 1 If (nlaps 50) Goto 200 If (npoin >200) Goto 200 !******************************************************* !******************************************************* If (nmats > 10) Goto 200 Goto 210 200 Write (6, 120) Stop Return 120 Format (/'SET DIMENSION EXCEEDED - CONTROL CHECK'/) 110 Format (16I5) 210 Continue End Subroutine contol !***************************************** !***************************************** !c inputd !***************************************** Subroutine inputd(coord, ifpre, lnods, matno, nconm, ncrit,& ndime, ndofn, nelem, ngaum, ngaus, nlaps,& nmats, nnode, npoin, nprev, nstre, ntype,& posgp, props, weigp) !C******************************************************* !C !C**** DYNPAK INPUT ROUTINE !C !C******************************************************* Dimension coord(npoin, 1),ifpre(ndofn, 1), weigp(1),& matno(1),lnods(nelem, 1), props(nmats, 1),& posgp(1),title(10) Read (5, 913) Write (6, 914) !C******************************************************* !C !C****READ THE FIRST DATA CARD,AND ECHIO IT IMMEDIATELY !C !C******************************************************* Read (5, 900) nvfix, ntype, nnode, nprop, ngaus, ndime, nstre, ncrit,& nprev, nconm, nlaps, ngaum, nrads Write (5, 901) npoin, nelem, nvfix, ntype, nnode, ndofn, nmats, nprop,& ngaus, ndime, nstre, ncrit, nprev, nconm, nlaps, ngaum,& nrads !C !C***** READ THE ELEMENT NODAL CONNECTIONS,AND THE PROPERTY NUMBERS !C Write (6, 902) Do ielem = 1, nelem Read (5, 900) numel, matno(numel),(lnods(numel,inode),inode=1, nnode) Write (13, 915) numel, (lnods(numel,inode), inode=1, nnode) Write (6, 903) numel, matno(numel), (lnods(numel,inode), inode=1, nnode) End Do !C !C*** ZERO ALL THE NODAL COORDINATES,PRIOR TO READING SOME OF THEM !C Do ipoin = 1, npoin Do idime = 1, ndime coord(ipoin, idime) = 0.0 End Do End Do !c !c ***** READ SOHE NODAL COORDINATES, FINISHING WITH THE LAST NODE OF ALL !c 200 Read (5, 905) ipoin, (coord(ipoin,idime), idime=1, ndime) Write (6, 906) ipoin, (coord(ipoin,idime), idime=1, ndime) If (ipoin/=npoin) Goto 200 !C !C*** INTERPOLATE COORDINATES OF MID-SIDE NODES !C !*************************************************************************** call nodxyr(coord, lnods, nelem, nnode, npoin, nrads, ntype) !C Write (6, 904) Write (13, 916)(ipoin, (coord(ipoin,idime),idime=1,ndime), ipoin=1, npoin) Write (13, 906)(ipoin, (coord(ipoin,idime),idime=1,ndime), ipoin=1, npoin) !C !C***READ THE FIXED VALUES !C Write (6, 907) Do ipoin = 1, npoin Do idofn = 1, ndofn ifpre(idofn, ipoin) = 0 End Do End Do Do ivfix = 1, nvfix Read (5, 908) ipoin, (ifpre(idofn,ipoin), idofn=1, ndofn) End Do Do ipoin = 1, npoin Write (6, 909) (ipoin, ifpre(idofn,ipoin), idofn=1, ndofn) End Do !C !C***READ THE AVAILABLE SELECTION OF ELEMENT PROPERTIES !C Write (6, 910) Do imats = 1, nmats Read (5, 900) numat Read (5, 917)(props(numat,iprop), iprop=1, nprop) Write (6,911) numat Write (6, 912)(props(numat,iprop), iprop=1, nprop) End Do !C !C****SET UP GANSSIAN INTERATION CONSTANTS !C Call gaussq(ngaus, posgp, weigp) Return 913 Format (10A4) 914 Format (//,5X, 10A4) 901 Format (/5X, 'CONTROL PARAMETERS'/& /5X, 'NPOIN=', I10, 5X, 'NELEM=', I10, 5X, 'NVFIX=', I10/& /5X, 'NTYPE=', I10, 5X, 'NNODE=', I10, 5X, 'NDOFN=', I10/& /5X, 'NMATS=', I10, 5X, 'NPROP=', I10, 5X, 'NGAUS=', I10/& /5X, 'NDIME=', I10, 5X, 'NSTRE=', I10, 5X, 'NCRIT=', I10/& /5X, 'NPREV=', I10, 5X, 'NCONM=', I10, 5X, 'NLAPS=', I10/& /5X, 'NGAUM=', I10, 5X, 'NRADS=', I10/) 900 Format (16I5) 902 Format (//5X, 'ELEMENT', 3X, 'PROPERTY', 6X, 'NODE NUMBERS') 903 Format (6X, I5, I9, 6X, 10I5) 915 Format (16I5) !C !C****READ SOME NODAL COORDIATES,FINISHING WITH THE LAST !C NODE OF ALL !C 904 Format (//5X, 'NODE', 9X, 'X', 9X, 'Y', 5X) 905 Format (I5, 6F10.5) 916 Format (I5, 2G15.6) 906 Format (5X, I5, 2F10.3) 907 Format (//5X, 'NODE', 2X, 'CODE') 908 Format (1X, I4, 3X, 2I1) 909 Format (6X, I5, 3X, 2I1) 910 Format (//5X, 'MATHRIAL PROPERTIES') 911 Format (/5X, 'MATERIAL NO', I5) 912 Format (/5X, 'YOUNG MODULUS', G12.4/5X, 'POISSON RATIO',G12.4/& 5X, 'THCKNESS' , G12.4/5X, 'MASS DENSITY', G12.4/& 5X, 'ALPHA TEMPR' , G12.4/5X, 'REFERENCE FO', G12.4/& 5X, 'HARDENING PAR', G12.4/5x, 'FRICT ANGLE', G12.4/& 5X, 'FLUDITY PAR' , G12.4/5X, 'EXP DELTA', G12.4/& 5X, 'NFLOW CODE' , G12.4) 917 Format (8E10.4) End Subroutine inputd !*************************************************************** !*************************************************************** !c $$$$ intime !*************************************************************** Subroutine intime(aalfa, acceh, accev, afact, azero, beeta,& bzero, delta, dtime, dtend, gaama, ifixd,& ifunc, intgr, kstep, miter, ndofn, nelem,& ngrqs, noutd, noutp, npoin, nprqd, nreqd,& nreqs, nstep, omega, tdisp, toler, veloc,& ipred) !C******************************************************* !C !C*** INITIAI VALUES AND TIME INTEGRATION DATA !C !C******************************************************* Dimension tdisp(1), acceh(1), nprod(1), intgr(1),& veloc(1), accev(1), ngrqs(1) !C******************************************************* !C !C**** READ TIME STEPPING AND SELECTIVE OUTPUT PARAMETERS !C !C******************************************************* Read (5, 902) nstep, noutd, noutp, nreqd, nreqs, nacce, ifunc,& ifixd, miter, kstep, ipred Read (5, 190) dtime, dtend, dtrec, aalfa, beeta, delta, gaama,& azero, bzero, omega, toler Write (6, 950) nstep, noutd, noutp, nreqd, nrrqs, nacce, ifunc,& ifixd, miter, kstep, ipreo Write (6, 960) dtime, dtend, dtrec, aalfa, beeta, delta, gaama, & azero, bzero, omega, toler !C !C*** SELECTED NODES AND GAUSS POINTS FOR OUTPUT !C Read (5, 902)(nprqd(ireqd), ireqd=1, nreqd) Read (5, 902)(ngrqs(ireqs), ireqs=1, nreqs) Write (6, 909) Write (6, 910)(nprqd(ireqd), ireqd=1, nreqd) Write (6, 911)(ngrqs(ireqs), ireqs=1, nreqs) !C !C** READ THE INDICATOR FOR EXPLICIT OR IMPLICIT ELEMENT !C Read (5, 902)(intgr(ielem), ielem=1, nelem) Write (6, 930) Write (6, 902)(intgr(ielem), ielem=1, nelem) !C !C*** INITIAL DISPLACEMENTS !C jpoin = 0 Do ipoin = 1, npoin Do idofn = 1, ndofn jpoin = jpoin + 1 tdisp(jpoin) = 0. veloc(jpoin) = 0. End Do End Do Write (6, 903) 200 Read (5, 904) ngash, xgash, ygash nposn = (ngash-1)*ndofn+1 tdisp(nposn) = xgash nposn = nposn+1 tdisp(nposn) = ygash Write (6, 905) ngash, xgash, ygash If (ngash/=npoin) Goto 200 !C !C***** INITIAL VELOCITIES !C Write (6, 906) 210 Read (5, 904) ngash, xgash, ygash nposn = (ngash-1)*ndofn + 1 veloc(nposn) = xgash nposn = nposn + 1 veloc(nposn) = ygash Write (6, 905) ngash, xgash, ygash If (nfash/=npoin) Goto 210 If (ifunc/=0) Goto 250 !C !C****READ ACCELEROGRAM DATA,X-DIREC FROM TAPE 7,Y-DIREC FROM TAPE 12 !C afact = dtrec/dtime If (ifixd-1) 220, 230, 240 220 Read (7, 907)(acceh(i), i=1, nacce) Write (6, 912) dtrec Write (6, 907)(acceh(i), i=1, nacce) Read (12, 907)(accev(i), i=1, nacce) Write (6, 913) dtrec Write (6, 907)(accev(i), i=1, nacce) Goto 250 230 Read (12, 907)(accev(i), i=1, nacce) Write (6, 913) dtrec Write (6, 907)(accev(i), i=1, nacce) Goto 250 240 Read (7, 907)(acceh(i), i=1, nacce) Write (6, 912) Write (6, 907)(acceh(i), i=1, nacce) 250 Continue Return 950 Format (/5X, 'TIME STEPPING PARAMETERS'/& /5X, 'NSTEP=', I5, 12X, 'NOUTD=', I5, 12X, 'NOUTP=', I5,/& /5X, 'NREQD=', I5, 12X, 'NREQS=', I5, 12X, 'NACCE=', I5,/& /5X, 'IFUNC=', I5, 12X, 'IFIXD=', I5, 12X, 'MITER =', I5,/& /5X, 'KSTEP=', I5, 12X, 'IPRED=', I5) 960 Format (/5X, 'DTIME =', G12.4, 5X, 'DTEND=', G12.4, 5X, 'DTREC=', G12.4/& /5X, 'AALFA=', G12.4, 5X, 'BEETA=', G12.4, 5X, 'DELTA=', G12.4/& /5X, 'GAAMA=', G12.4, 5X, 'AZERO=', G12.4, 5X, 'BZERO=', G12.4/& /5X, 'OMEGA=', G12.4, 5X, 'TOLER=', G12.4) 909 Format (//5X, 'SELECTIVE OUTPUT REQUESTED FOR FOLLOWING') 910 Format (/,5X, 'NODES', 10I5) 911 Format (5X, 'G.P.', 10I5) 902 Format (16I5) 190 Format (8F10.4) 930 Format (/5X, 'TYPE OF ELEMENT,IMPLICIT=1,EXPLICIT=2'/) 904 Format (I5, 2F10.5) 903 Format (//5X, 'NODE', 2X, 'INITIAL X-DISP.', 2X,& 'INITIAL Y-DISP.'/) 905 Format (I10, 2E18.5) 906 Format (//5X, 'NODE', 2X, 'INITIAL X-VELO.', 2X,& 'INITIAL Y-VELO.'/) 907 Format (7F10.3) 912 Format (/5X, 'HORIZONTAL ACCELERATION ORDINATES AT', F9.4,2X,'SEC'/) 913 Format (/5X, 'VERTICAL ACCELERATION ORDINATES AT', F9.4,2X,'SEC'/) End Subroutine intime !********************************************************* !********************************************************* ! c$ LOADPL !********************************************************* Subroutine loadpl(coord, force, lnods, matno, ndime, ndofn,& nelem, ngaus, nmats, nnode, npoin, nstre, & ntxpe, posgp, props, rload, strin, tempe,& weigp) !C********************************************************* !C !C***STANDARD LOAD ROUTINE !C !C********************************************************* Dimension coord(npoin, 1), gpcod(2, 9), posgp(1), stran(4),& lnods(nelem, 1), cartd(2, 9), weigp(1), stres(4),& props(nmats, 1), deriv(2, 9), tempe(1), noprs(3),& rload(nelem, 1), elcod(2, 9), matno(1), dgash(2),& strin( 4, 1), press(3, 2), shape(9), pgash(2),& dmatx( 4, 4), title( 10), point(2), force(1) twopi = 6.283185307179586 nevab = nnode*ndofn Do ielem = 1, nelem Do ievab = 1, nevab rload(ielem, ievab) = 0.0 End Do End Do Read (5, 901) title Write (6, 903) title !C********************************************************** !C !C**** READ DATA CONTRLLING LOADING TYPES TO BE INPUTTED !C !C********************************************************** Read (5, 919) iplod, igrav, iedge, itemp Write (6, 990) Write (6, 991) iplod, igrav, iedge, itemp !C !C*** READ NODAL POINT LOADS !C If (iplod==0) Goto 500 Write (6, 998) 20 Read (5, 931) lodpt, (point(idofn), idofn=1, ndofn) Write (6, 933) lodpt, (point(idofn), idofn=1, ndofn) !C !C*** ASSOCIATE THE NODAL POINT LOADS WITH AN ELEMENT !C Do ielem = 1, nelem Do inode = 1, nnode nloca = iabs(lnods(ielem,inode)) If (lodpt==nloca) Goto 40 End Do End Do 40 Do idofn = 1, ndofn ngash = (inode-1)*ndofn + idofn rload(ielem, ngash) = point(idofn) End Do If (lodptncode) ngash = 1 If (knode>ncode) mgash = 2 rload(neass, ngash) = rload(neass, ngash) + shape(kount)*& pxcom*dvolu rload(neass, mgash) = rload(neass, mgash) + shape(kount)*& pycom*dvolu End Do End Do End Do 700 Continue If (itemp==0) Goto 800 !C !C*** INITIALIZE AND INPUT THE NODAL TEMPERATURES !C Do ipoin = 1, npoin tempe(ipoin) = 0.0 End Do Write (6, 917) 180 Read (5, 916) nodpt, tempe(nodpt) Write (6, 916) nodpt, tempe(nodpt) If (nodpt

!*********************************************************************** !*********************************************************************** ! FEM法一维非线性弹塑性材料力学问题续 !*********************************************************************** ! MASTER UNIDIM organisation for one-dimensional(1D) ! nonlinear applications. !*********************************************************************** ! SINOMACH ! ! FORTRAN95,2013&2018 version ! ! ! PotsyYZ,XumZ,Johnny, Norck, ,SwanseaUK ! ! 20230407 !*********************************************************************** !C********************************************************************** !C !C *** PROGRM FOR THE 1-D SOLUTION OF NONLINEAR PROBLEMS !C (Continued I) !C !C********************************************************************** ! c$debug ! c$ master unidim !c*********************************************************** ! Implicit none ! call timestamp ( ) !c*********************************************************** Program Main !c*********************************************************************** Common /unim1/npoin, nelem, nboun, nload, nprop, nnode, iincs, iiter,& kresl, nchek, toler, nalgo, nsvab, ndofn, nincs, nevab,& niter, noutp, facto, pvalu Common /unim2/props(5, 4), coord(26), lnods(25, 2), ifpre(52),& fixed(52), tload(25, 4), rload(25, 4), eload(25, 4),& matno(25), stres(25, 2), plast(25), xdisp(52),& tdisp(26, 2), treac(26, 2), astif(52, 52), aslod(52),& react(52), fresv(1352), pefix(52), estif(4, 4) !c******************************************************** ! Open (5, File='data\input.dat', Status='unknown') ! Open (6, File='data\output.dat', Status='unknown') !c******************************************************** Call data Call inital Do 30 iincs = 1, nincs Call inclod Do 10 iiter = 1, niter Call nonal If (kresl==1) Call stiff1 Call assemb if(kresl==1) call greduc if(kresl==2) call resolv Call baksub !***************************************** !c$$$ ???? Call refor1/2/3 ! 20230406 !***************************************** Call refor1 Call monitr If(nchek==0) Goto 20 If(iiter==1 .And. noutp==1) Call result If(noutp==2) Call result 10 Continue Write(6,900) 20 Call result ! call timestamp ( ) 30 Continue ! c****************** ! Close (5) ! Close (6) ! c****************** ! call timestamp ( ) Stop 900 Format('0',5X,'SOLUTION NOT CONVERGED') End Program !************************************************************ !C*********************************************************** !c data !c*********************************************************** Subroutine data !C*********************************************************** !C !C*****INPUTS DATA DEFINING GEOMETRY,LOADING,BOUNDARY CONDITIONS...ETC. !C !C*********************************************************** Common /unim1/npoin, nelem, nboun, nload, nprop, nnode, iincs, iiter,& kresl, nchek, toler, nalgo, nsvab, ndofn, nincs, nevab,& niter, noutp, facto, pvalu Common /unim2/props(5, 4), coord(26), lnods(25, 2), ifpre(52),& fixed(52), tload(25, 4), rload(25, 4), eload(25, 4),& matno(25), stres(25, 2), plast(25), xdisp(52),& tdisp(26, 2), treac(26, 2), astif(52, 52), aslod(52),& react(52), fresv(1352), pefix(52), estif(4, 4) Dimension icode(2), value(2), title(18) Read (5, 965) Write (6, 965) Read (5, 900) npoin, nelem, nboun, nmats, nprop, nnode, nincs, nalgo, ndofn Write (6, 905) npoin, nelem, nboun, nmats, nprop, nnode, nincs, nalgo, ndofn nevab = ndofn*nnode nsvab = ndofn*npoin Write (6, 910) Do imats = 1, nmats Read (5, 915) jmats, (props(jmats,iprop), iprop=1, nprop) Write (6, 915) jmats, (props(jmats,iprop), iprop=1, nprop) End Do Write (6, 920) Do ielem = 1, nelem Read (5, 925) jelem, (lnods(jelem,inode), inode=1, nnode), matno(jelem) Write (6, 925) jelem, (lnods(jelem,inode), inode=1, nnode), matno(jelem) End Do Write (6, 930) Do ipoin = 1, npoin Read (5, 935) jpoin, coord(jpoin) Write (6, 935) jpoin, coord(jpoin) End Do Do isvab = 1, nsvab ifpre(isvab) = 0 pefix(isvab) = 0.0 End Do If (ndofn==1) Write (6, 940) If (ndofn==2) Write (6, 945) Do iboun = 1, nboun Read (5, 950) nodfx, (icode(idofn), value(idofn), idofn=1, ndofn) Write (6, 950) nodfx, (icode(idofn), value(idofn), idofn=1, ndofn) nposn = (nodfx-1)*ndofn Do idofn = 1, ndofn nposn = nposn + 1 ifpre(nposn) = icode(idofn) pefix(nposn) = value(idofn) End Do End Do Write (6, 955) Do ielem = 1, nelem !C !C FINITE ELEMENTS IN PLASTICITY !C Do ievab = 1, nevab rload(ielem, ievab) = 0.0 End Do End Do 70 Read (5, 960) jelem, (rload(jelem,ievab), ievab=1, nevab) If (jelem/=nelem) Goto 70 Do ielem = 1, nelem Write (6, 960) ielem, (rload(ielem,ievab), ievab=1, nevab) End Do Return 965 Format (18A4) 900 Format (9I5) 905 Format (//1X, 'NPOIN =', I5, 3X, 'NELEM =', I5, 3X, 'NBOUN =', I5,& 3X, 'NMATS =', I5//1X, 'NPROP =', I5, 3X, 'NNODE =', I5,& 3X, 'NINCS =', I5, 3X, 'NALGO =', I5//1X, 'NDOFN =', I5) 910 Format ('0', 5X, 'MATERIAL PROPERTIES') 915 Format (I10, 4F15.5) 920 Format ('0', 3X, 'EL NODES MAT.') 925 Format (4I5) 930 Format ('0', 4X, 'NODE', 5X, 'COORD') 935 Format (I10, F15.5) 940 Format ('0', 1X, 'RES.NODE', 2X, 'CODE', 3X, 'PRES.VALUES') 945 Format ('0', 1X, 'RES.NODE', 2X, 'CODE', 3X, 'PRES.VALUES', 2X,& 'CODE', 3X, 'PRES VALUES') 950 Format (I10, 2(I5,F15.5)) 955 Format ('0', 2X, 'ELEMENT', 10X, 'NODAL LOADS') 960 Format (I10, 5F15.5) End Subroutine data !c******************************************************* !c******************************************************* !c inital !c*********************************************************** Subroutine inital !C*********************************************************** !C !C**** INITIALIZES TO ZERO ALL ACCUNULATIVE ARRAYS !C !C*********************************************************** Common /unim1/npoin, nelem, nboun, nload, nprop, nnode, iincs, iiter,& kresl, nchek, toler, nalgo, nsvab, ndofn, nincs, nevab,& niter, noutp, facto, pvalu Common /unim2/props(5, 4), coord(26), lnods(25, 2), ifpre(52),& fixed(52), tload(25, 4), rload(25, 4), eload(25, 4),& matno(25), stres(25, 2), plast(25), xdisp(52),& tdisp(26, 2), treac(26, 2), astif(52, 52), aslod(52),& react(52), fresv(1352), pefix(52), estif(4, 4) Do ielem = 1, nelem plast(ielem) = 0.0 Do idofn = 1, ndofn stres(ielem, idofn) = 0.0 End Do Do ievab = 1, nevab eload(ielem, ievab) = 0.0 tload(ielem, ievab) = 0.0 End Do End Do Do ipoin = 1, npoin Do idofn = 1, ndofn tdisp(ipoin, idofn) = 0.0 treac(ipoin, idofn) = 0.0 End Do End Do Return End Subroutine inital !c************************************************************ !c************************************************************ !c inclod !c************************************************************ Subroutine inclod !C************************************************************** !C !C *** INPUTS DATA FOR CURRENT INCREMENT AND UPDATES LOAD VECTOR !C !C************************************************************** Common /unim1/npoin, nelem, nboun, nload, nprop, nnode, iincs, iiter,& kresl, nchek, toler, nalgo, nsvab, ndofn, nincs, nevab,& niter, noutp, facto, pvalu Common /unim2/props(5, 4), coord(26), lnods(25, 2), ifpre(52),& fixed(52), tload(25, 4), rload(25, 4), eload(25, 4),& matno(25), stres(25, 2), plast(25), xdisp(52),& tdisp(26, 2), treac(26, 2), astif(52, 52), aslod(52),& react(52), fresv(1352), pefix(52), estif(4, 4) Read (5, 900) niter, noutp, facto, toler Write (6, 905) iincs, niter, noutp, facto, toler Do ielem = 1, nelem Do ievab = 1, nevab eload(ielem, ievab) = eload(ielem, ievab) + rload(ielem, ievab)*facto tload(ielem, ievab) = tload(ielem, ievab) + rload(ielem, ievab)*facto End Do End Do Return 900 Format (2I5, 2F15.5) 905 Format ('0', 5X, 'IINCS =', I5, 3X, 'NITER =', I5, 3X, 'NOUTP=', I5,& 3X, 'FACTO =', E14.6, 3X, 'TOLER =', E14.6) End Subroutine inclod !C*********************************************************** !c************************************************************* !c nonal !C*********************************************************** Subroutine nonal !********************************************************** ! !******SETS INDICATOR TO DENTIFY TYPE OF SOLUTION ALGORITHM ! !********************************************************** Common /unim1/npoin, nelem, nboun, nload, nprop, nnode, iincs, iiter,& kresl, nchek, toler, nalgo, nsvab, ndofn, nincs, nevab,& niter, noutp, facto, pvalu Common /unim2/props(5, 4), coord(26), lnods(25, 2), ifpre(52),& fixed(52), tload(25, 4), rload(25, 4), eload(25, 4),& matno(25), stres(25, 2), plast(25), xdisp(52),& tdisp(26, 2), treac(26, 2), astif(52, 52), aslod(52),& react(52), fresv(1352), pefix(52), estif(4, 4) kresl = 2 If (nalgo==1) kresl = 1 If (nalgo==2) kresl = 1 If (nalgo==3 .And. iincs==1 .And. iiter==l) kresl = 1 If (nalgo==4 .And. iiter==1) kresl = 1 If (nalgo==5 .And. iincs==1 .And. iiter==l) kresl = 1 If (nalgo==5 .And. iiter==2) kresl = 1 If (iiter==1 .Or. nalgo==1) Goto 20 Do isvab = 1, nsvab fixed(isvab) = 0.0 End Do Return 20 Do isvab=1, nsvab fixed(isvab) = pefix(isvab)*facto End Do Return End Subroutine nonal !C*********************************************************** !c*********************************************************** ! c$debug ! c$large !********************************************************** Subroutine stiff1 ! c******************************************************************** ! c ! c ******Calculates element stiffness matrices. ! c ! c******************************************************************** Common /unim1/npoin, nelem, nboun, nload, nprop, nnode, iincs, iiter, & kresl, nchek, toler, nalgo, nsvab, ndofn, nincs, nevab, & niter, noutp, facto, pvalu Common /unim2/props(5, 4), coord(26), lnods(25, 2), ifpre(52),& fixed(52), tload(25, 4), rload(25, 4), eload(25, 4),& matno(25), stres(25, 2), plast(25), xdisp(52), & tdisp(26, 2), treac(26, 2), astif(52, 52), aslod(52),& react(52), fresv(1352), pefix(52), estif(4, 4) !******************************************* !c$ ? Rewind 1 ?? !******************************************* Do ielem = 1, nelem lprop = matno(ielem) sterm = props(lprop, 1) node1 = lnods(ielem, 1) node2 = lnods(ielem, 2) eleng = abs(coord(node1)-coord(node2)) averg = (tdisp(node1,1)+tdisp(node2,1))/2 !******************************************* !c$ ? fmult = sterm*varia/eleng ?? ! 20230407 !******************************************* fmult = sterm*1.25/eleng estif(1, 1) = fmult estif(1, 2) = -fmult estif(2, 1) = -fmult estif(2, 2) = fmult Write (1) estif End Do Return End Subroutine stiff1 !********************************************************** !c*********************************************************** !c assemb !C*********************************************************** Subroutine assemb !c*********************************************************** !c !c ***** ELEMENT ASSEMBLY ROUTINE !c !c************************************************************ Common /unim1/npoin, nelem, nboun, nload, nprop, nnode, iincs, iiter,& kresl, nchek, toler, nalgo, nsvab, ndofn, nincs, nevab,& niter, noutp, facto, pvalu Common /unim2/props(5, 4), coord(26), lnods(25, 2), ifpre(52),& fixed(52), tload(25, 4), rload(25, 4), eload(25, 4),& matno(25), stres(25, 2), plast(25), xdisp(52),& tdisp(26, 2), treac(26, 2), astif(52, 52), aslod(52),& react(52), fresv(1352), pefix(52), estif(4, 4) !c************************************************************ !c !c ELEMENT ASSEMBLY ROUTINE !c !c************************************************************ !c??? REWIND 1 ???? !c************************ Do isvab = 1, nsvab aslod(isvab) = 0.0 End Do If (kresl==2) Goto 30 Do isvab = 1, nsvab Do jsvab = 1, nsvab astif(isvab, jsvab) = 0.0 End Do End Do 30 Continue !c !c ASSEMBLE THE ELEMENT LOADS !c Do ielem = 1, nelem Read (1) estif Do inode = 1, nnode nodei = lnods(ielem, inode) Do idofn = 1, ndofn nrows = (nodei-1)*ndofn + idofn nrowe = (inode-1)*ndofn + idofn aslod(nr0ws) = aslod(nrows) + eload(ielem, nrowe) !c !c ASSEMBLE THE ELEMENT STIFFNESS MATRICES !c If (kresl==2) Goto 40 Do jnode = 1, nnode nodej = lnods(ielem, jnode) Do jdofn = 1, ndofn ncols = (nodej-1)*ndofn + jdofn ncole = (jnode-1)*ndofn + jdofn astif(nrows, ncols) = astif(nrows, ncols) + estif(nrowe, ncole) 40 End Do End Do End Do End Do End Do Return End Subroutine assemb !c*********************************************************** !C*********************************************************** !c &&& greduc !C*********************************************************** Subroutine greduc !C******************************************************************* !C !C **** GAUSSIAN REDUCTION ROUTINE !C !C******************************************************************* Common /unim1/npoin, nelem, nboun, nload, nprop, nnode, iincs, iiter,& kresl, nchek, toler, nalgo, nsvab, ndofn, nincs, nevab,& niter, noutp, facto, pvalu Common /unim2/props(5, 4), coord(26), lnods(25, 2), ifpre(52),& fixed(52), tload(25, 4), rload(25, 4), eload(25, 4),& matno(25), stres(25, 2), plast(25), xdisp(52),& tdisp(26, 2), treac(26, 2), astif(52, 52), aslod(52),& react(52), fresv(1352), pefix(52), estif(4, 4) !C !C GAUSSIAN REDUCTION ROUTINE !C kount = 0 neqns = nsvab Do ieqns = 1, neqns If (ifpre(ieqns)==1) Goto 40 !C !C REDUCE EQUATIONS !C pivot = astif(ieqns, ieqns) If (abs(pivot)toler) nchek = 1 20 pvalu = rcurr Write (6, 900) nchek, ratio Return 900 Format ('0', 5X, 'CONVERGENCE CODE =', I4, 3X, 'NORM OF RESIDUAL SUM& RATIO =', E14.6) End Subroutine monitr !*************************************************************** !************************************************************* ! result !************************************************************* Subroutine result !C*********************************************************** !C !C*****OUTPUTS DISPLACEMENT, REACTIONS AND STRESSES !C !C*********************************************************** Common /unim1/npoin, nelem, nboun, nload, nprop, nnode, iincs, iiter,& kresl, nchek, toler, nalgo, nsvab, ndofn, nincs, nevab,& niter, noutp, facto, pvalu Common /unim2/props(5, 4), coord(26), lnods(25, 2), ifpre(52),& fixed(52), tload(25, 4), rload(25, 4), eload(25, 4),& matno(25), stres(25, 2), plast(25), xdisp(52),& tdisp(26, 2), treac(26, 2), astif(52, 52), aslod(52),& react(52), fresv(1352), pefix(52), estif(4, 4) If(ndofn==1) Write(6,900) If(ndofn==2) Write(6,910) Do ipoin=1, npoin Write(6,920) ipoin, (tdisp(ipoin,idofn),treac(ipoin,idofn),idofn=1,ndofn) End Do If(ndofn==2) Write(6,930) If(ndofn==1) Write(6,940) Do ielem=1, nelem Write(6,950) ielem,(stres(ielem,idofn),idofn=1,ndofn), & plast(ielem) End Do Return 900 Format('0',5X,'NODE',4X,'DISPL.',12X,'REACTIONS') 910 Format('0',5X,'NODE',4X,'DISPL.',12X,'REACTION',& 7X,'DISPL.',12X,'REACTION') 920 Format(I10,2(E14.6,5X,E14.6)) 930 Format('0',2X,'ELE1ENT',12X,'STRESSES',12X,'PL.STRAIN') 940 Format('0',2X,'ELEMENT',5X,'STRESSES',5X,'PL.STRAIN') 950 Format(I10,3E14.6) End Subroutine result !**************************** !**************************** 2 1 3 9 14 15 16 10 5 4 3 1 5 10 16 17 18 11 7 6 4 1 12 19 23 24 25 20 14 13 5 1 14 20 25 26 27 21 16 15 6 1 16 21 27 28 29 22 18 17 7 1 23 30 34 35 36 31 25 24 8 1 25 31 36 37 38 32 27 26 9 1 27 32 38 39 40 33 29 28 10 1 34 41 45 46 47 42 36 35 11 1 36 42 47 48 49 43 38 37 12 1 38 43 49 50 51 44 40 39 !********************************** ! Output: !********************************** NPOIN = 0 NELEM = 3 NBOUN = 1 NMATS = 0 NPROP = 0 NNODE = 51 NINCS = 1 NALGO = 61 NDOFN = 71 0 MATERIAL PROPERTIES 0 EL NODES MAT. !***************************** !**************************** ! 20230407 !***************************** ! 广州 !

!********************************************************************* !********************************************************************* !********************************************************************* ! Modified and new routines ! Master BEML This routine is almost identical to routine BEAM dcscribed ! earlier. ! New routines for nonlayered elasto-plastic Timoshenko beam analysis ! Master BEAM The master calling routine BEAM simply organises the ! calling of the main routines as in one dimension(1D). ! !********************************************************************* ! SINOMACH ! ! FORTRAN95,2013&2018 version ! ! ! PotsyYZ, 周勇，Johnny, Frank, Mohammed,SwanseaUK ! ! 20230324 !********************************************************************* !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& !********************************************************************* ! c$debug ! c$BEML ! Program master beml TIMLAY !********************************************************************* ! Implicit none ! call timestamp ( ) !********************************************************************* Program main !c********************************************************** !c !c****Elasto-Plastic Layered Timoshenko Beam Program !c !c********************************************************** Common /unim1/npoin, nelem, nboun, nlayr, nprop, nnode, iincs, iiter,& kresl, nchek, toler, nalgo, nsvab, ndofn, nincs, nevab,& niter, noutp, facto Common /unim2/props(5, 25), coord(26), lnods(25, 2), ifpre(52),& fixed(52), tload(25, 4), rload(25, 4), eload(25, 4), matno(25),& stres(25, 2), plast(250), xdisp(52), tdisp(26, 2), treac(26, 2),& astif(52, 52), aslod(52), react(52), fresv(1352),& pefix(52), estif(4, 4), strsl(250, 2) Call data Call inital Do iincs = 1, nincs Call inclod Do iiter = 1, niter Call nonal If (kresl==1) Call stifbl Call assemb If (kresl==1) Call greduc If (kresl==2) Call resolv Call baksub Call rforbl Call conund If (nchek==0) Goto 20 If (iiter==1 .And. noutp==1) Call result If (noutp==2) Call result End Do Write (6, 900) Stop 20 Call result End Do Stop 900 Format ('0', 5X, 'SOLUTION NOT CONVERGED') End Program main !***************************************************************** !***************************************************************** Subroutine stifbl !c**************************************************************** !c !c ***** CALCULATES ELEMENT STIFFNESS MATRICES !c !c**************************************************************** Common /unim1/npoin, nelem, nboun, nlayr, nprop, nnode, iincs,& iiter, kresl, nchek, toler, nalgo, nsvab, ndofn, nincs,& nevab, niter, noutp, facto Common /unim2/props(5, 25), coord(26), lnods(25, 2), ifpre(52),& fixed(52), tload(25, 4), rload(25, 4), eload(25, 4), matno(25), & stres(25, 2), plast(250), xdisp(52), tdisp(26, 2), treac(26, 2),& astif(52, 52), aslod(52), react(52), fresv(1352), pefix(52), & estif(4, 4), strsl(250, 2) !c ??? REWIND 1 ??? Do ielem = 1, nelem lprop = matno(ielem) Call layer(ielem, eival, svalu) hards = props(lprop, 4) node1 = lnods(ielem, 1) node2 = lnods(ielem, 2) eleng = abs(coord(node2)-coord(node1)) valu1 = 0.5*svalu valu2 = svalu/eleng valu3 = eival/eleng valu4 = 0.25*svalu*eleng estif(1, 1) = valu2 estif(1, 2) = valu1 estif(1, 3) = -valu2 estif(1, 4) = valu1 estif(2, 2) = valu3 + valu4 estif(2, 3) = -valu1 estif(2, 4) = -valu3 + valu4 estif(3, 3) = valu2 estif(3, 4) = -valu1 estif(4, 4) = valu3 + valu4 Do istif = 1, 4 Do jstif = istif, 4 estif(jstif, istif) = estif(istif, jstif) End Do End Do Write (1) estif End Do Return End Subroutine stifbl !******************************************** !******************************************** Subroutine rforbl !c*********************************************** !c !c*** CALCULATES INTERNAL EQUIVALENT NODAL FORCES !c !c************************************************ Common /unim1/npoin, nelem, nboun, nlayr, nprop, nnode,& iincs, iiter, kresl, nchek, toler, nalgo, nsvab, ndofn,& nincs, nevab, niter, noutp, facto Common /unim2/props(2, 5), coord(26), lnods(25, 2),& ifpre(52), fixed(52), tload(25, 4), rload(25, 4),& eload(25, 4), matno(25), stres(25, 2), plast(250), & xdisp(52), tdisp(26, 2), treac(26, 2), astif(52, 52), & aslod(52), react(52), fresv(1352), pefix(52),& estif(4, 4), strsl(250, 2) Dimension stran(2) Do ielem = 1, nelem Do ievab = 1, nevab eload(ielem, ievab) = 0.0 End Do Do idofn = l, ndofn stres(ielem, idofn) = 0.0 End Do End Do klay r = 0 Do ielem = i, nelem lprop = matno(ielem) young = props(lprop, 1) shear = props(lprop, 2) yield = props(lprop, 3) hards = props(lprop, 4) thkto = props(lprop, 5) node1 = lnods(ielem, 1) node2 = lnods(ielem, 2) eleng = abs(coord(node2)-coord(node1)) wnod1 = xdisp(node1*ndofn-1) wnod2 = xdisp(node2*ndofn-1) thta1 = xdisp(node1*ndofn) thta2 = xdisp(node2*ndofn) stran(1) = (thta1-thta2)/eleng stran(2) = (wnod2-wnod1)/eleng - 0.5*(thta1+thta2) zmidl = -thkto/2.0 kount = 5 Do ilayr = 1, nlayr klayr = klayr + 1 kount = kount + 1 brdth = props(lprop, kount) kount = kount + 1 thick = props(lprop, kount) zmidl = zmidl + thick/2.0 stlin = young*stran(1)*zmidl stcur = strsl(klayr, 1) + stlin preys = yield + hards*abs(plast(klay)) If (abs(strsl(klayr,1))>=preys) Goto 20 escur = abs(stcur) - preys If (escur0.0 .And. stlinpvalu) nchek = 999 50 pvalu = ratio Write (6, 900) iiter, nchek, ratio Return 900 Format ('0', 5X, 'ITERATION NUMBER =', I5/'0',& 5X, 'CONVERGENCE CODE =', I4, 3X, 'NORM OF RESIDUAL SUM RATIO =', E14.6) End Program !C*********************************************************** !C*********************************************************** !c result !c*********************************************************** Subroutine result !C*********************************************************** !C !C*****OUTPUTS DISPLACEMENT, REACTIONS AND STRESSES !C !C*********************************************************** Common /unim1/npoin, nelem, nboun, nload, nprop, nnode, iincs, iiter,& kresl, nchek, toler, nalgo, nsvab, ndofn, nincs, nevab,& niter, noutp, facto, pvalu Common /unim2/props(5, 4), coord(26), lnods(25, 2), ifpre(52), & fixed(52), tload(25, 4), rload(25, 4), eload(25, 4),& matno(25), stres(25, 2), plast(25), xdisp(52), & tdisp(26, 5), treac(26, 2), astif(50), aslod(52), & react(52), fresv(1352), pefix(52), estif(4, 4) If (ndofn==1) Write (6, 900) If (ndofn==2) Write (6, 910) Do ipoin = 1, npoin Write (6, 920) ipoin, (tdisp(ipoin,idofn), treac(ipoin,idofn), idofn=1, ndofn) End Do If (ndofn==2) Write (6, 930) If (ndofn==1) Write (6.940) Do ielem = i, nelem Write (6, 950) ielem(stres(ielem,idofn), idofn=1, ndofn), plast(ielem) End Do Return 900 Format ('0', 5X, 'NODE', 4X, 'DISPL.', 12X, 'REACTIONS') 910 Format (Lh0, 5X, 'NODE', 4X, 'DISPL.', 12X, 'REACTION', 7X,& 'DISPL.', 12X, 'REACTION') 920 Format (I10, 2(E14.2,5X,E14.6)) 930 Format ('0', 2X, 'ELEMENT', 12X, 'STRESSES', 12X, 'PL.STRAIN') 940 Format ('0', 2X, 'ELEMENT', 5X, 'STRESSES', 5X, 'PL.STRAIN') End Subroutine result

!********************************************************************* !********************************************************************* !********************************************************************* ! Modified and new routines ! Master BEML This routine is almost identical to routine BEAM dcscribed ! earlier. ! New routines for nonlayered elasto-plastic Timoshenko beam analysis ! Master BEAM The master calling routine BEAM simply organises the ! calling of the main routines as in one dimension(1D). ! !********************************************************************* ! SINOMACH ! ! FORTRAN95,2013&2018 version ! ! ! PotsyYZ, 周勇，Johnny, Frank, Mohammed,SwanseaUK ! ! 20230324 !********************************************************************* !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& !********************************************************************* ! c$debug ! c$BEML ! Program master beml TIMLAY !********************************************************************* ! Implicit none ! call timestamp ( ) !********************************************************************* Program main !c********************************************************** !c !c****Elasto-Plastic Layered Timoshenko Beam Program !c !c********************************************************** Common /unim1/npoin, nelem, nboun, nlayr, nprop, nnode, iincs, iiter,& kresl, nchek, toler, nalgo, nsvab, ndofn, nincs, nevab,& niter, noutp, facto Common /unim2/props(5, 25), coord(26), lnods(25, 2), ifpre(52),& fixed(52), tload(25, 4), rload(25, 4), eload(25, 4), matno(25),& stres(25, 2), plast(250), xdisp(52), tdisp(26, 2), treac(26, 2),& astif(52, 52), aslod(52), react(52), fresv(1352),& pefix(52), estif(4, 4), strsl(250, 2) Call data Call inital Do iincs = 1, nincs Call inclod Do iiter = 1, niter Call nonal If (kresl==1) Call stifbl Call assemb If (kresl==1) Call greduc If (kresl==2) Call resolv Call baksub Call rforbl Call conund If (nchek==0) Goto 20 If (iiter==1 .And. noutp==1) Call result If (noutp==2) Call result End Do Write (6, 900) Stop 20 Call result End Do Stop 900 Format ('0', 5X, 'SOLUTION NOT CONVERGED') End Program main !***************************************************************** !***************************************************************** Subroutine stifbl !c**************************************************************** !c !c ***** CALCULATES ELEMENT STIFFNESS MATRICES !c !c**************************************************************** Common /unim1/npoin, nelem, nboun, nlayr, nprop, nnode, iincs,& iiter, kresl, nchek, toler, nalgo, nsvab, ndofn, nincs,& nevab, niter, noutp, facto Common /unim2/props(5, 25), coord(26), lnods(25, 2), ifpre(52),& fixed(52), tload(25, 4), rload(25, 4), eload(25, 4), matno(25), & stres(25, 2), plast(250), xdisp(52), tdisp(26, 2), treac(26, 2),& astif(52, 52), aslod(52), react(52), fresv(1352), pefix(52), & estif(4, 4), strsl(250, 2) !c ??? REWIND 1 ??? Do ielem = 1, nelem lprop = matno(ielem) Call layer(ielem, eival, svalu) hards = props(lprop, 4) node1 = lnods(ielem, 1) node2 = lnods(ielem, 2) eleng = abs(coord(node2)-coord(node1)) valu1 = 0.5*svalu valu2 = svalu/eleng valu3 = eival/eleng valu4 = 0.25*svalu*eleng estif(1, 1) = valu2 estif(1, 2) = valu1 estif(1, 3) = -valu2 estif(1, 4) = valu1 estif(2, 2) = valu3 + valu4 estif(2, 3) = -valu1 estif(2, 4) = -valu3 + valu4 estif(3, 3) = valu2 estif(3, 4) = -valu1 estif(4, 4) = valu3 + valu4 Do istif = 1, 4 Do jstif = istif, 4 estif(jstif, istif) = estif(istif, jstif) End Do End Do Write (1) estif End Do Return End Subroutine stifbl !******************************************** !******************************************** Subroutine rforbl !c*********************************************** !c !c*** CALCULATES INTERNAL EQUIVALENT NODAL FORCES !c !c************************************************ Common /unim1/npoin, nelem, nboun, nlayr, nprop, nnode,& iincs, iiter, kresl, nchek, toler, nalgo, nsvab, ndofn,& nincs, nevab, niter, noutp, facto Common /unim2/props(2, 5), coord(26), lnods(25, 2),& ifpre(52), fixed(52), tload(25, 4), rload(25, 4),& eload(25, 4), matno(25), stres(25, 2), plast(250), & xdisp(52), tdisp(26, 2), treac(26, 2), astif(52, 52), & aslod(52), react(52), fresv(1352), pefix(52),& estif(4, 4), strsl(250, 2) Dimension stran(2) Do ielem = 1, nelem Do ievab = 1, nevab eload(ielem, ievab) = 0.0 End Do Do idofn = l, ndofn stres(ielem, idofn) = 0.0 End Do End Do klay r = 0 Do ielem = i, nelem lprop = matno(ielem) young = props(lprop, 1) shear = props(lprop, 2) yield = props(lprop, 3) hards = props(lprop, 4) thkto = props(lprop, 5) node1 = lnods(ielem, 1) node2 = lnods(ielem, 2) eleng = abs(coord(node2)-coord(node1)) wnod1 = xdisp(node1*ndofn-1) wnod2 = xdisp(node2*ndofn-1) thta1 = xdisp(node1*ndofn) thta2 = xdisp(node2*ndofn) stran(1) = (thta1-thta2)/eleng stran(2) = (wnod2-wnod1)/eleng - 0.5*(thta1+thta2) zmidl = -thkto/2.0 kount = 5 Do ilayr = 1, nlayr klayr = klayr + 1 kount = kount + 1 brdth = props(lprop, kount) kount = kount + 1 thick = props(lprop, kount) zmidl = zmidl + thick/2.0 stlin = young*stran(1)*zmidl stcur = strsl(klayr, 1) + stlin preys = yield + hards*abs(plast(klay)) If (abs(strsl(klayr,1))>=preys) Goto 20 escur = abs(stcur) - preys If (escur0.0 .And. stlinpvalu) nchek = 999 50 pvalu = ratio Write (6, 900) iiter, nchek, ratio Return 900 Format ('0', 5X, 'ITERATION NUMBER =', I5/'0',& 5X, 'CONVERGENCE CODE =', I4, 3X, 'NORM OF RESIDUAL SUM RATIO =', E14.6) End Program !C*********************************************************** !C*********************************************************** !c result !c*********************************************************** Subroutine result !C*********************************************************** !C !C*****OUTPUTS DISPLACEMENT, REACTIONS AND STRESSES !C !C*********************************************************** Common /unim1/npoin, nelem, nboun, nload, nprop, nnode, iincs, iiter,& kresl, nchek, toler, nalgo, nsvab, ndofn, nincs, nevab,& niter, noutp, facto, pvalu Common /unim2/props(5, 4), coord(26), lnods(25, 2), ifpre(52), & fixed(52), tload(25, 4), rload(25, 4), eload(25, 4),& matno(25), stres(25, 2), plast(25), xdisp(52), & tdisp(26, 5), treac(26, 2), astif(50), aslod(52), & react(52), fresv(1352), pefix(52), estif(4, 4) If (ndofn==1) Write (6, 900) If (ndofn==2) Write (6, 910) Do ipoin = 1, npoin Write (6, 920) ipoin, (tdisp(ipoin,idofn), treac(ipoin,idofn), idofn=1, ndofn) End Do If (ndofn==2) Write (6, 930) If (ndofn==1) Write (6.940) Do ielem = i, nelem Write (6, 950) ielem(stres(ielem,idofn), idofn=1, ndofn), plast(ielem) End Do Return 900 Format ('0', 5X, 'NODE', 4X, 'DISPL.', 12X, 'REACTIONS') 910 Format (Lh0, 5X, 'NODE', 4X, 'DISPL.', 12X, 'REACTION', 7X,& 'DISPL.', 12X, 'REACTION') 920 Format (I10, 2(E14.2,5X,E14.6)) 930 Format ('0', 2X, 'ELEMENT', 12X, 'STRESSES', 12X, 'PL.STRAIN') 940 Format ('0', 2X, 'ELEMENT', 5X, 'STRESSES', 5X, 'PL.STRAIN') End Subroutine result !c******************************************************************** !c******************************************************************** !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ !********************************************************************* ! c$debug ! c$BEML !******************************************************************** !c$$$$ Input, Output: !******************************************************************** !******************************************************************** ! Modified and new routines ! Master BEML This routine is almost identical to routine BEAM dcscribed ! earlier. ! New routines for nonlayered elasto-plastic Timoshenko beam analysis ! Master BEAM The master calling routine BEAM simply organises the ! calling of the main routines as in one dimension(1D). !********************************************************************* ! SINOMACH ! ! FORTRAN95,2013&2018 version ! ! ! PotsyYZ, 周勇，Johnny, Frank, Mohammed,SwanseaUK ! ! 20230321 !********************************************************************* ! c$debug ! c$BEML !******************************************************************** !c$$$$ Input, Output: !******************************************************************** !******************************************************************** !******************************************************************** ! 广州 ! ! 20230324

（续） ******* !C !C REDUCE EQUATIONS !C pivot = astif(ieqns, ieqns) If (abs(pivot)coord(node1)) stran = (xdisp(node2)-(xdisp(node1))/eleng If(coord(node2)

感兴趣的楼主可进一步探讨。 ************************** !********************************************************************* !********************************************************************* !********************************************************************* ! MASTER UNVISC ! PROGRAM FOR THE 1-D SOLUTION OF NONLINEAR PROBLEMS ! ! FEM法则求解一维（1D）非线性弹塑性材料力学问题 ! 材料稳态热传导温度函数, 扩散材料保护膜与扩散浓度相关现象, ! 轴向杆系统特性，及取决于时间的杆系统弹塑性等 ! !********************************************************************* ! SINOMACH ! ! FORTRAN95,2013&2018 version ! ! ! PotsyYZ, 周勇，Johnny, Frank, Mohammed,SwanseaUK ! ! 20230323 !********************************************************************* !********************************************************************* ! c$debug ! c$ master unidim !c*********************************************************** ! Implicit none ! call timestamp ( ) !c*********************************************************** Program main !c*********************************************************** !c MASTER UNVISC !C*********************************************************** !C !C *** PROGRAM FOR THE 1-D SOLUTION OF NONLINEAR PROBLEMS !C !C*********************************************************** Common /unim1/npoin, nelem, nboun, nload, nprop, nnode, iincs,& istep, kresl, nchek, toler, nalgo, nsvab, ndofn, nincs, nevab, & nstep, noutp, facto, tauft, dtint, ftime, first, pvalu, dtime, ttime Common /unim2/props(5, 5), coord(26), lnods(25, 2), ifpre(52),& fixed(52), tl0ad(25, 4), rload(25, 4), eload(25, 4),& matno(25), stres(25, 2), plast(25), xdisp(52), tdisp(26, 2),& treac(26, 2), astif(52, 52), aslod(52), react(52),& fresv(1352), pefix(52), estif(4, 4), vivel(25) !c****************************************************** !c**** FINITE ELEMENTS IN PLASTICITY !c****************************************************** ttime = 0.0 Call data Call inital Call stunvp Do iincs = i, nincs Call inclod dtime = 0.0 Do istep = 1, nstep itime = ttime + dtime Call nonal Call assemb If (kresl==1) Call greduc If (kresl==2) Call resolv Call baksub Call incvp Call convp If (nchek==0) Goto 20 If (istep==1 .And. noutp==1) Call result If (noutp==2) Call result End Do Write (6, 900) Stop 20 Call result End Do Stop 900 Format ('0', 5X, 'STEADY STATE NOT ACHIEVED') End Program main !C*********************************************************** !C*********************************************************** !c data !c*********************************************************** !C*********************************************************** Subroutine data !C*********************************************************** !C !C*****INPUTS DATA DEFINING GEOMETRY,LOADING,BOUNDARY CONDITIONS...ETC. !C !C*********************************************************** Common /unim1/npoin, nelem, nboun, nload, nprop, nnode, iincs,& iiter, kresl, nchek, toler, nalgo, nsvab, ndofn, nincs, & nevab, niter, noutp, facto, pvalu Common /unim2/props(5, 4), coord(26), lnods(25, 2),& ifpre(52), fixed(52), tload(25, 4), rload(25, 4), load(25, 4),& matno(25), stres(25, 2), plast(25), xdisp(52), tdisp(26, 2), & treac(26, 2), astif(52, 52), aslod(52), react(52),& fresv(1352), pefix(52), estif(4, 4) Dimension icode(2), value(2), title(18) Read (5, 965) title Write (6, 965) title Read (5, 900) npoin, nelem, nboun, nmats, nprop, nnode, nincs, nalgo, ndofn Write (6, 905) npoin, nelem, nboun, nmats, nprop, nnode, nincs, nalgo, ndofn nevab = ndofn*nnode nsvab = ndofn*npoin Write (6, 910) Do imats = 1, nmats Read (5, 915) jmats, (props(jmats,iprop), iprop=1, nprop) Write (6, 915) jmats, (props(jmats,iprop), iprop=i, nprop) End Do Write (6, 920) Do ielem = 1, nelem Read (5, 925) jelem, (lnods(jelem,inode), inode=1, nnode), matno(jelem) Write (6, 925) jelem, (lnods(jelem,inode), inode=i, nnode), matno(jelem) End Do Write (6, 930) Do ipoin = i, npoin Read (5, 935) jpoin, coord(jpoin) Write (6, 935) jpoin, coord(jpoin) End Do Do isvab = 1, nsvab ifpre(isvab) = 0 pefix(isvab) = 0.0 End Do If (ndofn==1) Write (6, 940) If (ndofn==2) Write (6, 945) Do iboun = l, nboun Read (5, 950) nodfx, (icode(idofn), value(idofn), idofn=1, ndofn) Write (6, 950) nodfx, (icode(idofn), value(idofn), idofn=1, ndofn) nposn = (nodfx-1)*ndofn Do idofn = 1, ndofn nposn = nposn + 1 ifpre(nposn) = icode(idofn) pmfix(nposn) = value(idofn) End Do End Do Write (6, 955) Do ielem = i, nelem !C !C FINITE ELEMENTS IN PLASTICITY !C Do ievab = 1, nevab rload(ielem, ievab) = 0.0 End Do End Do 70 Read (5, 960) jelem, (rload(jelem,ievab), ievab=l, nevab) If (jelem/=nelem) Goto 70 Do ielem = 1, nelem Write (6, 960) ielem, (rload(ielem,ievab), ievab=1, nevab) End Do Return 965 Format (18A4) 900 Format (6, 9I5) 905 Format (//1X, 'NPOIN =', I5, 3X, 'NELEM =', I5,& 3X, 'NBOUN =', I5, 3X, 'NMATS =', I5//1X, 'NPROP =', I5,& 3X, 'NNODE =', I5, 3X, 'NINCS =', I5, 3X, & 'NALGO =', I5//1X, 'NDOFN =', I5) 910 Format ('0', 5X, 'MATERIAL PROPERTIES') 915 Format (I10, 4F15.5) 920 Format ('0', 3X, 'EL NODES MAT.') 925 Format (4I5) 930 Format ('0', 4X, 'NODE', 5X, 'COORD') 935 Format (I10, F15.5) 940 Format ('0', 1X, 'RES.NODE', 2X, 'CODE', 3X, 'PRES.VALUES') 945 Format ('0', 1X, 'RES.NODE', 2X, 'CODE', 3X, 'PRES.VALUES',& 2X, 'CODE', 3X, 'PRES VALUES') 950 Format (I10, 2(I5,F15.5)) 955 Format ('0', 2X, 'ELEMENT', 10X, 'NODAL LOADS') 960 Format (I10, 5F15.5) End Subroutine data !c************************************************************ !c************************************************************ !c stunvp !c************************************************************ Subroutine stunvp !C************************************************************ !C !C*** CALCULATES ELEMENT STIFFNESS MATRICES !C !C************************************************************ Common /unim1/npoin, nelem, nboun, nload, nprop, nnode, iincs, & istep, kresl, nchek, toler, nalgo, nsvab, ndofn, nincs, nevab,& nstep, noutp, facto, tauft, dtint, ftime, first, pvalu, dtime,& ttime Common /unim2/props(5, 5), coord(26), lnods(25, 2), ifpre(52), & fixed(52), tload(25, 4), rload(25, 4), eload(25, 4), matno(25),& stres(25, 2), plast(25), xdisp(52), tdisp(26, 2), treac(26, 2), & astif(52, 52), aslod(52), react(52), fresv(1352), pefix(52),& estif(4, 4), vivel(25) !C !C REWIND 1 !C Do ielem = 1, nelem lprop = matno(ielem) young = props(lprop, 1) xarea = props(lprop, 2) node1 = lnods(ielem, 1) node2 = lnods(ielem, 2) eleng = abs(coord(node1)-coord(node2)) fmult = young*xarea/eleng estif(1, 1) = fmult estif(1, 2) = -fmult estif(2, 1) = -fmult estif(2, 2) = fmult Write (1) estif End Do Return End Subroutine stunvp !C*********************************************************** !c*********************************************************** !c &&& Call assemb !C*********************************************************** Subroutine assemb !c*********************************************************** !c !c ELEMENT ASSEMBLY ROUTINE !c !c************************************************************ Common /unim1/npoin, nelem, nboun, nload, nprop, nnode, iincs, & iiter, kresl, nchek, toler, nalgo, nsvab, ndofn, nincs, nevab,& niter, noutp, facto, pvalu Common /unim2/props(5, 4), coord(26), lnods(25, 2), ifpre(52), fixed(52), tload(25, 4), rload(25, 4), eload(25, 4), matno(25), stres(25, 2), plast(25), xdisp(52), tdisp(26, 2), treac(26, 2), astif(52, 52), aslod(52), react(52), fresv(1352), pefix(52), estif(4, 4) !c************************************************************ !c !c ELEMENT ASSEMBLY ROUTINE !c !c************************************************************ !c??? REWIND 1 ???? !c************************ Do isvab = 1, nsvab aslod(isvab) = 0.0 End Do If (kresl==2) Goto 30 Do isvab = 1 nsvab Do jsvab = 1, nsvab astif(isvab, jsvab) = 0.0 End Do End Do 30 Continue !c !c ASSEMBLE THE ELEMENT LOADS !c Do ielem = 1, nelem Read (1) estif Do inode = 1, nnode nodei = lnods(ielem, inode) Do idofn = 1, ndofn nrows = (nodei-1)*ndofn + idofn nrowe = (inode-1)*ndofn + idofn aslod(nr0ws) = aslod(nrows) + eload(ielem, nrowe) !c !c ASSEMBLE THE ELEMENT STIFFNESS MATRICES !c If (kresl==2) Goto 40 Do jnode = 1, nnode nodej = lnods(ielem, jnode) Do jdofn = 1, ndofn ncols = (nodej-1)*ndofn + jdofn ncole = (jnode-1)*ndofn + jdofn astif(nrows, ncols) = astif(nrows, ncols) + estif(nrowe, ncole) 40 End Do End Do End Do End Do End Do Return End Subroutine assemb !c*********************************************************** !C*********************************************************** !c &&& greduc !C*********************************************************** Subroutine greduc !C******************************************************************* !C !C **** GAUSSIAN REDUCTION ROUTINE !C !C******************************************************************* Common /uniml/npoin, nelem, nboun, nload, nprop, nnode, iincs,& iiter, kresl, nchek, toler, nalgo, nsvab, ndofn, nincs, nevab,& niter, noutp, facto, pvalu Common /unim2/props(5, 4), coord(26), lnods(25, 2), ifpre(52),& fixed(52), tload(25, 4), rload(25, 4), eload(25, 4), matno(25),& stres(25, 2), plast(25), xdisp(52), tdisp(26, 2), treac(26, 2),& astif(52, 52), aslod(52), react(52), fresv(1352), pefix(52),& estif(4, 4) !C !C GAUSSIAN REDUCTION ROUTINE !C kount = o neqns = nsvab Do ieqns = l, neqns If (ifpre(ieqns)==l) Goto 40 !C !C REDUCE EQUATIONS !C pivot = astif(ieqns, ieqns) If (abs(pivot)coord(node1)) stran = (xdisp(node2)-(xdisp(node1))/eleng If(coord(node2)

第十章，尝试程序实践，欠佳，再争取一下： 有限元法应力变形分析大型球壳和水坝重力坝地震影响 !****************************************************************** !****************************************************************** !****************************************************************** ! c$ Master10Dynpak ! ! Stress and strain on surface of a ball, with about r=5.6m. ! ! 1. Ball shell shickness, thick=10mm ! ! 2. Waterdam, weighted, Earthquake ! !********************************************************************* ! SINOMACH ! ! FORTRAN95,2013&2018 version ! ! ! SwanseaUK ! 周勇 ! 20230317 !********************************************************************* !****************************************************************** !Program dynpak(input, tape5=input, tape4 .tape 10, tape12, tape3,& ! output, tape6=output, tape7, tape11, tape13) Program main !****************************************************************** ! !C*** DYNAMIC TRANSIENT ELASTO - -VISCOPLASTIC PRAGRAM ! !******************************************************************* ! ! Implicit none ! call timestamp ( ) Dimension acceh(600), accev(600), coord(200, 2), displ(400), & force(400), ifpre(2, 200), lnods(50, 9), matno(50),& intgr(50), nprqd(10), ngrqs(10),posgp(4), & props(10, 13), resid(400), rload(50, 18), strin(4, 450),& strsg(4, 450), tdisp(400), tempe(100), veloc(400), & vistn(4, 450),vivel(5, 450), weigp(4), ymass(400) !C Call contol(ndofn, nelem, nmats, npoin) !C Call inputd(coord, ifpre, lnods, matno, nconm, ncrit,& ndime, ndofn,nelem, ngaum, ngaus, nlaps,& nmats, nnode, npoin, nprev, nstre,ntype,& posgp, props, weigp) !C Call intime(aalfa, acceh, accev, afact, azero, beeta, & bzero, delta, dtime,dtend, gaama, ifixd,& ifunc, intgr, kstep, miter, ndofn, nelem, ngrqs,& noutd, noutp, npoin, nprqd, nreqd, & nreqs, nstep, omega, tdisp,& toler, veloc, ipred) !C Call prevos(force, ndofn, nelem, ngaus, npoin, nprev, & strin) !C Call loadpl(coord, force, lnods, matno, ndime, ndofn, nelem, ngaus,& nmats, nnode, npoin, nstre, ntype, posgp, props, rload, strin, tempe, meigp) !C Call lumass(coord, intgr, lnods, matno, nconm, ndime, ndofn, nelem, & ngaum, nmats, nnode, npoin, ntype, props, ymass) !C Call fixity(ifpre, ndofn, npoin, ymass) !C If (nprev/=0) & Call resvpl(coord, dtime, lnods, matno, ncrit, ndime,& ndofn, nelem, ngaus, nlaps, nnode, nmats,& npoin, nstre, ntype, posgp, props, resid, & rload, strin, strsg, tdisp, vistn, vivel,& weigp) !C Do istep = i, nstep !C Call explit(acceh, accev, afact, azero, aalfa, bzero,& dtime, dtend, force,ifixd, ifpre, ifunc, & istep, ndofn, npoin, omega, resid, tdisp, & veloc, ymass) !C Call resvpl(coord, dtime, lnods, matno, ncrit, ndime,& ndofn, nelem, ngaus,nlaps, nnode, nmats,& npoin, nstre, ntype, posgp, props, resid,& rload,strin, strsg, tdisp, vistn, vivel,& weigp) !C Call outdyn(displ, dtime, istep, ndofn, nelem, egaus,& ngrqs, noutd, noutp, npoin, nprqd, nreqd,& nreqs, ntype, strsg, tdisp, vivel) ! End Do ! ! c*** ! call timestamp ( ) Stop End Program main !****************************************************************** !****************************************************************** !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& !****************************************************************** subroutine blarge(bmatx, cartd, diacm, dlcod, gpcod, kgasp,& nlaps, nnode, ntype, shape) ! c******************************************************* ! c ! c***large displacement b matrix ! c ! c******************************************************** Dimension bmatx(4, 18), cartd(2, 9), djacm(2, 2), dlcod(2, 9), gpcod(2, 9), shape(9) ngash = 0 Do inode = 1, nnode mgash = ngash + 1 ngash = mgash + 1 bmatx(1, mgash) = cartd(1, inode)*djacm(1, 1) bmatx(1, ngash) = cartd(1, inode)*djacm(2, 1) bmatx(2, mgash) = cartd(2, inode)*djacm(1, 2) bmatx(2, ngash) = cartd(2, inoda)*djacm(2, 2) bmatx(3, mgash) = cartd(2, inode)*djacm(1, 1) + cartd(1, inode)*& diacm(1, 2) bmatx(3, ngash) = cartd(1, inode)*djacm(2, 2) + cartd(2, inode)*& djacm(2, 1) End Do If (ntype/=3) Return fmult = 1 If (nlaps 50) Goto 200 If (npoin >200) Goto 200 !******************************************************* !C$ ??? If (nmats > 10) Goto 200 ???? !******************************************************* If (nmats > 10) Goto 210 Goto 210 200 Write (6, 120) Stop Return 120 Format (/'SET DIMENSION EYCEEDED - CONTROL CHECK'/) 110 Format (16I5) 210 continue End Subroutine contol !******************************************************* ! 20230314 !******************************************************* !******************************************************* !c$ Subroutine explit() !******************************************************* Subroutine explit(acceh, accev, afact, azero, aalfa, bzero,& dtime, dtend,force, ifixd, ifpre, ifunc, & istep, ndofn, npoin, omega, resid, tdisp,& vetoc, ymass) !C******************************************************* !C !C*****TIME STEPPING ROUTINE !C !C******************************************************* Dimension ymass(1), acceh(1), tdisp(1), resid(1),& force(1), accev(1), veloc(1), ifpre(2, 1) cfact = 1.0 + 0.5*aalfa*dtime cfact = 1./cfact cons1 = 2.*cfact rcons = 1./cons1 cons2 = cons1 - 1 cons3 = dtime*dtime*cfact cons4 = -2.0*cons2*dtime If (istep>1) cons4 = cons2 nposn = 0 facts = functs(azero, bzero, dtend, dtime, ifunc, istep, omega) facth = functa(acceh, afact, dtend, dtime, ifunc, istep) factv = functa(accev, afact, dtend, dtime, ifunc, istep) Do ipoin = 1, npoin Do idofn = 1, ndofn factt = 0.0 If (ifunc/=0) Goto 200 If (ifixd==0 .And. idofn==1) factt = facth If (ifixd==0 .And. idofn==2) factt = factv If (ifixd==1 .And. idofn==2) factt = factv If (ifixd==2 .And. idofn==1) factt = facth If (ifpre(idofn, ipoin) == 0) Goto 200 factt = 0.0 facts = 1.0 200 Continue nposn = nposn + 1 dcurr = tdisp(nposn) resid(nposn) = resid(nposn) - force(nposn)*facts tdisp(nposn) = -resid(nposn)*cons3/ymass(nposn)& - factt*cons3 + dcurr*cons1 - veloc(nposn)*cons4 If (istep == 1) tdisp(nposn) = tdisp(nposn)*rcons veloc(nposn) = dcurr End Do End Do Return End Subroutine explit !************************************************************ !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& !************************************************************ ! 20230315 !************************************************************ !************************************************************ !************************************************************ Subroutine fixity(ifpre, ndofn, npoin, ymass) !C*********************************************************** !C !C**** DEALS WITH FIXED BOUNDARY NODES !C !C*********************************************************** Dimension ifpre(2, 1), ymass(1) ntotv = ndofn*npoin iposin = 0 !****************************************************** !c$$$$???? Do ipoin - 1, npoin ???? !****************************************************** ! 20230317 !****************************************************** Do ipoin = 1, npoin Do idofn = 1, ndofn iposn = iposn + 1 If (ifpre(idofn,ipoin)==1) ymass(iposn) = 1.E30 End Do End Do Write (6, 900) Write (6, 910)(itotv, ymass(itotv), itotv=1, ntotv) Return 900 Format (/5X, 'NODAL LUMPED MASSES'/) 910 Format (6(1X,I5,E13.5)) End Subroutine fixity !************************************************************** !************************************************************** !c$ !************************************************************** Subroutine flowvp(avect, kgaus, lprop, ncrit, nmats, props, & steff, vivel, yield) !C************************************************************* !C !C*** CALCULATES VISCOPLASTIC STRAIN RATE !C !C************************************************************** Dimension avect(4), props(nmats, 1), vivel(5, 1) If (steff==0.0) Goto 90 nstr1 = 4 eolor = 0.01 fdafm = props(lprop, 6) hards = props(lprop, 7) frict = props(lprop, 8) gamma = props(lprop, 9) delta = props(lprop, 10) nflow = props(lprop, 11) fricf = frict*0.017453292 If (ncrit==3) fdafm = fdatm*cos(frict) If (ncrit==4) fdatm = 6.0*fdatm*cos(frict)/(1.73205080757*(3.0-sin(frict))) If (hards>0.) fdatm = fdatm + vivel(5, kgaus)*hards If (fdatmdtend) Return If (ifunc==1) functs = 1.0 If (ifunc==2) argum = omega*jstep*dtime If (ifunc==2) functs = azero + bzero*sin(argum) Return End Function functs !***************************************** !c$ !***************************************** Subroutine inputd(coord, ifpre, lnods, matno, nconm, ncrit,& ndime, ndofn, nelem, ngaum, ngaus, nlaps,& nmats, nnode, npoin, nprev, nstre, ntype,& posgp, props, weigp) !C******************************************************* !C !C**** DYNPAK INPUT ROUTINE !C !C******************************************************* Dimension coord(npoin, 1), ifpre(ndofn, 1), wbigp(1),& matno(1), lnods(nelem, 1), props(nmats, 1),& posgp(1), title(10) Read (5, 913) title Write (6, 914) title !C******************************************************* !C !C****READ THE FIRST DATA CARD,AND ECHIO IT IMMEDIATELY !C !C******************************************************* Read (5, 900) nvfix, ntype, nnode, nprop, ngaus, ndime, nstre, ncrit, & nprev, nconm, nlaps, ngaum, nrads Write (5, 901) npoin, nelem, nvfix, ntype, nnode, ndofn, nmats, & nprop, ngaus, ndime, nstre, ncrit, nprev, nconm,& nlaps, ngaum, nrads !C !C***** READ THE ELEMENT NODAL CONNECTIONS,AND THE PROPERTY NUMBERS !C Write (6, 902) Do ielem = 1, nelem Read (5, 900) numel, matno(numel), (lnods(numel,inode),& inode=1, nnode) Write (13, 915) numel, (lnods(numel,inode), inode=1, nnode) Write (6, 903) numel, matno(numel), (lnods(numel,inode), inode=1, nnode) End Do !C !C*** ZERO ALL THE NODAL COORDINATES,PRIOR TO READING SOME OF THEM !C Do ipoin = 1, npoin Do idime = 1, ndime coord(ipoin, idime) = 0.0 End Do End Do 200 Read (5, 905) ipoin, (coord(ipoin,idime), idime=1, ndime) Write (6, 906) ipoin, (coord(ipoin,idime), idime=1, ndime) If (ipoin/=npoin) Goto 200 !C !C*** INTERPOLATE COORDINATES OF MID-SIDE NODES !C !*************************************************************************** !*************************************************************************** !c$$$$ ???? all nodxyr(coord, lnods, nelem, nnode, npoin, nrads, ntype) ???? !*************************************************************************** ! 20230317 !*************************************************************************** call nodxyr(coord, lnods, nelem, nnode, npoin, nrads, ntype) !C Write (6, 904) Write (13, 916)(ipoin, (coord(ipoin,idime),idime=1,ndime), ipoin=1, npoin) Write (13, 906)(ipoin, (coord(ipoin,idime),idime=1,ndime), ipoin=1, npoin) !C !C***READ THE FIXED VALUES !C Write (6, 907) Do ipoin = 1, npoin Do idofn = 1, ndofn ifpre(idofn, ipoin) = 0 End Do End Do Do ivfix = 1, nvfix Read (5, 908) ipoin, (ifpre(idofn,ipoin), idofn, ndofn) End Do Do ipoin = 1, npoin !********************************************************************** !c $$$ ??? Write (6, 909) ipoin, (ifpre(idofn,ipoin), idofn, ndofn)???? ! 20230317 !********************************************************************** Write (6, 909) ipoin, ifpre(idofn,ipoin), idofn, ndofn End Do !C !C***READ THE AVAILABLE SELECTION OF ELEMENT PROPERTIES !C !c Write (6, 910) Do imats = 1, nmats Read (5, 900) numat Read (6, 917)(props(numat,iprop), iprop=1, nprop) Write (6.911) numat Write (6, 912)(props(numat,iprop), iprop=1, nprop) End Do !C !C****SET UP GANSSIAN INTERATION CONSTANTS !C Call gaussq(ngaus, posgp, weigp) Return 913 Format (10A4) 914 Format (//, 5X, 10A4) 901 Format (/5X, 'CONTROL PARAMETERS'/& /5X, 'NPOIN=', I10, 5X, 'NELEM=', I10, 5X, 'NVFIX=', I10/& /5X, 'NTYPE=', I10, 5X, 'NNODE=', I10, 5X, 'NDOFN=', I10/& /5X, 'NMATS=', I10, 5X, 'NPROP=', I10, 5X, 'NGAUS=', I10/& /5X, 'NDIME=', I10, 5X, 'NSTRE=', I10, 5X, 'NCRIT=', I10/& /5X, 'NPREV=', I10, 5X, 'NCONM=', I10, 5X, 'NLAPS=', I10/& /5X, 'NGAUM=', I10, 5X, 'NRADS=', I10/) 900 Format (16I5) 902 Format (//5X, 'ELEMENT', 3X, 'PROPERTY', 6X, 'NODE NUMBERS') 903 Format (6X, I5, I9, 6X, 10I5) 915 Format (16I5) !C !C !C****READ SOME NODAL COORDIATES,FINISHING WITH THE LAST !C NODE OF ALL !C 904 Format (//5X, 'NODE', 9X, 'X', 9X, 'Y', 5X) 905 Format (I5, 6F10.5) !********************************************************** !c$$$ ??? 916 FORMAT(I ,2G15.6) ??? !********************************************************** ! 20230314 !********************************************************** 916 Format (I5, 2G15.6) 906 Format (5X, I5, 2F10.3) 907 Format (//5X, 'NODE', 2X, 'CODE') 908 Format (1X, I4, 3X, 2I1) 909 Format (6X, I5, 3X, 2I1) 910 Format (//5X, 'MATHRIAL PROPERTIES') 911 Format (/5X, 'MATERIAL NO', I5) !********************************************************** !c$$$ ??? 912 FORMAT (/5,'YOUNG MODULUS',G12.4/5X,'POISSON RATIO',G12.4/& ! ........................... ??? !********************************************************** ! 20230314 !********************************************************** 912 Format (/5X, 'YOUNG MODULUS', G12.4/5X, 'POISSON RATIO', G12.4/& 5X, 'THCKNESS', G12.4/5X, 'MASS DENSITY', G12.4/& 5X, 'ALPHA TEMPR', G12.4/5X, 'REFERENCE FO', G12.4/& 5X, 'HARDENING PAR', G12.4/5x, 'FRICT ANGLE', G12.4/& 5X, 'FLUDITY PAR', G12.4/5X, 'EXP DELTA', G12.4/& 5X, 'NFLOW CODE', G12.4) 917 Format (8E10.4) End Subroutine inputd !*************************************************************** !*************************************************************** ! 20230316 !*************************************************************** Subroutine intime(aalfa, acceh, accev, afact, azero, beeta, & bzero, delta, dtime, dtend, gaama, ifixd,& ifunc, intgr, kstep, miter, ndofn, nelem,& ngrqs, noutd, koutp, npoin, nprqd, nreqd, & nreqs, nstep, omega, tdisp, toler, veloc,& ipred) !C********************************************************************************** !C !C*** INITIAI VALUES AND TIME INTEGRATION DATA !C !C*********************************************************************************** Dimension tdisp(1), acceh(1), nprod(1), intgr(1),& veloc(1), accev(1), ngrqs(1) !C******************************************************* !C !C**** READ TIME STEPPING AND SELECTIVE OUTPUT PARAMETERS !C !C******************************************************* Read (5, 902) nstep, noutd, noutp, nreqd, nreqs, nacce, ifunc,& ifixd, miter, kstep, ipred Read (5.190) dtime, dtend, dtrec, aalfa, beeta, delta, gaama,& azero, bzero, omega, toler Write (6, 950) nstep, noutd, noutp, nreqd, nrrqs, nacce, ifunc,& ifixd, miter, kstep, ipreo Write (6, 960) dtime, dtend, dtrec, aalfa, beeta, delta, gaama, & azero, bzero, omega, toler !C !C*** SELECTED NODES AND GAUSS POINTS FOR OUTPUT !C Read (5, 902)(nprqd(ireqd), ireqd=1, nreqd) Read (5, 902)(ngrqs(ireqs), ireqs=1, nreqs) Write (6, 909) Write (6, 910)(nprqd(ireqd), ireqd=1, nreqd) Write (6, 911)(ngrqs(ireqs), ireqs=1, nreqs) !C !C** READ THE INDICATOR FOR EXPLICIT OR IMPLICIT ELEMENT !C Read (5, 902)(intgr(ielem), ielem=1, nelem) Write (6, 930) Write (6, 902)(intgr(ielem), ielem=1, nelem) !C !C*** INITIAL DISPLACEMENTS !C jpoin = 0 Do ipoin = 1, npoin Do idofn = 1, ndofn jpoin = jpoin + 1 tdisp(jpoin) = 0. veloc(jpoin) = 0. End Do End Do Write (6, 903) 200 Read (5, 904) ngash, xgash, ygash nposn = (ngash-1)*ndofn+1 tdisp(nposn) = xgash nposn = nposn+1 tdisp(nposn) = ygash Write (6, 905) ngash, xgash, ygash If (ngash/=npoin) Goto 200 !C !C***** INITIAL VELOCITIES !C Write (6, 906) 210 Read (5, 904) ngash, xgash, ygash nposn = (ngash-1)*ndofn + 1 veloc(nposn) = xgash nposn = nposn + 1 veloc(nposn) = ygash Write (6, 905) ngash, xgash, ygash If (nfash/=npoin) Goto 210 If (ifunc/=0) Goto 250 !C !C****READ ACCELEROGRAM DATA,X-DIREC FROM TAPE 7,Y-DIREC FROM TAPE 12 !C afact = dtrec/dtime If (ifixd-1) 220, 230, 240 220 Read (7, 907)(acceh(i), i=1, nacce) Write (6, 912) dtrec Write (6, 907)(acceh(i), i=1, nacce) Read (12, 907)(accev(i), i=1, nacce) Write (6, 913) dtrec Write (6, 907)(accev(i), i=1, nacce) Goto 250 230 Read (12, 907)(adcev(i), i=1, nacce) Write (6, 913) dtrec Write (6, 907)(adcev(i), i=1, nacce) Goto 250 240 Read (7, 907)(acceh(i), i=1, nacce) Write (6, 912) Write (6, 907)(acceh(i), i=1, nacce) 250 Continue Return 950 Format (/5X, 'TIME STEPPING PARAMETERS'/& /5X, 'NSTEP=', I5, 12X, 'NOUTD=', I5, 12X, 'NOUTP=', I5,/& /5X, 'NREQD=', I5, 12X, 'NREQS=', I5, 12X, 'NACCE=', I5,/& /5X, 'IFUNC=', I5, 12X, 'IFIXD=', I5, 12X, 'MITER =', I5,/& /5X, 'KSTEP=', I5, 12X, 'IPRED=', I5) 960 Format (/5X, 'DTIME =', G12.4, 5X, 'DTEND=', G12.4, 5X, 'DTREC=', G12.4/& /5X, 'AALFA=', G12.4, 5X, 'BEETA=', G12.4, 5X, 'DELTA=', G12.4/& /5X, 'GAAMA=', G12.4, 5X, 'AZERO=', G12.4, 5X, 'BZERO=', G12.4/& /5X, 'OMEGA=', G12.4, 5X, 'TOLER=', G12.4) 909 Format (//5X, 'SELECTIPE OUTPUT REQUESTED FOR FOLLOWING') 910 Format (/5X, 'NODES', 10I5) 911 Format (5X, 'G.P.', 10I5) 902 Format (16I5) 190 Format (8F10.4) 930 Format (/5X, 'TYPE OF ELEMENT,IMPLICIT=1,EXPLICIT=2'/) 904 Format (I5, 2F10.5) 903 Format (//5X, 'NODE', 2X, 'INITIAL X-DISP.', 2X,& 'INITIAL Y-DISP.'/) 905 Format (I10, 2E18.5) 906 Format (//5X, 'NODE', 2X, 'INITIAL X-VELO.', 2X,& 'INITIAL Y-VELO.'/) 907 Format (7F10.3) 912 Format (/5X, 'HORIZONTAL ACCELERATION ORDINATES AT', F9.4,& 2X,'SEC'/) 913 Format (/5X, 'VERTICAL ACCELERATION ORDINATES AT', F9.4,& 2X,'SEC'/) End Subroutine intime !****************************************************************** !****************************************************************** ! 20230316 !****************************************************************** Subroutine invar(devia, lprop, ncrit, nmats, props, sint3,& steff, stemp, theta, varj2, yield) !C******************************************************* !C !C*** STRESS INVARIANTS !C !C******************************************************* Dimension devia(4), props(nmats, 1), stemp(4) !C !C*** INVARIANTS !C root3 = 1.73205080757 smean = (stemp(1)+stemp(2)+stemp(4))/3.0 devia(1) = stemp(1) - smean devia(2) = stemp(2) - smean devia(3) = strmp(3) devia(4) = stemp(4) - smean varj2 = devia(3)*devia(3) + 0.5*devia(1)*devia(1) + & devia(2)*devia(2) + devia(4)*devia(4)) varj3 = devia(4)*(devia(4)*devia(4) - varj2) steff = sqrt(varj2) If (varj2 == 0.0 .Or. steff == 0.0) Goto 5 sint3 = -2.580762113*varj3/(varj2*steff) Goto 6 5 sint3 = 0.0 6 Continue If (sint3 < -1.0) sint3 = -1.0 If (sint3 > 1.0) sint3 = 1.0 theta = asin(sint3)/3.0 Goto (1, 2, 3, 4) ncrit !C*** TRESCA 1 yield = 2.0*cos(theta)*steff Return !C*** VON MISES 2 yield = root3*steff Return !C***MOHR-COULOMB 3 phira = props(lprop, 8)*0.017453292 snphi = sin(phira) !********************************************************* !c $$$$??? yield = smean*snphi + steff*(cos(theta) - sin(theta))*& ! snphi/root3 ???? ! !********************************************************* ! 20230317 !********************************************************* yield = smean*snphi + steff*(cos(theta) - sin(theta))*& snphi/root3 Return !C !C*** DRUCKER-PRAGER !C 4 phira = props(lprops, 8)*0.017453292 snphi = sin(phira) yield = 6.0*smean*snphi/(root3*(3.0-snphi)) + steff Return End Subroutine invar !************************************************************* !************************************************************* ! 20230316 !************************************************************* Subroutine jacobd(cartd, dlcod, djacm, ndime, nlaps, nnode) !C************************************************************* !C !C**** DEFORMATION JACOBIAN !C !C************************************************************** Dimension cartd(2, 9), dlcod(2, 9), djacm(2,2) If (nlaps>1) Goto 10 !C !C*** FOR SMALL DISPLACEMENT !C djacm(1, 1) = 1.0 djacm(2, 2) = 1.0 djacm(1, 2) = 0.0 djacm(2, 1) = 0.0 Return !C !C*** FOR LARGE DISPLACEMENT !C 10 Continue Do idime = 1, ndime Do jdime = 1, ndime djacm(idime, jdime) = 0.0 Do inode = 1, nnode djacm(idime, jdime) = djacm(idime, jdime) & + dlcod(idime, inode)*cartd(jdime, inode) End Do End Do End Do Return End Subroutine jacobd !********************************************************************* !********************************************************************* ! 20230316 !********************************************************************* !c$ LINGNL !********************************************************************* Subroutine lingnl(cartd, djacm, dmatx, eldis, gpcod, kgasp, & kgaus, ndofn, nlaps, nnode, nstre, ntype,& poiss, shape, stran, stres, strin) !C********************************************************************* !C !C*** ELASTIC STRAIN AND STRESSES !C !C********************************************************************** Dimension cartd(2, 9), stran(4), dmatx(4, 4), strin(4, 1),& eldis(2, 9), stres(4), djacm(2, 2), agash(2,2), & gpcod(2, 9), shape(9) !C !C*** CALCULATE STRAINS FROM DEFORMATION JACOBIAN !C If (nlaps

第七章，程序，欠佳： !********************************************************************* !********************************************************************* !********************************************************************* ! FEAM (Finite Element Analysis Method) Analysis for ! Static,Quasi-static,and Dynamic Forces of ! Stress, Shear and Strain on Surfaces of or within ! Loaded Materials with Deformation and ! Plastic and Elastic Solid with Ends and Joints in ! 1,2,3D, and M Dimensions in ! transient state or a fixed time of period. !********************************************************************* ! SINOMACH ! ! FORTRAN95,2013&2018 version ! ! ! PotsyYZ, 周勇，Earl, Coopper, Simzer ! ! 20230301 !********************************************************************* ! c$debug ! c$large ! Program master plast Program main ! c**************************************************************** ! c Program for the elasto - plastic analysis of plane stress, ! c plane strain end axisymmetric solids. ! c**************************************************************** ! ! Implicit none ! call timestamp ( ) Dimension asdis(300), coord(150, 2), eload(40, 18), estif(18, 18),& eqrhs(10), equat(80, 10), fixed(300), gload(80), gstif(3240),& iffix(300), lnods(40, 9), locel(18), matno(40), nacva(80),& namev(10), ndest(18), ndfro(40), nofix(30), noutp(2),& npivo(10), posgp(4), presc(30, 2), props(5, 7), rload(40, 18),& stfor(300), treac(30, 2), vecrv(80), weigp(4), strsg(4, 360),& tdisp(300), tload(40, 18), tofor(300), epstn(360), effst(360) ! c*** Open (5, File='data\input.dat', Status='unknown') Open (6, File='data\output.dat', Status='unknown') ! c*** ! c ! c***Preset variables associated with dynamic dimensioning. ! c Call dimen(mbufa, melem, mevab, mfron, mmats, mpoin, mstif, mtotg,& mtotv, mvfix, ndofn, nprop, nstre) ! c ! c***Call the subriutine which reads most of problem data. ! c Call input(coord, iffix, lnods, matno, melem, mevab, mfron, mmats,& mpoin, mtotv, mvfix, nalgo, ncrit, ndfro, ndofn, nelem, nevab,& ngaus, ngau2, nincs, nmats, nnode, nofix, npoin, nprop, nstre,& nstr1, ntotg, ntotv, ntype, nvfix, posgp, presc, props, weigp) ! c ! c***Call the subroutine which computes the consistent load vectors ! c for each element after reading the relevent input data. ! c Call loadps(coord, lnods, matno, melem, mmats, mpoin, nelem, nevab,& ngaus, nnode, npoin, nstre, ntype, posgp, props, rload, weigp, ndofn) ! c ! c***Initialise certain arrays. ! c Call zero(eload, melem, mevab, mpoin, mtotg, mtotv, ndofn, nelem, & nevab, ngaus, nstr1, ntotg, epstn, effst, ntotv, nvfix, strsg,& tdisp, tfact, tload, treac, mvfix) ! c ! c***Loop over each increment. ! c Do iincs = 1, nincs ! c ! c***Read data for current increment. ! c Call increm(eload, fixed, iincs, melem, mevab, miter, mtotv, & mvfix, ndofn, nelem, nevab, noutp, nofix, ntotv, nvfix, presc,& rload, tfact, tload, toler) ! c ! c***Loop over each iteration. ! c Do iiter = 1, miter ! c ! c***Call routine which selects solution alorithm variable kresl. ! c Call algor(fixed, iincs, iiter, kresl, mtotv, nalgo, ntotv) ! c ! c***check whether a new evaluation of the stiffness matrix is required ! c If (kresl==1) Call stiffp(coord, epstn, iincs, lnods, matno,& mevab, mmats, mpoin, mtotv, nelem, nevab, ngaus, nnode, nstre,& nstr1, posgp, props, weigp, melem, mtotg, strsg, ntype, ncrit) ! c ! c***Solve equations. ! c Call front(asdis, eload, eqrhs, equat, estif, fixed, iffix, & iincs, iiter, gload, gstif, locel, lnods, kresl, mbufa, melem,& mevab, mfron, mstif, mtotv, mvfix, nacva, namev, ndest, ndofn, & nelem, nevab, nnode, nofix, npivo, npoin, ntotv, tdisp, tload,& treac, vecrv) ! c ! c***Calculate residual forces. ! c Call residu(asdis, coord, effst, eload, facto, iiter, lnods, & lprop, matno, melem, mmats, mpoin, mtotg, mtotv, ndofn, nelem,& nevab, ngaus, nnode, nstr1, ntype, posgp, props, nstre,& ncrit, strsg, weigp, tdisp, epstn) ! c ! c***Check for convergence. ! c Call conver(eload, iiter, lnods, melem, mevab, mtotv, nchek, & ndofn, nelem, nevab, nnode, ntotv, pvalu, stfor, tload, tofor, toler) ! c ! c***Output results if required. ! c If (iiter==1 .And. noutp(1)>0) Call output(iiter, mtotg, mtotv,& mvfix, nelem, ngaus, nofix, noutp, npoin, nvfix, strsg, tdisp,& treac, epstn, ntype, nchek) ! c ! c***If solution has converged stop iterating and output results. ! c If (nchek==0) Goto 75 End Do ! c ! c*** ! c If (nalgo==2) Goto 75 Stop 75 Call output(iiter, mtotg, mtotv, mvfix, nelem, ngaus, nofix,& noutp, npoin, nvfix, strsg, tdisp, treac, epstn, ntype, nchek) End Do ! c*** Close (5) Close (6) ! c*** ! call timestamp ( ) Stop ! End Program master End Program main !***************************************************** !***************************************************** !c$ 20230227 !***************************************************** !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& !***************************************************** ! c$debug ! c$large Subroutine algor(fixed, iincs, iiter, kresl, mtotv, nalgo, ntotv) ! c****************************************************************** ! c ! c***This subroutine sets equation resolution index, krest. ! c ! c****************************************************************** Dimension fixed(mtotv) kresl = 2 If (nalgo==1 .And. iincs==1 .And. iiter==1) kresl = 1 If (nalgo==2) kresl = 1 If (nalgo==3 .And. iiter==1) kresl = 1 If (nalgo==4 .And. iincs==1 .And. iiter==1) kresl = 1 If (nalgo==4 .And. iiter==2) kresl = 1 If (iiter==1) Return Do itotv = 1, ntotv fixed(itotv) = 0.0 End Do Return End Subroutine algor !*************************************************************************** ! ! c$debug ! c$large Subroutine check1(ndofn, nelem, ngaus, nmats, nnode, npoin, nstre, ntype, & nvfix, ncrit, nalgo, nincs) ! ! c********************************************************************** ! c ! c***This subroutine checks the main control data. ! c ! c********************************************************************** Dimension neror(24) Do ieror = 1, 12 neror(ieror) = 0 End Do !End Do ! c ! c***Create the diagnostic messages. ! c If (npoin0) kount = 1 End Do If (kount==0) neror(23) = neror(23) + 1 kvfix = ivfix - 1 Do jvfix = 1, kvfix If (ivfix/=1 .And. nofix(ivfix)==nofix(jvfix)) neror(24) = neror(24) + 1 End Do End Do keror = 0 Do ieror = 13, 24 If (neror(ieror)==0) Goto 180 keror = 1 Write (6, 910) ieror, neror(ieror) Write (*, 910) ieror, neror(ieror) 180 Continue If (keror/=0) Goto 200 ! c ! c***Return all nodal connection numbers to positive values. ! c Do ielem = 1, nelem Do inode = 1, nnode lnods(ielem, inode) = iabs(lnods(ielem,inode)) End Do End Do Return 200 Call echo 900 Format (/'Check Why Node', I4, 'Never Appears') 905 Format (//'Maximum Front Width Encountered =', I5) 910 Format (//'***Diagnosis by check2, Error', I3, 6X, 'Associated Number', I5) End Do End Do End Do End Subroutine check2 !c !************************************************************************ Subroutine nodexy(coord, lnods, melem, mpoin, nelem, nnode) ! ! c********************************************************************** ! c ! c****This subroutine interpolates the mid side nodes of straight ! c sides of elements and the central node of 9 noded elements. ! c ! c********************************************************************** Dimension coord(mpoin, 2), lnods(melem, 9) If (nnode==4) Return ! c ! c***Loop over each element. ! c Do ielem = 1, nelem ! c ! c***Loop over each element edge. ! c nnod1 = 9 If (nnode==8) nnod1 = 7 Do inode = 1, nnod1, 2 If (inode==9) Goto 50 ! c ! c***Compute the node number of the first node. ! c nodst = lnods(ielem, inode) igash = inode + 2 If (igash>8) igash = 1 ! c ! c***Compute the node number of the last node. ! c nodfn = lnods(ielem, igash) midpt = inode + 1 ! c ! c***Compute the node number of the intermediate node. ! c nodmd = lnods(ielem, midpt) total = abs(coord(nodmd,1)) + abs(coord(nodmd,2)) ! c ! c***If the coordinates of the intermediate node are both ! c zero interpolate by a straight line. ! c If (total>0.0) Goto 20 kount = 1 10 coord(nodmd, kount) = (coord(nodst,kount)+coord(nodfn,kount))/2.0 kount = kount + 1 If (kount==2) Goto 10 20 End Do Goto 30 50 lnode = lnods(ielem, inode) total = abs(coord(lnode,1)) + abs(coord(lnode,2)) If (total>0.0) Goto 30 lnod1 = lnods(ielem, 1) lnod3 = lnods(ielem, 3) lnod5 = lnods(ielem, 5) lnod7 = lnods(ielem, 7) kount = 1 40 coord(lnode, kount) = (coord(lnod1,kount)+coord(lnod3,kount)+& coord(lnod5,kount)+coord(lnod7,kount))/4.0 kount = kount + 1 If (kount==2) Goto 40 30 End Do Return End Subroutine nodexy ! c !************************************************************************ Subroutine gaussq(ngaus, posgp, weigp) ! ! c******************************************************************* ! c ! c***Sets up the gauss - legendre integration constants ! c ! c******************************************************************* Dimension posgp(4), weigp(4) If (ngaus>2) Goto 4 posgp(1) = -0.577350269189626 weigp(1) = 1.0 Goto 6 4 posgp(1) = -0.774596669241483 posgp(2) = 0.0 weigp(1) = 0.5555555555555556 weigp(2) = 0.8888888888888889 6 kgaus = ngaus/2 Do igash = 1, kgaus jgash = ngaus + 1 - igash posgp(jgash) = -posgp(igash) weigp(jgash) = weigp(igash) End Do Return End Subroutine gaussq !************************************************* !c Subroutine sfr2(deriv, etasp, exisp, nnode, shape) ! ! c****************************************************************** ! c ! c***This subroutinue evaluates shape functions and their ! c derivertives for linear quadratic lagrangian and ! c serendipity isoparametric 2 - d elements. ! c ! c******************************************************************* Dimension deriv(2, 9), shape(9) s = exisp t = etasp If (nnode>4) Goto 10 st = s*t ! c ! c***Shape function for 4 noded element. ! c shape(1) = (1-t-s+st)*0.25 shape(2) = (1-t+s-st)*0.25 shape(3) = (1+t+s+st)*0.25 shape(4) = (1+t-s-st)*0.25 ! c ! c***Shape function derivertives. ! c deriv(1, 1) = (-1+t)*0.25 deriv(1, 2) = (+1-t)*0.25 deriv(1, 3) = (+1+t)*0.25 deriv(1, 4) = (-1-t)*0.25 deriv(2, 1) = (-1+s)*0.25 deriv(2, 2) = (-1-s)*0.25 deriv(2, 3) = (+1+s)*0.25 deriv(2, 4) = (+1-s)*0.25 Return 10 If (nnode>8) Goto 30 s2 = s*2.0 t2 = t*2.0 ss = s*s tt = t*t st = s*t sst = s*s*t stt = s*t*t st2 = s*t*2.0 ! c ! c***Shape functions for 8 noded element. ! c shape(1) = (-1.0+st+ss+tt-sst-stt)/4.0 shape(2) = (1.0-t-ss+sst)/2.0 shape(3) = (-1.0-st+ss+tt-sst+stt)/4.0 shape(4) = (1.0+s-tt-stt)/2.0 shape(5) = (-1.0+st+ss+tt+sst+stt)/4.0 shape(6) = (1.0+t-ss-sst)/2.0 shape(7) = (-1.0-st+ss+tt+sst-stt)/4.0 shape(8) = (1.0-s-tt+stt)/2.0 ! c ! c***Shape function derivertives. ! c deriv(1, 1) = (t+s2-st2-tt)/4.0 deriv(1, 2) = -s + st deriv(1, 3) = (-t+s2-st2+tt)/4.0 deriv(1, 4) = (1.0-tt)/2.0 deriv(1, 5) = (t+s2+st2+tt)/4.0 deriv(1, 6) = -s - st deriv(1, 7) = (-t+s2+st2-tt)/4.0 deriv(1, 8) = (-1.0+tt)/2.0 deriv(2, 1) = (s+t2-ss-st2)/4.0 deriv(2, 2) = (-1.0+ss)/2.0 deriv(2, 3) = (-s+t2-ss+st2)/4.0 deriv(2, 4) = -t - st deriv(2, 5) = (s+t2+ss+st2)/4.0 deriv(2, 6) = (1.0-ss)/2.0 deriv(2, 7) = (-s+t2+ss-st2)/4.0 deriv(2, 8) = -t + st Return 30 Continue ss = s*s st = s*t tt = t*t s1 = s + 1.0 t1 = t + 1.0 s2 = s*2.0 t2 = t*2.0 s9 = s - 1.0 t9 = t - 1.0 ! c ! c***Shape function for 9 noded element. ! c shape(1) = 0.25*s9*st*t9 shape(2) = 0.5*(1.0-ss)*t*t9 shape(3) = 0.25*s1*st*t9 shape(4) = 0.5*s*s1*(1.0-tt) shape(5) = 0.25*s1*st*t1 shape(6) = 0.5*(1.0-ss)*t*t1 shape(7) = 0.25*s9*st*t1 shape(8) = 0.5*s*s9*(1.0-tt) shape(9) = (1.0-ss)*(1.0-tt) !c !c***Shape function derivertives. !c deriv(1, 1) = 0.25*t*t9*(-1.0+s2) deriv(1, 2) = -st*t9 deriv(1, 3) = 0.25*(1.0+s2)*t*t9 deriv(1, 4) = 0.5*(1.0+s2)*(1.0-tt) deriv(1, 5) = 0.25*(1.0+s2)*t*t1 deriv(1, 6) = -st*t1 deriv(1, 7) = 0.25*(-1.0+s2)*t*t1 deriv(1, 8) = 0.5*(-1.0+s2)*(1.0-tt) deriv(1, 9) = -s2*(1.0-tt) deriv(2, 1) = 0.25*(-1.0+t2)*s*s9 deriv(2, 2) = 0.5*(1.0-ss)*(-1.0+t2) deriv(2, 3) = 0.25*s*s1*(-1.0+t2) deriv(2, 4) = -st*s1 deriv(2, 5) = 0.25*s*s1*(1.0+t2) deriv(2, 6) = 0.5*(1.0-ss)*(1.0+t2) deriv(2, 7) = 0.25*s*s9*(1.0+t2) deriv(2, 8) = -st*s9 deriv(2, 9) = -t2*(1.0-ss) !****************************************** !c$??? 20 CONTINUE ???? 20230227 !****************************************** 20 Continue Return End Subroutine sfr2 !*********************************************************** !c Subroutine jacob2(cartd, deriv, djacb, elcod, gpcod, ielem, kgasp, nnode, shape) ! ! c******************************************************************** ! c ! c***this subroutine evaluates the jacobian matrix and the ! c cartesian shape function derivertives. ! c ! c******************************************************************** Dimension cartd(2, 9), deriv(2, 9), elcod(2, 9), gpcod(2, 9), & shape(9), xjaci(2, 2), xjacm(2, 2) ! c ! c***Calculate coordinates of sampling point. ! c Do idime = 1, 2 gpcod(idime, kgasp) = 0.0 Do inode = 1, nnode gpcod(idime, kgasp) = gpcod(idime, kgasp) + & elcod(idime, inode)*shape(inode) End Do End Do ! c ! c***Create jacobian matrix xjacm. ! c Do idime = 1, 2 Do jdime = 1, 2 xjacm(idime, jdime) = 0.0 Do inode = 1, nnode xjacm(idime, jdime) = xjacm(idime, jdime) +& deriv(idime, inode)*elcod(jdime, inode) End Do End Do End Do ! c ! c***Calculate determinant and inverse of jacobian matrix ! c djacb = xjacm(1, 1)*xjacm(2, 2) - xjacm(1, 2)*xjacm(2, 1) If (djacb) 6, 6, 8 6 Write (*, 600) ielem Stop 8 Continue xjaci(1, 1) = xjacm(2, 2)/djacb xjaci(2, 2) = xjacm(1, 1)/djacb xjaci(1, 2) = -xjacm(1, 2)/djacb xjaci(2, 1) = -xjacm(2, 1)/djacb !c !c***Calculate cartesian derivertives !c Do idime = 1, 2 Do inode = 1, nnode cartd(idime, inode) = 0.0 Do jdime = 1, 2 cartd(idime, inode) = cartd(idime, inode) +& xjaci(idime, jdime)*deriv(jdime, inode) continue End Do End Do End Do Return 600 Format (//, 'PROGRAM HALTED IN SUBROUTINE JACOB2', /,& 11X, 'ZERO OR NEGATIVE AREA', /, 10X, 'ELEMENT NUMBER', I5) End Subroutine jacob2 !**************************************************************** ! c Subroutine bmatps(bmatx, cartd, nnode, shape, gpcod, ntype, kgasp) ! ! c********************************************************************** ! c ! c***This subroutine evaluates the strain - displacement matrix. ! c ! c********************************************************************** Dimension bmatx(4, 18), cartd(2, 9), shape(9), gpcod(2, 9) ngash = 0 Do inode = 1, nnode mgash = ngash + 1 ngash = mgash + 1 bmatx(1, mgash) = cartd(1, inode) bmatx(1, ngash) = 0.0 bmatx(2, mgash) = 0.0 bmatx(2, ngash) = cartd(2, inode) bmatx(3, mgash) = cartd(2, inode) bmatx(3, ngash) = cartd(1, inode) If (ntype/=3) Goto 10 bmatx(4, mgash) = shape(inode)/gpcod(1, kgasp) bmatx(4, ngash) = 0.0 10 End Do Return End Subroutine bmatps !************************************************************************ !c Subroutine modps(dmatx, lprop, mmats, ntype, props) ! ! c********************************************************************* ! c ! c***This subroutine evaluates the d - matrix. ! c ! c********************************************************************* Dimension dmatx(4, 4), props(mmats, 7) young = props(lprop, 1) poiss = props(lprop, 2) Do istr1 = 1, 4 Do jstr1 = 1, 4 dmatx(istr1, jstr1) = 0.0 End Do End Do If (ntype/=1) Goto 4 ! c ! c***d matrix for plane stress case. ! c const = young/(1.0-poiss*poiss) dmatx(1, 1) = const dmatx(2, 2) = const dmatx(1, 2) = const*poiss dmatx(2, 1) = const*poiss dmatx(3, 3) = (1.0-poiss)*const/2.0 Return 4 If (ntype/=2) Goto 6 ! c ! c***d Matrix for plane strain case. ! c const = young*(1.0-poiss)/((1.0+poiss)*(1.0-2.0*poiss)) dmatx(1, 1) = const dmatx(2, 2) = const dmatx(1, 2) = const*poiss/(1.0-poiss) dmatx(2, 1) = const*poiss/(1.0-poiss) dmatx(3, 3) = (1.0-2.0*poiss)*const/(2.0*(1.0-poiss)) Return 6 If (ntype/=3) Goto 8 !c !c***d Matrix for axisymmetric case. !c const = young*(1.0-poiss)/((1.0+poiss)*(1.0-2.0*poiss)) conss = poiss/(1.0-poiss) dmatx(1, 1) = const dmatx(2, 2) = const dmatx(3, 3) = const*(1.0-2.0*poiss)/(2.0*(1.0-poiss)) dmatx(1, 2) = const*conss dmatx(1, 4) = const*conss dmatx(2, 1) = const*conss dmatx(2, 4) = const*conss dmatx(4, 1) = const*conss dmatx(4, 2) = const*conss dmatx(4, 4) = const 8 Continue Return End Subroutine modps ! c Subroutine dbe(bmatx, dbmat, dmatx, mevab, nevab, nstre, nstr1) ! ! c********************************************************************** ! c ! c***This subroutine multiplies the d - matrix by b - matrix ! c ! c********************************************************************** Dimension bmatx(nstr1, mevab), dbmat(nstr1, mevab), dmatx(nstr1, nstr1) Do istre = 1, nstre Do ievab = 1, nevab dbmat(istre, ievab) = 0.0 Do jstre = 1, nstre dbmat(istre, ievab) = dbmat(istre, ievab) + & dmatx(istre, jstre)*bmatx(jstre, ievab) End Do End Do End Do Return End Subroutine dbe !*************************************************************** !*************************************************************** ! c$debug ! c$large Subroutine conver(eload, iiter, lnods, melem, mevab, mtotv, nchek, ndofn,& nelem, nevab, nnode, ntotv, pvalu, stfor, tload, tofor, toler) ! ! c********************************************************************** ! c ! c***this subroutine checks for convergence of the iteration process ! c ! c********************************************************************** Dimension eload(melem, mevab), lnods(melem, 9), stfor(mtotv), & tofor(mtotv), tload(melem, mevab) nchek = 0 resid = 0.0 retot = 0.0 remax = 0.0 Do itotv = 1, ntotv stfor(itotv) = 0.0 tofor(itotv) = 0.0 End Do Do ielem = 1, nelem kevab = 0 Do inode = 1, nnode locno = iabs(lnods(ielem,inode)) Do idofn = 1, ndofn kevab = kevab + 1 nposi = (locno-1)*ndofn + idofn stfor(nposi) = stfor(nposi) + eload(ielem, kevab) tofor(nposi) = tofor(nposi) + tload(ielem, kevab) End Do End Do End Do Do itotv = 1, ntotv refor = tofor(itotv) - stfor(itotv) resid = resid + refor*refor retot = retot + tofor(itotv)*tofor(itotv) agash = abs(refor) If (agash>remax) remax = agash End Do Do ielem = 1, nelem Do ievab = 1, nevab eload(ielem, ievab) = tload(ielem, ievab) - eload(ielem, ievab) End Do End Do resid = sqrt(resid) retot = sqrt(retot) ratio = 100.0*resid/retot If (ratio>toler) nchek = 1 If (iiter==1) Goto 20 If (ratio>pvalu) nchek = 999 20 pvalu = ratio Write (6, 30) nchek, ratio, remax Write (*, 30) nchek, ratio, remax Return 30 format('0', 3x, 'CONVERGENCE CODE =', i4, 3x, 'ratio=',& E14.6, 3X, ' maximum residual=', E14.6) End Subroutine conver !*************************************************************** !*************************************************************** ! c$debug ! c$large Subroutine dimen(mbufa,melem,mevab,mfron,mmats,mpoin,mstif,& mtotg,mtotv,mvfix,ndofn,nprop,nstre) ! ! c********************************************************************** ! c ! c***This subroutine presets variables associated with dynamic ! c dimensioning. ! c ! c********************************************************************** mbufa=10 melem=40 mfron=80 mmats=5 mpoin=150 mstif=(mfron*mfron-mfron)/2.0+mfron mtotg=melem*9 ndofn=2 mtotv=mpoin*ndofn mvfix=30 nprop=7 mevab=ndofn*9 Return End Subroutine dimen !*************************************************************** !*************************************************************** ! c$debug ! c$large Subroutine flowpl(avect,abeta,dvect,ntype,props,lprop,nstr1,mmats) ! ! c****************************************************************** ! c ! c****This subroutine evaluates the plastic d vector. ! c ! c****************************************************************** Dimension avect(4), dvect(4), props(mmats,7) young=props(lprop,1) poiss=props(lprop,2) hards=props(lprop,6) fmul1=young/(1.0+poiss) If(ntype==1) Goto 60 fmul2=young*poiss*(avect(1)+avect(2)+avect(4))/((1.0+poiss)*(1.0-2.0*poiss)) dvect(1)=fmul1*avect(1)+fmul2 dvect(2)=fmul1*avect(2)+fmul2 dvect(3)=0.5*avect(3)*young/(1.0+poiss) dvect(4)=fmul1*avect(4)+fmul2 Goto 70 60 fmul3=young*poiss*(avect(1)+avect(2))/(1.0-poiss*poiss) dvect(1)=fmul1*avect(1)+fmul3 dvect(2)=fmul1*avect(2)+fmul3 dvect(3)=0.5*avect(3)*young/(1.0+poiss) dvect(4)=fmul1*avect(4)+fmul3 70 denom=hards Do istr1=1, nstr1 denom=denom+avect(istr1)*dvect(istr1) End Do abeta=1.0/denom Return End Subroutine flowpl !*************************************************************** !*************************************************************** ! c$debug ! c$large Subroutine front(asdis,eload,eqrhs,equat,estif,fixed,iffix,iincs,& iiter,gload,gstif,locel,lnods,kresl,mbufa,melem,mevab,mfron,mstif,& mtotv,mvfix,nacva,namev,ndest,ndofn,nelem,nevab,nnode,nofix,npivo,& npoin,ntotv,tdisp,tload,treac,vecrv) ! ! c********************************************************************** ! c ! c***This subroutine undertakes equation solution by the frontal ! c method. ! c ! c********************************************************************** Dimension asdis(mtotv), eload(melem,mevab), eqrhs(mbufa), & equat(mfron,mbufa), estif(mevab,mevab), fixed(mtotv), iffix(mtotv),& npivo(mbufa), vecrv(mfron), gload(mfron), gstif(mstif), lnods(melem,9),& locel(mevab), nacva(mfron), namev(mbufa), ndest(mevab), nofix(mvfix), & noutp(2), tdisp(mtotv), tload(melem,mevab), treac(mvfix,ndofn) nfunc(i,j)=i+(j*j-j)/2 ! c ! c***Change the sign of the last appearance of each node. ! c If(iincs>1 .Or. iiter>1) Goto 455 Do ipoin=1, npoin klast=0 Do ielem=1, nelem Do inode=1, nnode If(lnods(ielem,inode)/=ipoin) Goto 120 klast=ielem nlast=inode 120 End Do End Do If(klast/=0) lnods(klast,nlast)=-ipoin End Do 455 Continue ! c ! c****Start by initializing everything that matters to zero. ! c Do ibufa=1, mbufa eqrhs(ibufa)=0.0 End Do Do istif=1, mstif gstif(istif)=0.0 End Do Do ifron=1, mfron gload(ifron)=0.0 vecrv(ifron)=0.0 nacva(ifron)=0 Do ibufa=1, mbufa equat(ifron,ibufa)=0.0 End Do End Do End Do End Do ! c ! c****And prepare for disc reading and writing operations. ! c nbufa=0 If(kresl>1) nbufa=mbufa !******************************************************* !c$? ! Rewind 1 ! Rewind 2 ! Rewind 3 ! Rewind 4 ! Rewind 8 ! 20230227 !******************************************************** ! c ! c***Enter main element assembly - reduction loop. ! c nfron=0 kelva=0 Do ielem=1, nelem If(kresl>1) Goto 400 kevab=0 Read(1) estif Do inode=1, nnode Do idofn=1, ndofn nposi=(inode-1)*ndofn+idofn locno=lnods(ielem,inode) If(locno>0) locel(nposi)=(locno-1)*ndofn+idofn If(locnonfron) nfron=ndest(kevab) 210 End Do Write(8) locel, ndest, nacva, nfron 400 If(kresl>1) Read(8) locel, ndest, nacva, nfron ! c ! c***Assemble element loads. ! c Do ievab=1, nevab idest=ndest(ievab) gload(idest)=gload(idest)+eload(ielem,ievab) ! c ! c****Assemble the element stiffnesses but not in resolution. ! c If(kresl>1) Goto 402 Do jevab=1, ievab jdest=ndest(jevab) ngash=nfunc(idest,jdest) ngish=nfunc(jdest,idest) If(jdest>=idest) gstif(ngash)=gstif(ngash)+estif(ievab,jevab) If(jdest0.0) Goto 235 Write(6,900) nikno, pivot Stop 235 Continue equat(ifron,nbufa)=0.0 ! c ! c***Enquire whether present variable is free or prescribed. ! c If(iffix(nikno)==0) Goto 250 ! c ! c***Deal with a prescribed deflection. ! c Do jfron=1, nfron gload(jfron)=gload(jfron)-fixed(nikno)*equat(jfron,nbufa) End Do Goto 280 ! c ! c***Eliminate a free variable - deal with the right hand side first. ! c 250 Do jfron=1, nfron gload(jfron)=gload(jfron)-equat(jfron,nbufa)*eqrhs(nbufa)/pivot ! c ! c***Now deal with the coefficients in core. ! c If(kresl>1) Goto 418 If(equat(jfron,nbufa)==0.0) Goto 270 nloca=nfunc(0,jfron) cureq=equat(jfron,nbufa) Do lfron=1, jfron ngash=lfron+nloca gstif(ngash)=gstif(ngash)-cureq*equat(lfron,nbufa)/pivot End Do 418 Continue 270 Continue 280 equat(ifron,nbufa)=pivot ! c ! c***Record the new vacant space, and reduce frontwidth if possible. ! c nacva(ifron)=0 Goto 290 ! c ! c***Complete the element loop in the forward elemination. ! c 300 End Do 290 If(nacva(nfron)/=0) Goto 310 nfron=nfron-1 If(nfron>0) Goto 290 310 End Do End Do If(kresl==1) Write(2) equat, eqrhs, npivo, namev !c ? Backspace 2 !c !c***Enter back - subsroutine phase, loop backwords through variables. !c Do ielva=1, kelva ! c ! c***Read a new block of equations - if needed. ! c If(nbufa/=0) Goto 412 Backspace 2 Read(2) equat, eqrhs, npivo, namev Backspace 2 nbufa=mbufa If(kresl==1) Goto 412 Backspace 4 Read(4) eqrhs Backspace 4 412 Continue ! c ! c***Preare to back - subroutine from the current eguation. ! c ifron=npivo(nbufa) nikno=namev(nbufa) pivot=equat(ifron,nbufa) If(iffix(nikno)/=0) vecrv(ifron)=fixed(nikno) If(iffix(nikno)==0) equat(ifron,nbufa)=0.0 ! c ! c***Back - subroutine in the current equation. ! c Do jfron=1, mfron eqrhs(nbufa)=eqrhs(nbufa)-vecrv(jfron)*equat(jfron,nbufa) End Do ! c ! c***Put the final values where they belong. ! c If(iffix(nikno)==0) vecrv(ifron)=eqrhs(nbufa)/pivot If(iffix(nikno)/=0) fixed(nikno)=-eqrhs(nbufa) nbufa=nbufa-1 asdis(nikno)=vecrv(ifron) ! c ! c***Add displacements to previous total values. ! c Do itotv=1, ntotv tdisp(itotv)=tdisp(itotv)+asdis(itotv) End Do ! c ! c***Store reactions for printing later. ! c kboun=1 Do ipoin=1, npoin nloca=(ipoin-1)*ndofn Do idofn=1, ndofn ngush=nloca+idofn If(iffix(ngush)>0) Goto 360 End Do Goto 370 360 Do idofn=1, ndofn ngash=nloca+idofn treac(kboun,idofn)=treac(kboun,idofn)+fixed(ngash) End Do kboun=kboun+1 370 End Do ! c ! c***Add reactions into the total load array. ! c Do ipoin=1, npoin Do ielem=1, nelem Do inode=1, nnode nloca=iabs(lnods(ielem,inode)) If(ipoin==nloca) Goto 720 End Do 720 Do idofn=1, ndofn ngash=(inode-1)*ndofn+idofn mgash=(ipoin-1)*ndofn+idofn tload(ielem,ngash)=tload(ielem,ngash)+fixed(mgash) End Do Return 900 Format('0',3X,'NEGATIVE OR ZERO PIVOT ENCOUNTERED FOR VARIABLE NO.',& I4,'OF VALUE',E17.6) End Do End Do End Do End Do End Subroutine front !*************************************************************** !*************************************************************** ! c$debug ! c$large Subroutine increm(eload,fixed,iincs,melem,mevab,miter,mtotv,mvfix,& ndofn,nelem,nevab,noutp,nofix,ntotv,nvfix,presc,rload,tfact,tload,toler) ! !c******************************************************************** !c !c***This subroutine increments the applied loading. !c !c******************************************************************** Dimension eload(melem,mevab), fixed(mtotv), iffix(mtotv), noutp(2),& nofix(mvfix), presc(mvfix,ndofn), rload(melem,mevab),& tload(melem,mevab) Write(6,900) iincs Write(*,900) iincs Write(*,*) ' FACTO, TOLER,MITER,NOUTP(1),NOUTP(2)' !c$$$ Read(5,950) facto, toler, miter, noutp(1), noutp(2) ! c write(5, 950) facto, toler, miter, noutp(1), noutp(2) tfact=tfact+facto Write(6,960) tfact, toler, miter, noutp(1), noutp(2) Write(*,960) tfact, toler, miter, noutp(1), noutp(2) Do ielem=1, nelem Do ievab=1, nevab eload(ielem,ievab)=eload(ielem,ievab)+rload(ielem,ievab)*facto tload(ielem,ievab)=tload(ielem,ievab)+rload(ielem,ievab)*facto End Do End Do ! c ! c***Interpret fixity data in vector form. ! c Do itotv=1, ntotv fixed(itotv)=0.0 End Do Do ivfix=1, nvfix nloca=(nofix(ivfix)-1)*ndofn Do idofn=1, ndofn ngash=nloca+idofn fixed(ngash)=presc(ivfix,idofn)*facto End Do End Do Return 960 Format('0',5X,'LOAD FACTOR =',F10.5,5X,'CONVERGENCE TOLERANCE=',& F10.5,//,5X,'MAX. NO. OF ITERATIONS =',I5,//'INITIAL OUTPUT PARAMETER =',& I5,5X,'FINAL OUTPUT PARAMETER=',I5) 900 Format('0',5X,'INCREMENT NUMBER ',I5) 950 Format(2F10.5,3I5) End Subroutine increm !*************************************************************** !*************************************************************** ! c$debug ! c$large Subroutine input(coord,iffix,lnods,matno,melem,mevab,mfron,mmats,& mpoin,mtotv,mvfix,nalgo,ncrit,ndfro,ndofn,nelem,nevab,ngaus,ngau2,& nincs,nmats,nnode,nofix,npoin,nprop,nstre,nstr1,ntotg,ntotv,ntype,& nvfix,posgp,presc,props,weigp) !c !c********************************************************************* !c !c****This subroutine accepts most of the input data. !c !c********************************************************************* Dimension coord(mpoin,2), iffix(mtotv), lnods(melem,9), matno(melem),& ndfro(melem), nofix(mvfix), posgp(4), presc(mvfix,ndofn),& props(mmats,nprop), title(80), weigp(4) !************************************************ !c$?? !************************************************ ! 20230224 !************************************************ ! Rewind 1 ! Rewind 2 ! Rewind 3 ! Rewind 4 ! Rewind 8 !************************************************ Write(*,*) 'TITLE' !c$ Read(5,920) title !c Write(5, 920) title Write(6,920) title Write(*,920) title !c !c***Read the first data card, and echo it immediately. !c Write(*,*) nalgo, ncrit, nincs, nstre !!! Read(5,900) npoin, nelem, nvfix, ntype, nnode, nmats,& ngaus, nalgo, ncrit, nincs, nstre ! c write(5, 900) npoin, nelem, nvfix, ntype, nnode, nmats, ngaus, ! c nalgo, ncrit, nincs, nstre nevab=ndofn*nnode nstr1=nstre+1 If(ntype==3) nstr1=nstre ntotv=npoin*ndofn ngau2=ngaus*ngaus ntotg=nelem*ngau2 Write(6,901) npoin, nelem, nvfix, ntype, nnode, nmats,& ngaus, nevab, nalgo, ncrit, nincs, nstre !c$$$ !c$$$ Write(*,901) npoin, nelem, nvfix, ntype, nnode, nmats,& ngaus, nevab, nalgo, ncrit, nincs, nstre Call check1(ndofn,nelem,ngaus,nmats,nnode,npoin,nstre,& ntype,nvfix,ncrit,nalgo,nincs) ! c ! c***Read the element nodal connections, and the property numbers. ! c Write(6,902) Write(*,902) Do ielem=1, nelem Write(*,*) 'NUMEL, MATNO(NUMEL),(LNODS(NUMEL,INODE),INODE=1,NNODE)' !!! Read(5,903) numel, matno(numel), (lnods(numel,inode),inode=1,nnode) ! c Write(5, 903) numel, matno(numel), (lnods(numel, inode), inode = 1, nnode) Write(6,903) numel, matno(numel), (lnods(numel,inode),inode=1,nnode) Write(*,903) numel, matno(numel), (lnods(numel,inode),inode=1,nnode) End Do ! c ! c***Zero all the nodal coordinates, prior to reading some of them. ! c Do ipoin=1, npoin Do idime=1, 2 coord(ipoin,idime)=0.0 End Do End Do End Do End Do ! c ! c***Read some nodal coordinates, finishing with the node of all. ! c Write(6,904) Write(*,904) 6 Write(*,*) 'IPOIN,(COORD(IPOIN,IDIME),IDIME=1,2)' !!! Read(5,905) ipoin, (coord(ipoin,idime),idime=1,2) ! c Write(5, 905) ipoin, (coord(ipoin, idime), idime = 1, 2). If(ipoin/=npoin) Goto 6 ! c ! c***Interpolate coordinates of mid - side nodes. ! c Call nodexy(coord,lnods,melem,mpoin,nelem,nnode) Do ipoin=1, npoin Write(6,906) ipoin, (coord(ipoin,idime),idime=1,2) Write(*,906) ipoin, (coord(ipoin,idime),idime=1,2) End Do ! c ! c***Read the fixed values. ! c Write(6,907) Write(*,907) Do ivfix=1, nvfix Write(*,*) 'NOFIX(IVFIX),IFPRE,(PRESC(IVFIX,IDOFN),IDOFN=1,NDOFN)' !!! Read(5,908) nofix(ivfix), ifpre, (presc(ivfix,idofn),idofn=1,ndofn) ! c write(5, 908) nofix(ivfix), ifpre, (presc(ivfix, idofn), idofn = 1, ndofn) Write(6,908) nofix(ivfix), ifpre, (presc(ivfix,idofn),idofn=1,ndofn) Write(*,908) nofix(ivfix), ifpre, (presc(ivfix,idofn),idofn=1,ndofn) nloca=(nofix(ivfix)-1)*ndofn ifdof=10**(ndofn-1) Do idofn=1, ndofn ngash=nloca+idofn If(ifpre1.0) sint3=1.0 Goto 20 10 sint3=0.0 20 Continue If(sint31.0) sint3=1.0 theta=asin(sint3)/3.0 Goto(1,2,3,4) ncrit ! c***Tresca. 1 yield=2.0*cos(theta)*steff Return ! c***Von mises. 2 yield=root3*steff Return ! c***Mohr - coulomb. 3 phira=props(lprop,7)*0.017453292 snphi=sin(phira) yield=smean*snphi+steff*(cos(theta)-sin(theta)*snphi/root3) Return ! c***Drucker - prager. 4 phira=props(lprop,7)*0.017453292 snphi=sin(phira) yield=6.0*smean*snphi/(root3*(3.0-snphi))+steff ! (see (eqn.7.18)) Return End Subroutine invar !************************************************ !************************************************ ! c$debug ! c$large Subroutine loadps(coord,lnods,matno,melem,mmats,mpoin,nelem,& nevab,ngaus,nnode,npoin,nstre,ntype,posgp,props,rload,weigp,ndofn) ! c******************************************************************** ! c ! c****This subroutine evaluates the consistent nodal forces for each ! c element. ! c ! c******************************************************************** Dimension cartd(2,9), coord(mpoin,2), deriv(2,9), dgash(2),& dmatx(4,4), elcod(2,9), lnods(melem,9), matno(melem),& noprs(4), pgash(2), point(2), posgp(4), press(4,2),& props(mmats,7), rload(melem,18), shape(9), stran(4),& stres(4), title(12), weigp(4), gpcod(2,9) twopi=6.283185308 Do ielem=1, nelem Do ievab=1, nevab rload(ielem,ievab)=0.0 End Do End Do Write(*,*) ' TITLE' !!! Read(5,901) title ! c write(5, 901) title Write(6,903) title ! c ! c***Read data controlling loading types to be inputted. ! c Write(*,*) ' IPLOD,IGRAV,IEDGE' !!! Read(5,918) iplod, igrav, iedge ! c write(5, 918) iplod, igrav, iedge Write(6,919) iplod, igrav, iedge Write(*,919) iplod, igrav, iedge ! c ! c***Read nodal point loads. ! c If(iplod==0) Goto 500 20 Write(*,*) 'LODPT,(POINT(IDOFN),IDOFN=1,2)' !!! Read(5,931) lodpt, (point(idofn),idofn=1,2) ! c write(5, 931) lodpt, (point(idofn), idofn=1, 2) Write(6,931) lodpt, (point(idofn),idofn=1,2) Write(*,931) lodpt, (point(idofn),idofn=1,2) ! c ! c***Associate the nodal point loads with an element. ! c Do ielem=1, nelem Do inode=1, nnode nloca=iabs(lnods(ielem,inode)) If(lodpt==nloca) Goto 40 End Do End Do 40 Do idofn=1, 2 ngash=(inode-1)*2+idofn rload(ielem,ngash)=point(idofn) End Do If(lodptncode) ngash=1 If(knode>ncode) mgash=2 rload(neass,ngash)=rload(neass,ngash)+shape(kount)*pxcom*dvolu rload(neass,mgash)=rload(neass,mgash)+shape(kount)*pycom*dvolu End Do End Do End Do 700 Continue Write(6,907) Write(*,907) Do ielem=1, nelem Write(6,905) ielem, (rload(ielem,ievab),ievab=1,nevab) Write(*,905) ielem, (rload(ielem,ievab),ievab=1,nevab) End Do Return 901 Format(2X,80A) 903 Format(2X,80A) 918 Format(3I5) 919 Format(1X,'IPLOD=',I5,5X,'IGRAV=',I5,5X,'IEDGE=',I5) 931 Format(I5,2F10.3) 906 Format(2F10.3) 911 Format('0','GRAVITY ANGLE =',F10.3,' GRAVITY CONSTANT =',F10.3) 932 Format(I5) 912 Format('0',5X,'NO. OF LOADED EDGES =',I5) 915 Format('0',5X,'LIST OF LOADED EDGES AND APPLIED LOADS') 902 Format(4I5) 913 Format(I10,5X,3I5) 914 Format(6F10.3) 907 Format('0',5X,' TOTAL NODAL FORCES FOR EACH ELEMENT') 905 Format(1X,I3,2X,8F9.2/(5X,8F9.2)) End Do End Subroutine loadps !************************************************ !c$ ! 20230228 !************************************************ ! c$debug ! c$large Subroutine output(iiter,mtotg,mtotv,mvfix,nelem,ngaus,nofix,& noutp,npoin,nvfix,strsg,tdisp,treac,epstn,ntype,nchek) ! ! c******************************************************************* ! c ! c****This subroutine output displacements reaction and stresses. ! c ! c******************************************************************* Dimension nofix(mvfix), noutp(2), strsg(4,mtotg), strsp(3),& tdisp(mtotv), treac(mvfix,2), epstn(mtotg) koutp=noutp(1) If(iiter>1) koutp=noutp(2) If(iiter==1 .And. nchek==0) koutp=noutp(2) ! c ! c***Output displacements. ! c If(koutp

内容很感兴趣，尝试运行实践程序，比较困难。

对学习有限元非常有帮助

下载学习

我只能说，灰常不错的一本书啊！！！！

下下来看看，谢谢楼主！