balph2.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/.
!||====================================================================
!||    balph2        ../engine/source/ale/bimat/balph2.F
!||--- called by ------------------------------------------------------
!||    bforc2        ../engine/source/ale/bimat/bforc2.F
!||--- calls      -----------------------------------------------------
!||    idp_free      ../engine/source/system/machine.F
!||    idp_lock      ../engine/source/system/machine.F
!||--- uses       -----------------------------------------------------
!||    message_mod   ../engine/share/message_module/message_mod.F
!||====================================================================
      SUBROUTINE balph2(PM,ALPH1,ALPH2,VOLT,FILL,
     .                  SIG1,EINT1,VOLO1,RHON1,FLUX1,FLU11,OFF1,
     .                  SIG2,EINT2,VOLO2,RHON2,FLUX2,FLU12,OFF2,
     .                  SIGT,EINTT,RHOT ,TEMPT,L_TEMP,
     .                  BFRACT,L_BFRAC,PLAST,L_PLAS,VOLN, BULKT, L_BULK,NEL,
     .                  AIRE, AIRES, 
     .                  D1, D2, D3, D4, D5, D6,
     .                  D1S, D2S, D3S, D4S, D5S, D6S,
     .                  MAT, NC1, NC2, NC3, NC4,
     .                  DALPH1, DALPH2)
      USE message_mod
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
#include      "comlock.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      "com04_c.inc"
#include      "com08_c.inc"
#include      "param_c.inc"
#include      "vect01_c.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER L_TEMP, L_PLAS, L_BFRAC,L_BULK,NEL
      INTEGER MAT(*), NC1(*), NC2(*), NC3(*), NC4(*)
      my_real
     .   PM(NPROPM,*), ALPH1(*), ALPH2(*), VOLT(*), FILL(NUMNOD,*),
     .   SIG1(NEL,6), EINT1(*), VOLO1(*), RHON1(*), FLUX1(4,*), FLU11(*),
     .   OFF1(*), SIG2(NEL,6), EINT2(*), VOLO2(*), RHON2(*), FLUX2(4,*),
     .   FLU12(*), OFF2(*), SIGT(NEL,6), EINTT(*), RHOT(*), TEMPT(*),
     .   PLAST(*), BFRACT(*),VOLN(MVSIZ),BULKT(*),
     .   AIRE(*), AIRES(*), 
     .   d1(*), d2(*), d3(*), d4(*), d5(*), d6(*),
     .   d1s(*), d2s(*), d3s(*), d4s(*), d5s(*), d6s(*),
     .   dalph1(*), dalph2(*)
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I, MX1, MX2
      my_real
     .   ALPH1N(MVSIZ), ALPH2N(MVSIZ), ALPN1,
     .   ALPN2, ALPN3, ALPN4, ALPD1, ALPD2, ALPD3, ALPD4, ALPD, ALPN,
     .   EPS, SOFF1, SOFF2, ALPHN, ALPHNN, C11, C12, ALPHT, DA,
     .   IALPH1, IALPH2
C-----------------------------------------------
      DO I=1,nel
        sigt(i,1)=zero
        sigt(i,2)=zero
        sigt(i,3)=zero
        sigt(i,4)=zero
        sigt(i,5)=zero
        sigt(i,6)=zero
        eintt(i) =zero
        rhot(i)  =zero  
        volt(i)  =voln(i)
        aires(i) =aire(i)
        d1s(i)   =d1(i)
        d2s(i)   =d2(i)
        d3s(i)   =d3(i)
        d4s(i)   =d4(i)
        d5s(i)   =d5(i)
        d6s(i)   =d6(i)
      ENDDO
 
      IF(l_temp>0)THEN
        DO i=1,nel
          tempt(i)=zero     
        ENDDO
      ENDIF
 
      IF(l_plas>0)THEN
        DO i=1,nel
          plast(i)=zero     
        ENDDO
      ENDIF
      
      IF(l_bfrac>0)THEN
        DO i=1,nel
          bfract(i)=zero     
        ENDDO
      ENDIF         
      
      IF(l_bulk>0)THEN
        DO i=1,nel
          bulkt(i)=zero     
        ENDDO
      ENDIF            
