m17law.F

Go to the documentation of this file.

Copyright>        OpenRadioss
Copyright>        Copyright (C) 1986-2025 Altair Engineering Inc.
Copyright>
Copyright>        This program is free software: you can redistribute it and/or modify
Copyright>        it under the terms of the GNU Affero General Public License as published by
Copyright>        the Free Software Foundation, either version 3 of the License, or
Copyright>        (at your option) any later version.
Copyright>
Copyright>        This program is distributed in the hope that it will be useful,
Copyright>        but WITHOUT ANY WARRANTY; without even the implied warranty of
Copyright>        MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
Copyright>        GNU Affero General Public License for more details.
Copyright>
Copyright>        You should have received a copy of the GNU Affero General Public License
Copyright>        along with this program.  If not, see <https://www.gnu.org/licenses/>.
Copyright>
Copyright>
Copyright>        Commercial Alternative: Altair Radioss Software
Copyright>
Copyright>        As an alternative to this open-source version, Altair also offers Altair Radioss
Copyright>        software under a commercial license.  Contact Altair to discuss further if the
Copyright>        commercial version may interest you: https://www.altair.com/radioss/.
!||====================================================================
!||    m17law          ../engine/source/materials/mat/mat017/m17law.F
!||--- called by ------------------------------------------------------
!||    mmain           ../engine/source/materials/mat_share/mmain.f90
!||--- uses       -----------------------------------------------------
!||    eos_param_mod   ../common_source/modules/mat_elem/eos_param_mod.F90
!||====================================================================
      SUBROUTINE m17law(
     1   PM,      OFF,     SIG,     EINT,
     2   RHO,     RK,      RE,      VORTI,
     3   WXX,     WYY,     WZZ,     VOLN,
     4   MAT,     VIS,     D1,      D2,
     5   D3,      D4,      D5,      D6,
     6   NEL,     ALOGEY,  SSP,     RHO0,
     7   TMU,     AMU,     AMU2,    PSH,
     8   PC,      ESPE,    C1,      C2,
     9   C3,      C4,      C5,      C6,
     A   DF,      DPDM,    JPOR,    EOS_STRUCT)
C-----------------------------------------------
C   M o d u l e s
C-----------------------------------------------
        USE eos_param_mod , ONLY : eos_param_
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
C-----------------------------------------------
C   G l o b a l   P a r a m e t e r s
C-----------------------------------------------
#include      "mvsiz_p.inc"
C-----------------------------------------------
C   C o m m o n   B l o c k s
C-----------------------------------------------
#include      "com01_c.inc"
#include      "com08_c.inc"
#include      "param_c.inc"
#include      "scr14_c.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER, INTENT(IN) :: JPOR
      INTEGER MAT(*),NEL
      my_real
     .   PM(NPROPM,*), OFF(*), SIG(NEL,6), EINT(*), RHO(*), RK(*), RE(*),
     .   VORTI(*),WXX(*),WYY(*),WZZ(*),VOLN(MVSIZ),VIS(*),
     .   D1(*), D2(*), D3(*), D4(*), D5(*), D6(*), ALOGEY(*), SSP(*),
     .   RHO0(*), TMU(*), AMU(*), AMU2(*), PSH(*), PC(*), ESPE(*),
     .   C1(*), C2(*), C3(*), C4(*), C5(*), C6(*), DF(*), DPDM(*)
      TYPE(EOS_PARAM_),INTENT(INOUT) :: EOS_STRUCT
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I, MX
      my_real
     .   VIS2(MVSIZ), 
     .   DAV(MVSIZ),
     .   YP0, CMU, AX, E, A, XMU,
     .   XM, XK, XE, YPLUS, RK2T, FAC,
     .   rho0_1, vis_1, pc_1, c1_1, c2_1,
     .   c3_1, c4_1, c5_1, c6_1, psh_1     
C-----------------------------------------------
C
      ! SAUVEGARDE VORTICITE SI ANIM
      IF((anim_e(10)==1 .OR. anim_se(10)==1).AND. dt1/=0.)THEN
C      IF(ANIM_E(10)==1.AND. DT1/=0.)THEN
        fac=four/dt1
        IF(n2d==0)THEN
          DO 5 i=1,nel
 5        vorti(i)=fac*sqrt(wxx(i)**2+wyy(i)**2+wzz(i)**2)
        ELSE
          DO 6 i=1,nel
 6        vorti(i)=fac*wzz(i)
        ENDIF
      ENDIF
C
      mx     =mat(1)
C
      rho0_1=pm( 1,mx)
      vis_1 =pm(24,mx)
      pc_1  =pm(37,mx)
      c1_1  = eos_struct%UPARAM(1)
      c2_1  = eos_struct%UPARAM(2)
      c3_1  = eos_struct%UPARAM(3)
      c4_1  = eos_struct%UPARAM(4)
      c5_1  = eos_struct%UPARAM(5)
      c6_1  = eos_struct%UPARAM(6)
      psh_1 =pm(88,mx)
C
      DO 10 i=1,nel
      rho0(i)=rho0_1
      vis(i) =vis_1
      pc(i)  =pc_1
      c1(i)  =c1_1
      c2(i)  =c2_1
      c3(i)  =c3_1
      c4(i)  =c4_1
      c5(i)  =c5_1
      c6(i)  =c6_1
      psh(i) =psh_1
   10 CONTINUE
C
      mx     =mat(1)
      IF(jpor/=2)THEN
      DO 20 i=1,nel
      yp0=pm(51,mx)
      cmu=pm(81,mx)
      ax =pm(47,mx)
      e  =pm(48,mx)
      a  =pm(49,mx)
      xmu      =rho(i)*vis(i)
      xm       =rho(i)*voln(i)
      xk       =rk(i)/xm
      xe       = max(em15,re(i)/xm)
      yplus    =cmu*xk**2/(ax*xe*vis(i))
      yplus    = max(yplus,yp0)
      alogey(i)= a * log(e*yplus)
      vis(i)   =xmu*ax*yplus/alogey(i)
      tmu(i)   =vis(i)-xmu
   20 CONTINUE
      ELSE
        DO 21 i=1,nel
         vis(i)=zero
         tmu(i)=zero 
   21   CONTINUE
      ENDIF
C
      DO 30 i=1,nel
      vis2(i)=two*vis(i)
   30 dav(i) =-(d1(i)+d2(i)+d3(i))*third
C------------------------------
C     DEVIATORIC STRESSES
C------------------------------
      DO 40 i=1,nel
      sig(i,1)=vis2(i)*(d1(i)+dav(i))
      sig(i,2)=vis2(i)*(d2(i)+dav(i))
      sig(i,3)=vis2(i)*(d3(i)+dav(i))
      sig(i,4)=vis(i) *d4(i)
      sig(i,5)=vis(i) *d5(i)
   40 sig(i,6)=vis(i) *d6(i)
C
      DO 50 i=1,nel
      df(i)  =rho0(i)/rho(i)
      amu(i) =one/df(i)-one
      amu2(i)= max(zero,amu(i))**2
   50 espe(i)=df(i)*eint(i)/voln(i)
C
      DO 60 i=1,nel
      rk2t=two*rk(i)/(three*voln(i))
      dpdm(i) = dpdm(i)
     .     +(c5(i)+c6(i)*amu(i))*df(i)*df(i)*rk2t + rk2t*df(i)
   60 CONTINUE
C
      DO 70 i=1,nel
   70 ssp(i)=sqrt(abs(dpdm(i))/rho0(i))
C
      RETURN
      END