1 : C DDSIRANK SOURCE AM 12/02/28 21:15:09 7291 2 : SUBROUTINE ddsigRank(sigma,d2gt,i1,i2,i3,i4,i5,i6,lcp, 3 : . parahot3,idimpara3) 4 : 5 : c This subroutine calculates the second derivative of the Rankine 6 : c plastic function with respect to sigma 7 : c gt = ft = S1 = largest principal stress 8 : 9 : IMPLICIT REAL*8 (A-B,D-H,O-Z) 10 : implicit integer (I-K,M,N) 11 : implicit logical (L) 12 : implicit character*10 (C) 13 : 14 : **** dimension sigma(i6),dgt(i6),d2gt(i6,i6) 15 : dimension sigma(i6),dgt(6),d2gt(i6,i6) 16 : **** dimension rcossigmapr(i3,i3),P1(3,3),pi1(3) 17 : dimension rcossigmapr(3,3),P1(3,3),pi1(3) 18 : dimension sigmacp(3),vloc(3),vloc2(1,3),vloc1(1) 19 : dimension vloc3(3,3),d2gtcp(3,3),parahot3(idimpara3) 20 : dimension sigmaim(6),sigmail(6),rhi(6),dgtm(6),dgtl(6) 21 : 22 : 23 : 24 : ***** 25 : * MESSAGES D'ERREUR ( SUPPRESSION DES AUTOMATIC OBJECTS) 26 : IF(I3.GT.3) PRINT *, ' DDSIGRANK - ERREUR I2 = ', I2, ' > 2 ' 27 : IF(I6.GT.6) PRINT *, ' DDSIGRANK - ERREUR I6 = ', I6, ' > 6 ' 28 : ***** 29 : 30 : 31 : if (lcp) then 31 : c In plane stress, the second derivative can be calculated analytically 32 : do jloc=1,6 33 : do iloc=1,6 34 : d2gt(iloc,jloc) = 0.d0 35 : enddo 36 : enddo 37 : P1(1,1) = 0.5d0 38 : P1(2,2) = 0.5d0 39 : P1(1,2) = -0.5d0 40 : P1(2,1) = -0.5d0 41 : do iloc=1,2 42 : P1(iloc,3) = 0.d0 43 : P1(3,iloc) = 0.d0 44 : enddo 45 : P1(3,3) = 2.d0 46 : pi1(1) = 1.d0 47 : pi1(2) = 1.d0 48 : pi1(3) = 0.d0 49 : sigmacp(1) = sigma(1) 50 : sigmacp(2) = sigma(2) 51 : sigmacp(3) = sigma(4) 51 : call mulAB(P1,sigmacp,vloc,3,3,3,1) 51 : call mulATB(sigmacp,vloc,vloc1,3,1,3,1) 52 : rloc = 0.5d0 * vloc1(1) 53 : psi1 = dsqrt(rloc) 53 : call mulATB(sigmacp,P1,vloc2,3,1,3,3) 53 : call mulAB(vloc,vloc2,vloc3,3,1,1,3) 54 : 55 : test = 0.5d0*(sigmacp(1)+sigmacp(2)) 56 : psi1lim=max(1.E-7*test,1.d-3) 57 : if (dabs(psi1).le.psi1lim) then 58 : continue 59 : c d2gt = 0 60 : else 61 : do jloc=1,3 62 : do iloc=1,3 63 : d2gtcp(iloc,jloc)=(P1(iloc,jloc)/(2.0d0*psi1)) 64 : . - (vloc3(iloc,jloc)/(4.0d0*psi1*psi1*psi1)) 65 : enddo 66 : enddo 67 : do iloc=1,2 68 : d2gt(1,iloc) = d2gtcp(1,iloc) 69 : d2gt(2,iloc) = d2gtcp(2,iloc) 70 : d2gt(4,iloc) = d2gtcp(3,iloc) 71 : d2gt(iloc,4) = d2gtcp(iloc,3) 72 : end do 73 : d2gt(4,4) = d2gtcp(3,3) 74 : endif 75 : else 76 : c In full 3D, numerical differentiation 77 : do jloc=1,6 78 : do iloc=1,6 79 : d2gt(iloc,jloc) = 0.d0 80 : enddo 81 : enddo 82 : rh = 1.0d-5 83 : do iloc=1,6 84 : c ith component of sigma => ith column of d2gt 85 : c calculation of sigma + hi * ei and sigma - hi * ei 86 : do jloc=1,6 87 : rhi(jloc) = 0.d0 88 : enddo 89 : rhi(iloc) = 1.0d0 * rh 90 : do jloc=1,6 91 : sigmaim(jloc) = sigma(jloc) + rhi(jloc) 92 : sigmail(jloc) = sigma(jloc) - rhi(jloc) 93 : enddo 94 : c Calculation of q(x+hi*ei) = dft(x+hi*ei)/dsig 94 : call dsigRank(sigmaim,dgtm,3,6,parahot3,idimpara3, 95 : . rkappait,lcp) 95 : call dsigRank(sigmail,dgtl,3,6,parahot3,idimpara3, 96 : . rkappait,lcp) 97 : do jloc=1,6 98 : d2gt(jloc,iloc) = (dgtm(jloc) - dgtl(jloc))/(2.0d0*rh) 99 : enddo 100 : enddo 101 : c Calculation of f(x+hi*ei) 102 : endif 103 : 104 : RETURN 105 : END 106 : 107 :
© Cast3M 2003 - All rights reserved.
Disclaimer