chsti3.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/.
!||====================================================================
!||    chsti3       ../engine/source/elements/shell/coque/chsti3.F
!||--- called by ------------------------------------------------------
!||    cforc3       ../engine/source/elements/shell/coque/cforc3.F
!||    cforc3_crk   ../engine/source/elements/xfem/cforc3_crk.F
!||====================================================================
      SUBROUTINE chsti3(JFT ,JLT,THK,HOUR,OFF ,PX1,
     2                  PX2,PY1,PY2 ,SIGY,IXC,DT1C,
     3                  SSP,RHO,STI ,Z2 ,EANI,STIR,
     4                  SHF,THK0,THK02,VISCMX,G  ,A11 ,
     5   H1  ,H2  ,H3 ,YM   ,NU  , ALPE, 
     6   VHX ,VHY ,VX1 ,VX2 ,VX3 ,VX4 ,VY1 ,
     7   VY2 ,VY3 ,VY4,VZ1  ,VZ2  ,VZ3 ,VZ4 ,AREA,
     8   H11 ,H12 ,H13 ,H21 ,H22 ,H23  ,H31  ,H32  ,H33  ,
     9   B11 ,B12 ,B13 ,B14 ,B21 ,B22  ,B23  ,B24  ,
     A   RX1 ,RX2 ,RX3 ,RX4 ,RY1 ,RY2  ,RY3  ,RY4,
     B   IPARTC ,PARTSAV,
     C   IHBE   ,NFT   ,ISMSTR ,SRH3 ,IGTYP ,
     D   IGMAT  ,A11R  ,NEL)
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      "com06_c.inc"
#include      "param_c.inc"
#include      "remesh_c.inc"
#include      "scr02_c.inc"
#include      "scr14_c.inc"
#include      "scr16_c.inc"
#include      "scr17_c.inc"
#include      "sms_c.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER IXC(NIXC,*),IPARTC(*), JFT, JLT,
     .              IHBE   ,NFT     ,ISMSTR,IGTYP,IGMAT,NEL
C     REAL
      my_real
     .   THK(*), HOUR(NEL,5), OFF(*),EANI(*),PARTSAV(NPSAV,*),
     .   PX1(*), PX2(*), PY1(*), PY2(*), SIGY(*),DT1C(*),STIR(*),
     .   SSP(MVSIZ), RHO(MVSIZ),STI(MVSIZ),Z2(MVSIZ), SHF(MVSIZ),
     .   THK0(MVSIZ),VISCMX(MVSIZ), ALPE(MVSIZ), A11R(*)
C     REAL
      my_real
     .   B11(MVSIZ), B12(MVSIZ), B13(MVSIZ), B14(MVSIZ),
     .   B21(MVSIZ), B22(MVSIZ), B23(MVSIZ), B24(MVSIZ),
     .   h11(mvsiz), h12(mvsiz), h13(mvsiz), h14(mvsiz),
     .   h21(mvsiz), h22(mvsiz), h23(mvsiz), h24(mvsiz),
     .   h31(mvsiz), h32(mvsiz), h33(mvsiz), h34(mvsiz),
     .   rx1(mvsiz), rx2(mvsiz), rx3(mvsiz), rx4(mvsiz),
     .   ry1(mvsiz), ry2(mvsiz), ry3(mvsiz), ry4(mvsiz),
     .   vx1(mvsiz), vx2(mvsiz), vx3(mvsiz), vx4(mvsiz),
     .   vy1(mvsiz), vy2(mvsiz), vy3(mvsiz), vy4(mvsiz),
     .   vz1(mvsiz), vz2(mvsiz), vz3(mvsiz), vz4(mvsiz),
     .   area(mvsiz),vhx(mvsiz),
     .   h1(mvsiz),  h2(mvsiz),  h3(mvsiz), vhy(mvsiz),
     .   ym(mvsiz), nu(mvsiz), thk02(mvsiz),srh3(*)
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I, MX, J
C     REAL
      my_real
     .   GAMA1(MVSIZ), GAMA2(MVSIZ), GAMA3(MVSIZ), GAMA4(MVSIZ),
     .   H3L(MVSIZ), H3Q(MVSIZ),
     .   HG1(MVSIZ), HG2(MVSIZ) ,HP1(MVSIZ), HP2(MVSIZ) , 
     .   G(MVSIZ) ,  B1(MVSIZ),  B2(MVSIZ),  
     .   A11(MVSIZ) ,EHOU(MVSIZ),
     .   PX1V, PX2V, PY1V, PY2V, SSPHVL,ZZ2,VV, 
     .   SHFPR3, FAC, BB, HOUR1(MVSIZ), HOUR2(MVSIZ), HOUR3(MVSIZ), 
     .   EHOURT, PX1VP,PX2VP,PY1VP,PY2VP, FMAX(MVSIZ),
     .   HOUR1P,HOUR2P,HOUR3P,FMAX1,FMAXP1,FMAX2,FMAXP2,FMAXP,
     .   scale(mvsiz)