C-----------------------------------------------
C     MULTIMATERIAUX TAUX DE REMPLISSAGE
C-----------------------------------------------
       DO i=1,nel
         alph1n(i)=zero      
         alph2n(i)=zero        
         alpn1= max(zero,fill(nc1(i),1))
         alpn2= max(zero,fill(nc2(i),1))
         alpn3= max(zero,fill(nc3(i),1))
         alpn4= max(zero,fill(nc4(i),1))
         alpd1=abs(fill(nc1(i),1))
         alpd2=abs(fill(nc2(i),1))
         alpd3=abs(fill(nc3(i),1))
         alpd4=abs(fill(nc4(i),1))
         alpd=alpd1+alpd2+alpd3+alpd4
         alpn=alpn1+alpn2+alpn3+alpn4
         IF (alpd>em20) alph1n(i)=alpn/alpd
       ENDDO
 
       IF(jmult>1)THEN
         DO i=1,nel
           alpn1= max(zero,fill(nc1(i),2))
           alpn2= max(zero,fill(nc2(i),2))
           alpn3= max(zero,fill(nc3(i),2))
           alpn4= max(zero,fill(nc4(i),2))
           alpd1=abs(fill(nc1(i),2))
           alpd2=abs(fill(nc2(i),2))
           alpd3=abs(fill(nc3(i),2))
           alpd4=abs(fill(nc4(i),2))
           alpd=alpd1+alpd2+alpd3+alpd4
           alpn=alpn1+alpn2+alpn3+alpn4
           IF(alpd>em20)alph2n(i)=alpn/alpd  
         ENDDO
       ENDIF
 
       eps=em15
       soff1=zero
       soff2=zero
 
      IF(jmult==1)THEN
       ialph1 = 0
       DO i=1,nel
         dalph1(i)=zero
         IF(alph1n(i)/=zero.AND.alph1n(i)/=one)THEN
          alphn=(volo1(i)-dt1*half*(flu11(i)
     .    +flux1(1,i)+flux1(2,i)+flux1(3,i)+flux1(4,i)))/voln(i)
          alphnn=alphn*(one+(d1(i)+d2(i)+d3(i))*dt1)
          IF((sig1(i,1)+sig1(i,2)+sig1(i,3))>zero)THEN
           alphn= min(alphn,alphnn)
          ELSE
           alphn= max(alphn,alphnn)
          ENDIF
          IF(alphn<=zep99.AND.alphn>zero.AND.rhon1(i)<=zero)THEN
           IF(rhon1(i)/=zero)THEN
#include "lockon.inc"
            WRITE(6,*)' ***NEGATIVE RHO****ALPH1,RHON1'
            WRITE(6,*)i+nft,alph1n(i),rhon1(i)
