sboltlaw.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/.
!||====================================================================
!||    sboltlaw         ../engine/source/elements/solid/solide/sboltlaw.F
!||--- called by ------------------------------------------------------
!||    mmain            ../engine/source/materials/mat_share/mmain.F90
!||--- calls      -----------------------------------------------------
!||    mdtsph           ../engine/source/materials/mat_share/mdtsph.F
!||    mqviscb          ../engine/source/materials/mat_share/mqviscb.F
!||--- uses       -----------------------------------------------------
!||    glob_therm_mod   ../common_source/modules/mat_elem/glob_therm_mod.f90
!||====================================================================
      SUBROUTINE sboltlaw(
     1   PM,      SIG,     MAT,     D1,
     2   D2,      D3,      D4,      D5,
     3   D6,      NEL,     RHO,     BPRELD,
     4   EINT,    QOLD,    VOL0,    STIFN,
     5   DT2T,    NELTST,  ITYPTST, OFFG,
     6   GEO,     PID,     MUMAX,   NGL,
     7   SSP,     DVOL,    AIRE,    VNEW,
     8   VD2,     DELTAX,  VIS,     PNEW,
     9   PSH,     QNEW,    SSP_EQ,  SOLD1,
     A   SOLD2,   SOLD3,   SOLD4,   SOLD5,
     B   SOLD6,   MSSA,    DMELS,   CONDE,
     C   AMU,     VOL_AVG, DTEL,    G_DT,
     D   OFF,     IPM,     RHOREF,  RHOSP,
     E   VOL0DP,  ISMSTR,  JSPH,    JTUR,
     F   ITY,     JTHE,    JSMS,    NPG ,
     G   glob_therm)
C-----------------------------------------------
C   M o d u l e s
C-----------------------------------------------
      use glob_therm_mod
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      "com08_c.inc"
#include      "param_c.inc"
#include      "scr05_c.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER, INTENT(IN) :: JSMS
      INTEGER, INTENT(IN) :: JTUR
      INTEGER, INTENT(IN) :: ITY
      INTEGER, INTENT(IN) :: JTHE
      INTEGER, INTENT(IN) :: ISMSTR
      INTEGER, INTENT(IN) :: JSPH,NPG
C
      INTEGER NEL
      my_real
     .   PM(NPROPM,*), SIG(NEL,6), RHO(MVSIZ),
     .   D1(*), D2(*), D3(*), D4(*),D5(*), D6(*),
     .   BPRELD(NEL,*), OFF(*)
      INTEGER NELTST,ITYPTST,PID(*),G_DT
      INTEGER MAT(*),NGL(*),IPM(NPROPMI,*)
      my_real
     .   DT2T
 
      my_real
     .   eint(*), qold(*),
     .   vol0(*),stifn(*), offg(*),geo(npropg,*),mumax(*)
      my_real
     .   vnew(*), vd2(*), deltax(*), ssp(*), aire(*), vis(*), 
     .   psh(*), pnew(*),qnew(*) ,ssp_eq(*), dvol(*),
     .   sold1(*), sold2(*), sold3(*), sold4(*), sold5(*), sold6(*),
     .   mssa(*), dmels(*),conde(*),amu(*),vol_avg(*),dtel(*),
     .   rhoref(*)  ,rhosp(*) 
      DOUBLE PRECISION VOL0DP(*) 
      type (glob_therm_) ,intent(inout) :: glob_therm
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I, MX,IBID !, J
      my_real
     .   reduc ,reduc1 , t1, t2 , ee, nu, ggt, treps, p, tp 
      my_real
     .   c1(mvsiz), lmbd(mvsiz), gg(mvsiz),  
     .   rho0(mvsiz), rho0_1,c1_1 
C-----------------------------------------------
      my_real
     .   df,dpdm,facq0,
     .   e1, e2, e3, e4, e5, e6, einc,
     .   bid1, bid2, bid3, dta,  ym, dpdmp