C-----------------------------------------------
C
      ehourt = zero
C
      IF(ismstr/=3.AND.ihbe/=0)THEN
       DO i=jft,jlt
         px1v    = px1(i)*vhx(i)
         px2v    = px2(i)*vhx(i)
         py1v    = py1(i)*vhy(i)
         py2v    = py2(i)*vhy(i)
         gama1(i)= off(i)*( one - px1v-py1v)
         gama3(i)= off(i)*( one + px1v+py1v)
         gama2(i)= off(i)*(-one - px2v-py2v)
         gama4(i)= off(i)*(-one + px2v+py2v)   
       ENDDO
      ELSE
       DO i=jft,jlt
         gama1(i)= off(i)
         gama3(i)= off(i)
         gama2(i)= -off(i)
         gama4(i)= -off(i)
       ENDDO
      ENDIF
C
      DO i=jft,jlt
       hp1(i) =h1(i)
       hp2(i) =h2(i)
       h3(i) =srh3(i)
       shfpr3 = shf(i)/(three*(one+nu(i)))
       fac    = fourth*rho(i)*thk02(i)
       h3l(i) = h3(i)*zep072169*fac*area(i)
       h3q(i) = h3(i)*h3l(i)*hundred
       h3l(i) = h3l(i)*ssp(i)
C
       fac    = ym(i)*thk0(i)*dt1c(i)*one_over_8
       h1(i)  = hp1(i) * fac
       b1(i) = px1(i)*px1(i)+py1(i)*py1(i)
       b2(i) = px2(i)*px2(i)+py2(i)*py2(i)
       h2(i)  = hp2(i) * fac*thk02(i)*shfpr3/(b1(i)+b2(i))
      ENDDO
C-----------------------------------
C     TRIANGLES
C-----------------------------------
      DO i=jft,jlt
       IF(ixc(4,i)==ixc(5,i))THEN
         h3q(i)=zero
         h3l(i)=zero
         h1(i)=zero
         h2(i)=zero   
       ENDIF
      ENDDO
C-----------------------------------
C     STIFFNESS - DT
C-----------------------------------
       DO i=jft,jlt
          scale(i) = zero
      ENDDO
      IF(nodadt/=0.OR.idt1sh==1.OR.idtmins==2)THEN
        DO i=jft,jlt
          scale(i) = +max(gama1(i)*gama1(i),gama2(i)*gama2(i),
     .              gama3(i)*gama3(i),gama4(i)*gama4(i)) *
     .    dt1c(i)*max(h1(i),h2(i),h3l(i)) /
     .              max(dt1c(i)*dt1c(i),em20)
          sti(i)=sti(i) + scale(i) 
        ENDDO
C     
        !!IF(IDTMINS==2)IGMAT = 1
        IF(igtyp == 52 .OR.
     .    ((igtyp == 11 .OR. igtyp == 17 .OR. igtyp == 51) 
     .                                       .AND. igmat > 0)) THEN
          IF(nadmesh==0)THEN
             DO i=jft,jlt
              IF (off(i)==zero) THEN
               sti(i) = zero
               stir(i) = zero          
             ELSE
               vv = viscmx(i) * viscmx(i)  * alpe(i)
               fac = max(b1(i),b2(i))/(area(i) * vv)
               sti(i) =  sti(i) + fac * thk0(i) * a11(i) 
               stir(i) = fac *  a11r(i)*thk0(i)*thk02(i)*one_over_12 +
     .                   fac *  a11(i) *thk0(i)*area(i) *one_over_9 +
     .                   scale(i)*(thk02(i)* one_over_12 +  area(i)*one_over_9)
             ENDIF 
             ENDDO
         ELSE
            DO i=jft,jlt
             IF (off(i)==zero) THEN
              sti(i) = zero
              stir(i) = zero          
            ELSE
              vv = viscmx(i) * viscmx(i)  * alpe(i)
               fac = max(b1(i),b2(i))/(area(i) * vv)
               sti(i) =  sti(i) + fac * thk0(i) * a11(i) 
               stir(i) = fac *  a11r(i)*thk0(i)*thk02(i)* one_over_12 +
     .                           scale(i)*thk02(i)* one_over_12
            ENDIF 
            ENDDO
         END IF
        ELSE ! IGTYP
          IF(nadmesh==0)THEN
             DO i=jft,jlt
              IF (off(i)==zero) THEN
               sti(i) = zero
               stir(i) = zero          
             ELSE
               vv = viscmx(i) * viscmx(i)  * alpe(i)
               sti(i) =  sti(i) + max(b1(i),b2(i)) 
     .                      * thk0(i) * a11(i) / (area(i) * vv)
               stir(i) = sti(i) * (thk02(i)* one_over_12 + area(i) * one_over_9)