#include "lockoff.inc"
            rhon1(i)=zero
           ENDIF
           alphn=zero
          ENDIF
          dalph1(i)=alphn-alph1n(i)
          alph1(i)= min(one,alphn)
          alph1(i)= max(zero,alph1(i))
         ELSEIF(alph1n(i)==one.AND.alph1(i)<zep999)THEN
          alphn=(volo1(i)-dt1*half*(flu11(i)
     .    +flux1(1,i)+flux1(2,i)+flux1(3,i)+flux1(4,i)))/voln(i)
          alphnn=alphn*(one+(d1(i)+d2(i)+d3(i))*dt1)
          alphn= max(alphn,alphnn,alph1(i))
          IF(rhon1(i)<=zero)THEN
           rhon1(i)=zero
           alphn=zero
          ENDIF
          alph1(i)= min(one,alphn)
         ELSEIF(alph1n(i)==zero.AND.alph1(i)>em3)THEN
          alph1(i)=alph1n(i)
         ELSE
          alph1(i)=alph1n(i)
         ENDIF
       ENDDO !next I
 
       DO i=1,nel
         IF(alph1(i)<eps)THEN
          alph1(i)=zero
          off1(i) =zero
          volo1(i)=zero
          eint1(i)=zero
         ELSE
          off1(i)=one
         ENDIF
         soff1=soff1+off1(i)
         IF(alph1(i)==zero)THEN
         dalph1(i)=zero
          ialph1=1
         ELSEIF(alph1(i)==one)THEN
          dalph1(i)=zero
          ialph1=1
         ENDIF
       ENDDO
 
       IF(ialph1==1) THEN
         CALL idp_lock(2)
         DO i=1,nel
           IF(alph1(i)==zero)THEN
             fill(nc1(i),1)= min(zero,fill(nc1(i),1))
             fill(nc2(i),1)= min(zero,fill(nc2(i),1))
             fill(nc3(i),1)= min(zero,fill(nc3(i),1))
             fill(nc4(i),1)= min(zero,fill(nc4(i),1))
           ELSEIF(alph1(i)==one)THEN
             fill(nc1(i),1)= max(zero,fill(nc1(i),1))
             fill(nc2(i),1)= max(zero,fill(nc2(i),1))
             fill(nc3(i),1)= max(zero,fill(nc3(i),1))
             fill(nc4(i),1)= max(zero,fill(nc4(i),1))
           ENDIF
         ENDDO
         CALL idp_free(2)
       ENDIF
 
      ELSEIF(jmult==2)THEN
       ialph1 = 0
       ialph2 = 0
       DO i=1,nel
         dalph1(i)=zero
         dalph2(i)=zero
         mx1=nint(pm(21,mat(i)))
         c11 =pm(32,mx1)
         mx2=nint(pm(22,mat(i)))
         c12 =pm(32,mx2)
         IF(alph1n(i)/=zero.AND.alph1n(i)/=one)THEN
          alphn=(volo1(i)-dt1*half*(flu11(i)
     .    +flux1(1,i)+flux1(2,i)+flux1(3,i)+flux1(4,i)))/voln(i)
          alphnn=alphn*(one+(d1(i)+d2(i)+d3(i))*dt1)
          IF((sig1(i,1)+sig1(i,2)+sig1(i,3))>zero)THEN
           alphn= min(alphn,alphnn)
          ELSE
           alphn= max(alphn,alphnn)
          ENDIF
          IF(alphn<=zep99.AND.alphn>zero.AND.rhon1(i)<=zero)THEN
           IF(rhon1(i)/=zero)THEN
#include "lockon.inc"
            WRITE(6,*)' ***NEGATIVE RHO****ALPH1,RHON1'
            WRITE(6,*)i+nft,alph1n(i),rhon1(i)
#include "lockoff.inc"
            rhon1(i)=zero
           ENDIF
           alphn=zero
          ENDIF
          dalph1(i)=alphn-alph1n(i)
          alph1(i)= min(one,alphn)
          alph1(i)= max(zero,alph1(i))
         ELSEIF(alph1n(i)==one.AND.alph1(i)<zep999)THEN
          alphn=(volo1(i)-dt1*half*(flu11(i)
     .    +flux1(1,i)+flux1(2,i)+flux1(3,i)+flux1(4,i)))/voln(i)
          alphnn=alphn*(one + (d1(i)+d2(i)+d3(i))*dt1)
          alphn= max(alphn,alphnn,alph1(i))
          IF(rhon1(i)<=zero)THEN
           rhon1(i)=zero
           alphn=zero
          ENDIF
          alph1(i)= min(one,alphn)
         ELSEIF(alph1n(i)==zero.AND.alph1(i)>em3)THEN
          alph1(i)=alph1n(i)
         ELSE
          alph1(i)=alph1n(i)
         ENDIF
         IF(alph2n(i)/=zero.AND.alph2n(i)/=one)THEN
          alphn=(volo2(i)-dt1*half*(flu12(i)
     .    +flux2(1,i)+flux2(2,i)+flux2(3,i)+flux2(4,i)))/voln(i)
          alphnn=alphn*(one+(d1(i)+d2(i)+d3(i))*dt1)
          IF((sig2(i,1)+sig2(i,2)+sig2(i,3))>one)THEN
           alphn= min(alphn,alphnn)
          ELSE
           alphn= max(alphn,alphnn)
          ENDIF
          IF(alphn<=zep99.AND.alphn>zero.AND.rhon2(i)<=zero)THEN
           IF(rhon2(i)/=zero)THEN
#include "lockon.inc"
            WRITE(6,*)' ***NEGATIVE RHO****ALPH2,RHON2'
            WRITE(6,*)i+nft,alph2n(i),rhon2(i)