C-----------------------------------------------
      facq0 = one
      reduc1 = em04 !EM03 !ONE 
 
      DO i=1,nel
       t1 = bpreld(i,1)+zep4*(bpreld(i,2)-bpreld(i,1))
       t2 = bpreld(i,1)+zep7*(bpreld(i,2)-bpreld(i,1))
       reduc = reduc1
       tp = tt-t1
       IF(tp > zero ) THEN  
         reduc = min(reduc1*(one-tp/(t2-t1))+tp/(t2-t1), one)
       ENDIF
C       
       mx = mat(i)
       rho0_1 = pm(1, mx)      
       nu = pm(21, mx)
       ee = pm(20, mx)*reduc
       gg(i) = pm(22, mx)*reduc
       c1_1  = pm(32, mx)   
       lmbd(i) = dt1*ee*nu/((one+nu)*(one-two*nu))
       c1(i) = c1_1*reduc 
       rho0(i) =rho0_1  
C
       !SSP(I)=SQRT((ONEP333*GG(I)+C1(I))/RHO0(I))
       ssp(i)=sqrt((onep333*pm(22, mx)+c1_1)/rho0(i))
      ENDDO
C
      IF (jsph==0)THEN
       CALL mqviscb(
     1   pm,      off,     rho,     bid1,
     2   bid2,    ssp,     bid3,    stifn,
     3   dt2t,    neltst,  ityptst, aire,
     4   offg,    geo,     pid,     vnew,
     5   vd2,     deltax,  vis,     d1,
     6   d2,      d3,      pnew,    psh,
     7   mat,     ngl,     qnew,    ssp_eq,
     8   vol0,    mssa,    dmels,   ibid,
     9   facq0,   conde,   dtel,    g_dt,
     a   ipm,     rhoref,  rhosp,   nel,
     b   ity,     ismstr,  jtur,    jthe,
     c   jsms,    npg   ,  glob_therm)
      ELSE
       CALL mdtsph(
     1   pm,      off,     rho,     bid1,
     2   bid2,    bid3,    stifn,   dt2t,
     3   neltst,  ityptst, offg,    geo,
     4   pid,     mumax,   ssp,     vnew,
     5   vd2,     deltax,  vis,     d1,
     6   d2,      d3,      pnew,    psh,
     7   mat,     ngl,     qnew,    ssp_eq,
     8   g_dt,    dtel,    nel,     ity,
     9   jtur,    jthe)
      ENDIF
 
      dta =half*dt1
C       
      DO i=1,nel
        pnew(i)  = c1(i)*amu(i)                         !
        treps    = d1(i)+d2(i)+d3(i)
        ggt      = dt1*gg(i)
        sig(i,1) = sig(i,1)+two*ggt*d1(i)+lmbd(i)*treps
        sig(i,2) = sig(i,2)+two*ggt*d2(i)+lmbd(i)*treps
        sig(i,3) = sig(i,3)+two*ggt*d3(i)+lmbd(i)*treps
        sig(i,4) = sig(i,4)+ggt*d4(i)
        sig(i,5) = sig(i,5)+ggt*d5(i)
        sig(i,6) = sig(i,6)+ggt*d6(i)
        qold(i)  = qnew(i)
      ENDDO
      !write(*,'("***** sboltlaw ")')
      !DO I=1,NEL
      !  write(*,'(i5,6e10.3)')i,sig(i,1:6)
      !enddo