c     ..                + 0.5 * SHF(I) * AREA(I) * G(I)/A11(I)) 
             ENDIF 
             ENDDO
         ELSE
            DO i=jft,jlt
             IF (off(i)==zero) THEN
              sti(i) = zero
              stir(i) = zero          
            ELSE
              vv = viscmx(i) * viscmx(i)  * alpe(i)
              sti(i) =  sti(i) + max(b1(i),b2(i)) 
     .                     * thk0(i) * a11(i) / (area(i) * vv)
              stir(i) = sti(i) * thk02(i)* one_over_12
            ENDIF 
            ENDDO
        END IF
       endif! igtyp
      ENDIF ! 
C-----------------------------------
C     ANTI-HOURGLASS MEMBRANE FORCES
C-----------------------------------
      DO i=jft,jlt
        hg1(i)=vx1(i)*gama1(i)+vx2(i)*gama2(i)
     .        +vx3(i)*gama3(i)+vx4(i)*gama4(i)
        hg2(i)=vy1(i)*gama1(i)+vy2(i)*gama2(i)
     .        +vy3(i)*gama3(i)+vy4(i)*gama4(i)
      ENDDO
 
      DO i=jft,jlt
        hour(i,1)=hour(i,1)+hg1(i)*h1(i)
        hour(i,2)=hour(i,2)+hg2(i)*h1(i)
        fmax(i) = half*sigy(i)*thk0(i)*sqrt(area(i))*hp1(i)
        hour(i,1)=sign( min( abs(hour(i,1)),fmax(i) ) , hour(i,1))
        hour(i,2)=sign( min( abs(hour(i,2)),fmax(i) ) , hour(i,2))
        hour1(i)=hour(i,1)
        h11(i)=hour1(i)*gama1(i)
        h12(i)=hour1(i)*gama2(i)
        h13(i)=hour1(i)*gama3(i)
        h14(i)=hour1(i)*gama4(i)
        hour2(i)=hour(i,2)
        h21(i)=hour2(i)*gama1(i)
        h22(i)=hour2(i)*gama2(i)
        h23(i)=hour2(i)*gama3(i)
        h24(i)=hour2(i)*gama4(i)
        ehou(i) = hour1(i)*hg1(i) + hour2(i)*hg2(i)
      ENDDO
C
C----------------------------------
C     ANTI-HOURGLASS BENDING FORCES
C----------------------------------
      DO i=jft,jlt
        hg1(i)=vz1(i)*gama1(i)+vz2(i)*gama2(i)
     .        +vz3(i)*gama3(i)+vz4(i)*gama4(i)
      ENDDO
      DO i=jft,jlt
       hour(i,3)=hour(i,3)+hg1(i)*h2(i)
       fmax(i) = fourth*sigy(i)*thk0(i)*thk0(i)*hp2(i)
       hour(i,3)=sign( min( abs(hour(i,3)),fmax(i) ) , hour(i,3))
       hour3(i) =hour(i,3)
       h31(i)=hour3(i)*gama1(i)
       h32(i)=hour3(i)*gama2(i)
       h33(i)=hour3(i)*gama3(i)
       h34(i)=hour3(i)*gama4(i)
       ehou(i) = ehou(i) + hour3(i)*hg1(i)
      ENDDO
C
C------------------------
C     ANTI-HOURGLASS MOMENTS
C------------------------
      DO i=jft,jlt
       hg1(i)=rx1(i)-rx2(i)+rx3(i)-rx4(i)
       hg2(i)=ry1(i)-ry2(i)+ry3(i)-ry4(i)
      ENDDO
C
      DO i=jft,jlt
       hour(i,4)=hg1(i)*(h3l(i)+h3q(i)*abs(hg1(i)))
       hour(i,5)=hg2(i)*(h3l(i)+h3q(i)*abs(hg2(i)))
       ehou(i) = ehou(i) + hour(i,4)*hg1(i) + hour(i,5)*hg2(i)
       ehou(i) = dt1c(i) * ehou(i)*off(i)
       ehourt = ehourt + ehou(i)
       b11(i)= hour(i,4)*off(i)
       b12(i)=-hour(i,4)*off(i)
       b13(i)= hour(i,4)*off(i)
       b14(i)=-hour(i,4)*off(i)
       b21(i)= hour(i,5)*off(i)
       b22(i)=-hour(i,5)*off(i)
       b23(i)= hour(i,5)*off(i)
       b24(i)=-hour(i,5)*off(i)
      ENDDO
C
      DO i=jft,jlt
         mx = ipartc(i)
         partsav(8,mx)=partsav(8,mx) + ehou(i)
      ENDDO
C
!$OMP ATOMIC
        ehour = ehour + ehourt
C
      DO i=jft,jlt
        eani(nft+numels+i) = eani(nft+numels+i)+ehou(i)
      ENDDO
C
      RETURN
      END