#include "lockoff.inc"
            rhon2(i)=zero
           ENDIF
           alphn=zero
          ENDIF
          dalph2(i)=alphn-alph2n(i)
          alph2(i)= min(one,alphn)
          alph2(i)= max(zero,alph2(i))
         ELSEIF(alph2n(i)==one.AND.alph2(i)<=zep999)THEN
          alphn=(volo2(i)-dt1*half*(flu12(i)
     .    +flux2(1,i)+flux2(2,i)+flux2(3,i)+flux2(4,i)))/voln(i)
          alphnn=alphn*(one + (d1(i)+d2(i)+d3(i))*dt1)
          alphn= max(alphn,alphnn,alph2(i))
          IF(rhon2(i)<=zero)THEN
           rhon2(i)=zero
           alphn=zero
          ENDIF
          alph2(i)= min(one,alphn)
         ELSEIF(alph2n(i)==zero.AND.alph2(i)>em3)THEN
          alph2(i)=alph2n(i)
         ELSE
          alph2(i)=alph2n(i)
         ENDIF
         alpht=alph1(i)+alph2(i)
         IF(alpht>one)THEN
           da=(alpht-one)/(c11*alph2(i)+c12*alph1(i))
           alph1(i)=alph1(i)*(one -c12*da)
           alph2(i)=alph2(i)*(one -c11*da)
           dalph1(i)=dalph1(i)-c12*da
           dalph2(i)=dalph2(i)-c11*da
         ENDIF
       ENDDO
 
       DO i=1,nel
         IF(alph1(i)<eps)THEN
          alph1(i)=zero
          off1(i) =zero
          volo1(i)=zero
          eint1(i)=zero
         ELSE
          off1(i) =one
         ENDIF
         soff1=soff1+off1(i)
         IF(alph2(i)<eps)THEN
          alph2(i)=zero
          off2(i) =zero
          volo2(i)=zero
          eint2(i)=zero
         ELSE
          off2(i)=one
         ENDIF
         soff2=soff2+off2(i)
         IF(alph1(i)==zero)THEN
          dalph1(i)=zero
          ialph1 = 1
         ELSEIF(alph1(i)==one)THEN
          dalph1(i)=zero
          ialph1 = 1
         ENDIF
         IF(alph2(i)==zero)THEN
          dalph2(i)=zero
          ialph2 = 1
         ELSEIF(alph2(i)==one)THEN
          dalph2(i)=zero
          ialph2 = 1
         ENDIF
       ENDDO
 
       IF(ialph1==1) THEN
         CALL idp_lock(2)
         DO i=1,nel
           IF(alph1(i)==zero)THEN
             fill(nc1(i),1)= min(zero,fill(nc1(i),1))
             fill(nc2(i),1)= min(zero,fill(nc2(i),1))
             fill(nc3(i),1)= min(zero,fill(nc3(i),1))
             fill(nc4(i),1)= min(zero,fill(nc4(i),1))
           ELSEIF(alph1(i)==one)THEN
             fill(nc1(i),1)= max(zero,fill(nc1(i),1))
             fill(nc2(i),1)= max(zero,fill(nc2(i),1))
             fill(nc3(i),1)= max(zero,fill(nc3(i),1))
             fill(nc4(i),1)= max(zero,fill(nc4(i),1))
           ENDIF
         ENDDO
         CALL idp_free(2)
       ENDIF
 
       IF(ialph2==1) THEN
         CALL idp_lock(2)
         DO i=1,nel
           IF(alph2(i)==zero)THEN
             fill(nc1(i),2)= min(zero,fill(nc1(i),2))
             fill(nc2(i),2)= min(zero,fill(nc2(i),2))
             fill(nc3(i),2)= min(zero,fill(nc3(i),2))
             fill(nc4(i),2)= min(zero,fill(nc4(i),2))
           ELSEIF(alph2(i)==one)THEN
             fill(nc1(i),2)= max(zero,fill(nc1(i),2))
             fill(nc2(i),2)= max(zero,fill(nc2(i),2))
             fill(nc3(i),2)= max(zero,fill(nc3(i),2))
             fill(nc4(i),2)= max(zero,fill(nc4(i),2))
           ENDIF
         ENDDO
         CALL idp_free(2)
       ENDIF
      ENDIF
 
      RETURN
      END