lecm51__check_initial_state.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/.
!||====================================================================
!||    lecm51__check_initial_state   ../starter/source/materials/mat/mat051/lecm51__check_initial_state.F
!||--- called by ------------------------------------------------------
!||    hm_read_mat51                 ../starter/source/materials/mat/mat051/hm_read_mat51.F
!||--- calls      -----------------------------------------------------
!||    ancmsg                        ../starter/source/output/message/message.F
!||--- uses       -----------------------------------------------------
!||    message_mod                   ../starter/share/message_module/message_mod.F
!||====================================================================
      SUBROUTINE lecm51__check_initial_state(   
     .   AV,  R0,     C0,   C1,   C2, C3, C4,   C5,
     .   E0,  PM,   RHO0, RHOR, IEXP, PEXT,IFLG,
     .   A1,     A2,  PLA,
     .   ID,  TITR,
     .   SSP1,SSP2, SSP3, SSP4,
     .   LC1 ,LC2 ,  LC3,  LC4,
     .   P0a          )
      USE message_mod
      USE names_and_titles_mod , ONLY : nchartitle
C-----------------------------------------------
C   D e s c r i p t i o n
C-----------------------------------------------
C   Check initial state of the different material
C
C   1 INIVOL FRAC in [0 1]
C   2 AV1+AV2+AV3+AV4 = 1  (tol 0.1%)
C   3 Negative or null density
C   5 Unbalanced initial Pressure
C   4 Global Density not consistent with material 
C                                       densities
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      my_real   
     . av(4),r0(4), c0(4), c1(4),c2(4), c3(4),
     . a1(4),a2(4),
     . c4(4),c5(4), e0(4), pm(4), pext,  rho0, rhor,
     . ssp(4),lc(4),
     . ssp1,ssp2,ssp3,ssp4,
     . lc1,lc2,lc3,lc4,
     . p0a(4)
      INTEGER IEXP, IFLG,PLA(4)
      INTEGER ID
      CHARACTER(LEN=NCHARTITLE)::TITR
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER i,j, TEST1,TEST2,TEST3,TEST4,TEST5,
     .        test6, imax,ipla(4)
     
      my_real sum, p0b, err
      
      character*8 chain
      character*47 chain3
      character*64 chain1
