C SIGMAPR3 SOURCE CB215821 16/04/21 21:18:26 8920 subroutine sigmapr3ETC(sigma,S1,S2,S3,rcossigmapr,lcp) c ****************************************************************** c * * c * Principal stresses ( Batoz, Vol. 1., 1.2.5. p 51) * c * ================================================= * c * * c ****************************************************************** c * sigma : vector of the stresses in the global axis c : Sxx, Syy, Szz, Sxy, Sxz, Syz c * S1 : largest principal stress c * S2 : intermediate principal stress c * S3 : smallest principal stress c * rcossigmapr : ( lx ly lz ) c ( mx my mz ) c ( nx ny nz ) c for S1 S2 S3 IMPLICIT REAL*8 (A-B,D-H,O-Z) implicit integer (I-K,M,N) implicit logical (L) implicit character*10 (C) character*20 cloc1,cloc2 dimension sigma(6) dimension rcossigmapr(3,3) i1 = 1 i2 = 2 i3 = 3 i4 = 4 i5 = 5 i6 = 6 i7 = 7 i8 = 8 i9 = 9 i10=10 i11=11 i12=12 i13=13 i14=14 i15=15 i16=16 i17=17 i18=18 i19=19 i20=20 i21=21 i22=22 i23=23 i24=24 i25=25 i26=26 i27=27 i28=28 i29=29 i30=30 i31=31 i32=32 i33=33 i34=34 i35=35 i36=36 i37=37 i38=38 i39=39 r0 = 0. r1 = 1. r2 = 2. r3 = 3. r6 = 6. precision = 1.0d-3 precision2 = precision*precision precision3 = precision2*precision pi = ACOS(-r1) racine2 = SQRT(r2) S1 = r0 S2 = r0 S3 = r0 do jloc=i1,i3 do iloc=i1,i3 RCOSSIGMAPR(iloc,jloc) = r0 end do end do Sxx = sigma(i1) Syy = sigma(i2) Szz = sigma(i3) Sxy = sigma(i4) Sxz = sigma(i5) Syz = sigma(i6) if (lcp) then c PLANE STRESS sxx = sigma(1) syy = sigma(2) sxy = sigma(4) diff = sxx-syy rloc = diff*diff+4.*sxy*sxy cloc1 = 'rloc' cloc2 = 'sigmapr3ETCtg' call check_sqrt(rloc,cloc1,cloc2,lerror) c! **** ********** if (lerror) return S1 = (SQRT(rloc) + sxx+syy) / r2 S2 = (-SQRT(rloc) + sxx+syy) / r2 S3 = 0.d0 if (S1.lt.S2) then rloc = S1 S1 = S2 S2 = rloc endif lfirststress = .true. Smax = S1 1 continue if (.not.lfirststress) Smax = S2 c Main direction of Smax c ----------------- c condition for the case Smax = Sxx = Syy (then deltax=deltay=deltaz=0) if ( (ABS(Sxx-Smax).le.r1).and.(ABS(Syy-Smax).le.r1).and. . (ABS(Szz-Smax).gt.r1) ) then if(lfirststress) then rnx = (SQRT(r2))/r2 else rnx = - (SQRT(r2))/r2 endif rny = (SQRT(r2))/r2 rnz = r0 c condition for the case Smax = Syy = Szz (then deltax=deltay=deltaz=0) else if ( (ABS(Syy-Smax).le.r1).and.(ABS(Szz-Smax).le.r1) . .and.(ABS(Sxx-Smax).gt.r1) ) then if(lfirststress) then rny = (SQRT(r2))/r2 else rny = - (SQRT(r2))/r2 endif rnz = (SQRT(r2))/r2 rnx = r0 c condition for the case Smax = Sxx = Szz (then deltax=deltay=deltaz=0) else if ( (ABS(Sxx-Smax).le.r1).and.(ABS(Szz-Smax).le.r1) . .and.(ABS(Syy-Smax).gt.r1) ) then if(lfirststress) then rnx = (SQRT(r2))/r2 else rnx = - (SQRT(r2))/r2 endif rnz = (SQRT(r2))/r2 rny = r0 c condition for the case Smax = Sxx = Syy = Szz (then deltax=deltay=deltaz=0) else if ( (ABS(Sxx-Smax).le.r1).and.(ABS(Syy-Smax).le.r1) . .and.(ABS(Szz-Smax).le.r1).and.(ABS(Sxy).le.r1).and. . (ABS(Sxz).le.r1).and.(ABS(Syz).le.r1) ) then rnx = (SQRT(r3))/r3 rny = (SQRT(r3))/r3 rnz = (SQRT(r3))/r3 else deltax = (Szz-Smax)*(Syy-Smax)-Syz*Syz deltay = (Sxx-Smax)*(Szz-Smax)-Sxz*Sxz deltaz = (Sxx-Smax)*(Syy-Smax)-Sxy*Sxy deltamax = ABS(deltax) ldeltay = .false. ldeltaz = .false. if (ABS(deltay).gt.deltamax) then deltamax = ABS(deltay) ldeltay = .true. endif if (ABS(deltaz).gt.deltamax) then deltamax = ABS(deltaz) ldeltay = .false. ldeltaz = .true. endif if (ldeltaz) then c Main orientation OZ aa = ( Sxy*Syz - Sxz*(Syy-Smax) ) / deltaz bb = ( Sxz*Sxy - Syz*(Sxx-Smax) ) / deltaz rloc = SQRT( aa*aa + bb*bb + r1 ) rnz = r1/rloc rnx = aa/rloc rny = bb/rloc else if (ldeltay) then c Main orientation OY aa = ( Sxz*Syz - Sxy*(Szz-Smax) ) / deltay bb = ( Sxz*Sxy - Syz*(Sxx-Smax) ) / deltay rloc = SQRT( aa*aa + bb*bb + r1 ) rny = r1/rloc rnx = aa/rloc rnz = bb/rloc else c Main orientation OX aa = ( Sxz*Syz - Sxy*(Szz-Smax) ) / deltax bb = ( Syz*Sxy - Sxz*(Syy-Smax) ) / deltax rloc = SQRT( aa*aa + bb*bb + r1 ) rnx = r1/rloc rny = aa/rloc rnz = bb/rloc endif endif c The 3 following lines to try to solve the problem of unsymetry c in routine STRESS if (ABS(rnx).lt.precision3) rnx = r0 if (ABS(rny).lt.precision3) rny = r0 if (ABS(rnz).lt.precision3) rnz = r0 if(lfirststress) then c go and calculate for S2 c ----------------------- lfirststress = .false. rcossigmapr(i1,i1) = rnx rcossigmapr(i2,i1) = rny rcossigmapr(i3,i1) = rnz go to 1 else c S2 has been calculated c ---------------------- rcossigmapr(i1,i2) = rnx rcossigmapr(i2,i2) = rny rcossigmapr(i3,i2) = rnz endif c S1 and S2 have been calculated c ------------------------------ c calculate S3. c ------------- rcossigmapr(i1,i3) = rcossigmapr(i2,i1)*rcossigmapr(i3,i2) . - rcossigmapr(i3,i1)*rcossigmapr(i2,i2) rcossigmapr(i2,i3) = rcossigmapr(i3,i1)*rcossigmapr(i1,i2) . - rcossigmapr(i1,i1)*rcossigmapr(i3,i2) rcossigmapr(i3,i3) = rcossigmapr(i1,i1)*rcossigmapr(i2,i2) . - rcossigmapr(i2,i1)*rcossigmapr(i1,i2) else c FULL 3D c Hydrostatic stress c ------------------ rI1 = Sxx+Syy+Szz Sh = rI1/r3 SumS = rI1+Sxy+Sxz+Syz c 2nd invariant of the deviator J c ------------------------------- rJ2 = ( (Sxx-Syy) * (Sxx-Syy) . + (Sxx-Szz) * (Sxx-Szz) . + (Syy-Szz) * (Syy-Szz) ) / r6 . + Sxy*Sxy + Sxz*Sxz + Syz*Syz c octahedral shear stress c ----------------------- toc = r2*rJ2/r3 cloc1 = 'toc' cloc2 = 'Sigmapr3' toc = SQRT(toc) smallS = SumS*precision2 if (ABS(toc).gt.ABS(smallS)) then c 3th invariant of the deviator J c ------------------------------- rI2 = Sxx*Syy + Sxx*Szz + Syy*Szz . - Sxy*Sxy - Sxz*Sxz - Syz*Syz rI3 = Sxx*Syy*Szz + r2*Sxy*Syz*Sxz . - Sxx*Syz*Syz - Syy*Sxz*Sxz - Szz*Sxy*Sxy rJ3 = rI3 - rI1*rI2/r3 + r2*rI1*rI1*rI1/27d0 c angle omega c ----------- rloc = ( racine2*rJ3 / (toc*toc*toc) ) if (rloc.ge.1d0) then omega = r0 else if (rloc.le.(-1d0)) then omega = ( ACOS( (-1d0) ) )/r3 else omega = ( ACOS( racine2*rJ3 / (toc*toc*toc) ) )/r3 if ((omega.lt.r0).or.(omega.gt.pi/r3)) then print *,'Subroutine sigmapr3 - ERROR in omega' print *,'omega = ',omega stop endif endif c principal stresses c ------------------ S1 = racine2*toc*COS(omega-r2*pi/r3)+Sh S2 = racine2*toc*COS(omega+r2*pi/r3)+Sh S3 = racine2*toc*COS(omega) +Sh c sort the principal stresses : S1 > S2 > S3 if (S2.lt.S3) then rloc = S2 S2 = S3 S3 = rloc endif if (S1.lt.S2) then rloc = S1 S1 = S2 S2 = rloc endif if (S2.lt.S3) then rloc = S2 S2 = S3 S3 = rloc endif lfirststress = .true. Smax = S1 2 continue if (.not.lfirststress) Smax = S2 c Main direction of Smax c ----------------- c condition for the case Smax = Sxx = Syy (then deltax=deltay=deltaz=0) if ( (ABS(Sxx-Smax).le.r1).and.(ABS(Syy-Smax).le.r1).and. . (ABS(Szz-Smax).gt.r1) ) then if(lfirststress) then rnx = (SQRT(r2))/r2 else rnx = - (SQRT(r2))/r2 endif rny = (SQRT(r2))/r2 rnz = r0 c condition for the case Smax = Syy = Szz (then deltax=deltay=deltaz=0) else if ( (ABS(Syy-Smax).le.r1).and.(ABS(Szz-Smax).le.r1) . .and.(ABS(Sxx-Smax).gt.r1) ) then if(lfirststress) then rny = (SQRT(r2))/r2 else rny = - (SQRT(r2))/r2 endif rnz = (SQRT(r2))/r2 rnx = r0 c condition for the case Smax = Sxx = Szz (then deltax=deltay=deltaz=0) else if ( (ABS(Sxx-Smax).le.r1).and.(ABS(Szz-Smax).le.r1) . .and.(ABS(Syy-Smax).gt.r1) ) then if(lfirststress) then rnx = (SQRT(r2))/r2 else rnx = - (SQRT(r2))/r2 endif rnz = (SQRT(r2))/r2 rny = r0 c condition for the case Smax = Sxx = Syy = Szz (then deltax=deltay=deltaz=0) else if ( (ABS(Sxx-Smax).le.r1).and.(ABS(Syy-Smax).le.r1) . .and.(ABS(Szz-Smax).le.r1).and.(ABS(Sxy).le.r1).and. . (ABS(Sxz).le.r1).and.(ABS(Syz).le.r1) ) then rnx = (SQRT(r3))/r3 rny = (SQRT(r3))/r3 rnz = (SQRT(r3))/r3 else deltax = (Szz-Smax)*(Syy-Smax)-Syz*Syz deltay = (Sxx-Smax)*(Szz-Smax)-Sxz*Sxz deltaz = (Sxx-Smax)*(Syy-Smax)-Sxy*Sxy deltamax = ABS(deltax) ldeltay = .false. ldeltaz = .false. if (ABS(deltay).gt.deltamax) then deltamax = ABS(deltay) ldeltay = .true. endif if (ABS(deltaz).gt.deltamax) then deltamax = ABS(deltaz) ldeltay = .false. ldeltaz = .true. endif if (ldeltaz) then c Main orientation OZ aa = ( Sxy*Syz - Sxz*(Syy-Smax) ) / deltaz bb = ( Sxz*Sxy - Syz*(Sxx-Smax) ) / deltaz rloc = SQRT( aa*aa + bb*bb + r1 ) rnz = r1/rloc rnx = aa/rloc rny = bb/rloc else if (ldeltay) then c Main orientation OY aa = ( Sxz*Syz - Sxy*(Szz-Smax) ) / deltay bb = ( Sxz*Sxy - Syz*(Sxx-Smax) ) / deltay rloc = SQRT( aa*aa + bb*bb + r1 ) rny = r1/rloc rnx = aa/rloc rnz = bb/rloc else c Main orientation OX aa = ( Sxz*Syz - Sxy*(Szz-Smax) ) / deltax bb = ( Syz*Sxy - Sxz*(Syy-Smax) ) / deltax rloc = SQRT( aa*aa + bb*bb + r1 ) rnx = r1/rloc rny = aa/rloc rnz = bb/rloc endif endif c The 3 following lines to try to solve the problem of unsymetry c in routine STRESS if (ABS(rnx).lt.precision3) rnx = r0 if (ABS(rny).lt.precision3) rny = r0 if (ABS(rnz).lt.precision3) rnz = r0 if(lfirststress) then c go and calculate for S2 c ----------------------- lfirststress = .false. rcossigmapr(i1,i1) = rnx rcossigmapr(i2,i1) = rny rcossigmapr(i3,i1) = rnz go to 2 else c S2 has been calculated c ---------------------- rcossigmapr(i1,i2) = rnx rcossigmapr(i2,i2) = rny rcossigmapr(i3,i2) = rnz endif c S1 and S2 have been calculated c ------------------------------ c calculate S3. c ------------- rcossigmapr(i1,i3) = rcossigmapr(i2,i1)*rcossigmapr(i3,i2) . - rcossigmapr(i3,i1)*rcossigmapr(i2,i2) rcossigmapr(i2,i3) = rcossigmapr(i3,i1)*rcossigmapr(i1,i2) . - rcossigmapr(i1,i1)*rcossigmapr(i3,i2) rcossigmapr(i3,i3) = rcossigmapr(i1,i1)*rcossigmapr(i2,i2) . - rcossigmapr(i2,i1)*rcossigmapr(i1,i2) else S1 = sigma(i1) S2 = sigma(i2) S3 = sigma(i3) c case where S1 = S2 = S3 rcossigmapr(1,1)=1.0d0 rcossigmapr(2,2)=1.0d0 rcossigmapr(3,3)=1.0d0 endif endif return end