m6law.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/.
!||====================================================================
!||    m6law           ../engine/source/materials/mat/mat006/m6law.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 m6law(
     1   PM,      OFF,     SIG,     EINT,
     2   RHO,     RK,      RE,      VORTI,
     3   VNEW,    RHO0,    VIS,     PSH,
     4   D1,      D2,      D3,      D4,
     5   D5,      D6,
     6   PC,      AMU,     AMU2,    ESPE,
     7   SSP,     DF,      TMU,     MAT,
     8   WXX,     WYY,     WZZ,     NEL,
     9   DPDM,    JTUR,    JPOR,    EOS_STRUCT)
C-----------------------------------------------
C   M o d u l e s
C-----------------------------------------------
        USE eos_param_mod , ONLY : eos_param_
C-----------------------------------------------
#include      "implicit_f.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) :: JTUR
      INTEGER, INTENT(IN) :: JPOR
      INTEGER MAT(*),NEL
      my_real
     .   PM(NPROPM,*), OFF(*), SIG(NEL,6), EINT(*), RHO(*), RK(*), RE(*),
     .   VORTI(*)
      my_real
     .   VNEW(*), RHO0(*), D1(*), D2(*), D3(*), D4(*), D5(*), D6(*),
     .   PC(*), VIS(*),
     .   amu(*), amu2(*), espe(*), ssp(*), df(*),
     .   tmu(*), psh(*), wxx(*), wyy(*), wzz(*), rho0_1,
     .   pc_1, tmu_1, psh_1
      my_real, INTENT(IN) :: dpdm(nel)
      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
     .   DAV, RK2T, VIS2, FAC
C-----------------------------------------------
 
      !SAUVEGARDE VORTICITE SI ANIM
      IF((anim_e(10)==1 .OR. anim_se(10)==1).AND. dt1/=zero)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
 
      !EOS COEFFICIENTS
      mx     =mat(1)
      rho0_1=pm( 1,mx)
      pc_1  =pm(37,mx)
      tmu_1 =pm(81,mx)
      psh_1=pm(88,mx)
      DO 10 i=1,nel
      rho0(i)=rho0_1
      vis(i) =pm(24,mx)*rho(i)
      pc(i)  =pc_1
      tmu(i) =tmu_1
      psh(i)=psh_1
   10 CONTINUE
 
      !IF TURBULENT
      IF(jtur/=0)THEN
        IF(jpor/=2)THEN
          DO i=1,nel
            tmu(i)=tmu(i)*rk(i)*rk(i) / max(em15,(re(i)*vnew(i)))
            vis(i)=vis(i)+tmu(i)
          ENDDO
        ELSE
          DO i=1,nel
            tmu(i)=zero
          ENDDO
        ENDIF
      ELSE
        DO i=1,nel
          rk(i)=zero
        ENDDO
      ENDIF
 
      !OUTPUT
      DO i=1,nel
       vis2 = two*vis(i)
       dav  = -(d1(i)+d2(i)+d3(i))*third
       !------------------------------!
       !   VISCOUS STRESS TENSOR      !
       !------------------------------!
       sig(i,1)=vis2*(d1(i)+dav)
       sig(i,2)=vis2*(d2(i)+dav)
       sig(i,3)=vis2*(d3(i)+dav)
       sig(i,4)=vis(i) *d4(i)
       sig(i,5)=vis(i) *d5(i)
       sig(i,6)=vis(i) *d6(i)
       amu2(i)= max(zero,amu(i))**2
       espe(i)=df(i)*eint(i) / max(em15,vnew(i))
       rk2t = two*rk(i) / max(em15,three*vnew(i))
       
       !------------------------------!
       !         SOUND SPEED          !
       !------------------------------!
       ssp(i)=sqrt(abs(dpdm(i))/rho0(i))
      ENDDO
 
      RETURN
      END