C
C - To avoid discontinuous element pressure when switching material law 
C   cf Law1, Law2, etc <=> Pressure is computed using total formulation P = C1*(rho/rho0-1)
      IF(ismstr == 1)THEN
        IF(iresp==0)THEN
          DO i=1,nel
           IF(tt >= (bpreld(i,2)-two*dt1)) THEN
             p = -third*(sig(i,1)+sig(i,2)+sig(i,3))
             rho(i)=rho0_1*(one+p/c1_1) ! rho will be used in next cycles to compute 
                                       ! rho_new = rho_old - rho0 dV
                                       ! AMU = rho/rho0 - 1 in double precision
           ENDIF
          ENDDO
        ELSE ! Single precision
          DO i=1,nel
           IF(tt >= (bpreld(i,2)-two*dt1)) THEN
             p = -third*(sig(i,1)+sig(i,2)+sig(i,3))
             rho(i)=rho0_1*(one+p/c1_1)
             vol0dp(i) = vol0(i)*(one+p/c1_1) ! VOL0DP will be used in next cycle to compute
                                             ! AMU = VOL0DP/VOLDP - 1 in single precision, VOLDP=V0 in small strain
           ENDIF
          ENDDO
        END IF
      ELSEIF(ismstr==2)THEN
        IF(iresp==0)THEN
          DO i=1,nel
            IF(offg(i) > one)THEN
              IF(tt >= (bpreld(i,2)-two*dt1)) THEN
                p = -third*(sig(i,1)+sig(i,2)+sig(i,3))
                rho(i)=rho0_1*(one+p/c1_1) ! rho will be used in next cycles to compute 
                                          ! rho_new = rho_old - rho0 dV
                                          ! AMU = rho/rho0 - 1 in double precision
              ENDIF
            ELSE
              p = -third*(sig(i,1)+sig(i,2)+sig(i,3))
              vol0(i)   = vnew(i)*(one+p/c1_1)
              rho(i)    = rho0_1*vol0(i)/vnew(i) ! cf rho = rho0 * V0 / V at each time step
                                                 ! AMU = rho/rho0 - 1 in double precision
 
            END IF
          ENDDO
        ELSE ! Single precision
          DO i=1,nel
            IF(offg(i) > one)THEN
              IF(tt >= (bpreld(i,2)-two*dt1)) THEN
                p = -third*(sig(i,1)+sig(i,2)+sig(i,3))
                rho(i)    = rho0_1*(one+p/c1_1)
                vol0dp(i) = vol0(i)*(one+p/c1_1) ! VOL0DP will be used in next cycle to compute
                                                ! AMU = VOL0DP/VOLDP - 1 in single precision, VOLDP=V0 in small strain
               ENDIF
            ELSE
              p = -third*(sig(i,1)+sig(i,2)+sig(i,3))
              vol0(i)   = vnew(i)*(one+p/c1_1)
              vol0dp(i) = vol0(i)  ! cf rho = rho0 * VOL0DP / VOLDP
                                   ! AMU = VOL0DP/VOLDP - 1 in single precision, VOLDP = V in large strain
              rho(i)    = rho0_1*vol0(i)/vnew(i)
            END IF
          ENDDO
        END IF
      ELSEIF(ismstr<=4)THEN
        IF(iresp==0)THEN
          DO i=1,nel
            p = -third*(sig(i,1)+sig(i,2)+sig(i,3))
            vol0(i)=vnew(i)*(one+p/c1_1)
            rho(i)= rho0_1*vol0(i)/vnew(i) ! cf rho = rho0 * V0 / V at each time step
                                           ! AMU = rho/rho0 - 1 in double precision
          ENDDO
        ELSE ! Single precision
          DO i=1,nel
            p = -third*(sig(i,1)+sig(i,2)+sig(i,3))
            vol0(i)   = vnew(i)*(one+p/c1_1)
            rho(i)    = rho0_1*vol0(i)/vnew(i)
            vol0dp(i) = vol0(i)  ! cf rho = rho0 * VOL0DP / VOLDP
                                 ! AMU = VOL0DP/VOLDP - 1 in single precision, VOLDP = V in large strain
          ENDDO
        END IF
      END IF
C
      DO i=1,nel
        e1=d1(i)*(sold1(i)+sig(i,1))
        e2=d2(i)*(sold2(i)+sig(i,2))
        e3=d3(i)*(sold3(i)+sig(i,3))
        e4=d4(i)*(sold4(i)+sig(i,4))
        e5=d5(i)*(sold5(i)+sig(i,5))
        e6=d6(i)*(sold6(i)+sig(i,6))
        einc= zero !VOL_AVG(I)*(E1+E2+E3+E4+E5+E6)*DTA - HALF*DVOL(I)*(QOLD(I)+QNEW(I))
        eint(i)=(eint(i)+einc*off(i)) / max(em15,vol0(i))
      ENDDO
C
      RETURN
      END