C===============================================|
      IPLA(1)=pla(1)
      ipla(2)=pla(2)
      ipla(3)=pla(3)
      ipla(4)=pla(4) 
 
      imax=iexp+3
      !IMAX number of phases to be checked 3 or 4
      
      ssp(1:4) = (/ssp1,ssp2,ssp3,ssp4/)
      lc(1:4) = (/ lc1, lc2, lc3, lc4/)      
      
      !0: warning or error not detected
      !1: test failed, message in listing file.
      test1=0    !AV in [0,1]
      test2=0    !Sum(AV(i))=1
      test3=0    !Rho not defined
      test4=0    !Unbalanced Pressure  t=0              
      test5=0    !Global density = mean(rho_i)
      test6=0    !C1,C4,C5 >=0,  C1>0 phase 4
 
      IF(iflg==6)THEN
        sum=0
        DO i=1,imax
          sum=sum+av(i)*r0(i)
        END DO
        IF(sum/=zero)THEN
          rho0=sum
          rhor=rho0
        ELSE
          rho0=one !will be computed in engine, not need for sound speed. This formulation has SSP=EM20
          rhor=one
        ENDIF        
        RETURN
      ENDIF
      !TEST1 (quel que soit iflg)
      ! Initial Volume fraction checked : must be between 0 and 1
      ! Otherwise : error
      DO i=1,imax
        IF ((av(i)<zero).OR.(av(i)>one)) THEN
         chain='SUBMAT-0'
         write(chain(8:8),'(i1)')i
         chain1='INITIAL VOLUMETRIC FRACTION MUST BE BETWEEN 0 AND 1 FOR '//chain(1:8)
         CALL ancmsg(msgid=99,msgtype=msgerror,anmode=aninfo,i1=id,c1=titr,c2=chain1)
         test1=1
         !exit
        END IF
      END DO
 
 
      !TEST2 (quel que soit iflg)
      ! Initial Volume fraction : Sum must be 1
      ! Otherwise : warning
      ! if TEST1/=0 it is useless to do this test)
      sum=zero
      IF (test1==0) THEN
        DO i=1,imax
          sum=sum+av(i)  ! AV(I) in [0,1], SUM=0 <=> AV(I)=0 , i=1..4
        END DO
        IF(sum==zero)THEN
         av(1)=1 
         CALL ancmsg(msgid=98,msgtype=msgwarning,anmode=aninfo_blind_1,i1=id,c1=titr,
     .   c2='INITIAL VOLUMETRIC FRACTIONS ARE NULL , SUBMAT-1 SET TO 100%')
        ELSEIF (abs(1-sum)>em03) THEN
         CALL ancmsg(msgid=99,msgtype=msgerror,anmode=aninfo_blind_1,i1=id,c1=titr,
     .   c2='SUM OF INITIAL VOLUMETRIC FRACTION IS NOT EQUAL TO 1')
         test2=1
        END IF
      END IF
 
 
      !TEST3(whatever is iflg)
      !Negative Density or Missing density material
      ! No negative density
      ! If AV>0 THEN R0 must be >0
      DO i=1,imax
        IF (   ((r0(i)<zero).OR.((av(i)>zero).AND.(r0(i)<=zero)))  ) THEN
         chain='SUBMAT-0'
         write(chain(8:8),'(i1.1)')i
         !print *,"chain(1:8) =", chain(1:8)
         chain3='NULL OR NEGATIVE INITIAL DENSITIES FOR '//chain
         CALL ancmsg(msgid=99,msgtype=msgerror,anmode=aninfo,i1=id,c1=titr,c2=chain3)
         test3=1
         !exit
        END IF
      END DO
   
      
      !TEST4 (whatever is iflg)
      !Unbalanced Pressure, material are checked two by two (Cnp couples)
      !IF (TEST1+TEST2+TEST3==0) THEN !If no error
        sum=zero
        IF(imax==4)p0a(4) = c0(4)
        DO i=1,imax-1
        if (r0(i)==zero)cycle !if empty submaterial
        p0a(i)=c0(i)+c4(i)*e0(i)
          DO j=i,imax
          IF (r0(j)==zero)cycle !if empty submaterial
          !comparing Pi(t=0) with Pj(t=0)
            p0b=c0(j)+c4(j)*e0(j)
            err=abs(p0a(i)-p0b)
          IF (err<=em20) err=zero  !avoid epsilon issue if P0a ou P0b is zero
          IF (err/=zero) THEN  ! JWL under pressure may detonate.
            IF (err/max(abs(p0a(i)),abs(p0b))>em06) THEN
                IF (j==4) THEN
                  test4=1
                  CALL ancmsg(msgid=98,msgtype=msgwarning,anmode=aninfo_blind_1,
     .                        i1=id,c1=titr,c2='PRESSURES ARE UNBALANCED WITH JWL MATERIAL')
                  exit
              ELSE
                  test4=1
                  CALL ancmsg(msgid=98,msgtype=msgwarning,anmode=aninfo_blind_1,
     .                        i1=id,c1=titr,c2='INITIAL PRESSURES ARE UNBALANCED')
                EXIT
            END IF
            END IF
            END IF
          END DO !next J
          IF(test4==1)EXIT
        END DO !next I
      !END IF
      
      
      !TEST5 (whatever is iflg)
      !Check the consistency between global density and materials densities.
      !IF(TEST1+TEST2+TEST3==0)THEN !If no error
        sum=0
        DO i=1,imax
          sum=sum+av(i)*r0(i)
        END DO
        test5=1
        rho0=sum
        rhor=rho0
      !END IF
 
 
      !TEST6 (whatever is iflg)
      ! Check if C1>=0, C4>=0, C5>=0
      ! Otherwise : error
      DO i=1,3
        IF ((c1(i)<zero)) THEN
         chain='SUBMAT-0'
         write(chain(8:8),'(i1.1)')i
         chain1='POLYNOMIAL COEFFICIENT C1 MUST BE POSITIVE OR NULL FOR '//chain
         CALL ancmsg(msgid=99,msgtype=msgerror,anmode=aninfo,i1=id,c1=titr,c2=chain1)
         test6=1
         !exit
        END IF
 
        IF ((c2(i)<zero)) THEN
         chain='SUBMAT-0'
         write(chain(8:8),'(i1.1)')i
         chain1='POLYNOMIAL COEFFICIENT C2 IS NEGATIVE FOR '//chain
         CALL ancmsg(msgid=98,msgtype=msgwarning,anmode=aninfo_blind_1,i1=id,c1=titr,c2=chain1)
         test6=1
         !exit
        END IF
        
        IF ((c3(i)<zero)) THEN
         chain='SUBMAT-0'
         write(chain(8:8),'(i1.1)')i
         chain1='POLYNOMIAL COEFFICIENT C3 IS NEGATIVE FOR '//chain
         CALL ancmsg(msgid=98,msgtype=msgwarning,anmode=aninfo_blind_1,i1=id,c1=titr,c2=chain1)
         test6=1
         !exit
        END IF        
 
        IF ((c4(i)<zero)) THEN
         chain='SUBMAT-0'
         write(chain(8:8),'(i1.1)')i
         chain1='POLYNOMIAL COEFFICIENT C4 MUST BE POSITIVE OR NULL FOR '//chain
         CALL ancmsg(msgid=99,msgtype=msgerror,anmode=aninfo_blind_1,i1=id,c1=titr,c2=chain1)
         test6=1
         !exit
        END IF        
 
        IF ((c5(i)<zero)) THEN
         chain='SUBMAT-0'
         write(chain(8:8),'(i1.1)')i
         chain1='POLYNOMIAL COEFFICIENT C5 MUST BE POSITIVE OR NULL FOR '//chain
         CALL ancmsg(msgid=99,msgtype=msgerror,anmode=aninfo_blind_1,i1=id,c1=titr,c2=chain1)
         test6=1
         !exit
        END IF                        
      END DO !next I
      
          
      !check pressure cut-off pressure for perfect gas
      IF(iflg<=1)THEN
        DO i=1,3
        IF(ipla(i)>0)cycle
        p0a(i) = pext+c0(i)+c4(i)*e0(i)
          IF(((pext+c0(i))==c1(i)).AND.(c2(i)==zero).AND.
     .    (c3(i)==zero).AND.(c4(i)==c5(i)).AND.
     .    (c4(i)>zero).AND.(p0a(i)>zero))THEN
           !we have a perfect gas
           IF(pm(i)/=zero .AND. pext/=zero)THEN
             IF(abs(pm(i)+pext)>em06)THEN
               !PM(I)=-PEXT 
               chain='SUBMAT-0'
               write(chain(8:8),'(i1.1)')i           
               chain1=
     .         chain//' IS A PERFECT GAS:MINIMUM PRESSURE SHOULD BE -PEXT       '
               CALL ancmsg(msgid=98,msgtype=msgwarning,anmode=aninfo_blind_1,i1=id,c1=titr,c2=chain1)
             END IF
           ENDIF
         END IF
        END DO !next I
      END IF !(IFLG<=1)
      
      !check Drucker Prager Yield Surface
      IF(iflg<=1)THEN
        DO i=1,3
          chain='SUBMAT-0'
          write(chain(8:8),'(i1.1)')i           
          !
          IF(a1(i) < zero .AND. a2(i) == zero)THEN                               
            chain1 = chain//': INVERTED YIELD SURFACE. CHECK A1 SIGN.                '
            CALL ancmsg(msgid=829,msgtype=msgwarning,anmode=aninfo,i1=51,i2=id,c1='WARNING',c2=titr,c3=chain)                                              
          ENDIF                                                                  
          IF(a2(i) < zero)THEN                                                   
            chain1 = chain//': UNTIPYCAL YIELD SURFACE. CHECK A2 SIGN.               '        
            CALL ancmsg(msgid=829,msgtype=msgwarning,anmode=aninfo,i1=51,i2=id,c1='WARNING',c2=titr,c3=chain)                                              
          ENDIF                                                                  
          !
        END DO !next I
      END IF !(IFLG<=1)      
      
      RETURN
      END
C======================================================================|