thermexpc.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/.
!||====================================================================
!||    thermexpc      ../engine/source/materials/mat_share/thermexpc.F
!||--- called by ------------------------------------------------------
!||    cmain3         ../engine/source/materials/mat_share/cmain3.F
!||--- calls      -----------------------------------------------------
!||    finter         ../engine/source/tools/curve/finter.f
!||--- uses       -----------------------------------------------------
!||    elbufdef_mod   ../common_source/modules/mat_elem/elbufdef_mod.F90
!||====================================================================
      SUBROUTINE thermexpc(ELBUF_STR,
     1                       JFT  ,JLT  ,FORTH  ,FOR    ,EINT  ,
     2                       OFF  ,ETH  ,THK0   ,EXX    , EYY  ,  
     3                       PM   ,NPT  ,AREA   ,A1     ,A2    , 
     4                       MAT  ,MTN  ,EINTTH ,DIR    ,IR    ,
     5                       IS   ,NLAY ,THK    ,NEL    ,IGTYP ,
     6                       NPF  , TF  ,IPM   , TEMPEL , DTEMP ,
     7                       THKLY ,POSLY ,MOM, MATLY)
C-----------------------------------------------
C   M o d u l e s
C-----------------------------------------------
      USE elbufdef_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      "param_c.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER JFT, JLT,NPT,MTN,NSG,IR,IS,NLAY,
     .     MAT(MVSIZ),NEL,MATLY(*),IPM(NPROPMI,*),IGTYP,
     .     NPF(*)
C     REAL
      my_real
     .   FOR(NEL,5), FORTH(NEL,2), EINT(JLT,2), EINTTH(*),
     .   OFF(*),EINTT(MVSIZ),DIR(NEL,2),THK(*), PM(NPROPM,*),
     .   thkly(*),dtemp(*),tempel(*),tf(*),posly(mvsiz,*),
     .   mom(nel,3)
      my_real
     .   exx(mvsiz), eyy(mvsiz), eth(mvsiz),thk0(mvsiz) ,
     .   area(mvsiz),pla(mvsiz),a1(mvsiz), a2(mvsiz)
      TYPE(elbuf_struct_), TARGET :: ELBUF_STR
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER  I, MX, J,IORTH,IPT,ILAY,II(5),NPTT,IT,
     .         IORTH_LAY,IFUNC_ALPHA,J1,J2,J3,JJ,IPT_ALL
C     REAL
      my_real
     .   dtinv, amu, fact, visc,m12,anu12,anu21,s1,s2,
     .   fscal_alpha,alpha,deintth,df,wmc,aa,thklay,kxx,kyy
      my_real
     .   nu(mvsiz),e(mvsiz),a12(mvsiz),
     .   e11(mvsiz), e22(mvsiz),
     .   b3(mvsiz), degmb(mvsiz), degfx(mvsiz),degmbth(mvsiz),
     .   ezz(mvsiz),einf(mvsiz),p(mvsiz),degfxth(mvsiz),
     .   p1(mvsiz),p2(mvsiz),ethke(mvsiz),zi2(mvsiz),sumalz(mvsiz)
      TYPE(buf_lay_) ,POINTER :: BUFLY
      TYPE(l_bufel_) ,POINTER :: LBUF
C-----------------------------------------------    
      my_real FINTER 
      EXTERNAL FINTER
C-----------------------------------------------  
      DO I=1,5
        ii(i) = nel*(i-1)
      ENDDO
C-------- 1st Step : Elasticity matrix-----------------------
      iorth_lay = 0
      iorth= -1! not activated
      IF(igtyp == 11 .OR. igtyp == 17 .OR. igtyp == 51 .OR. igtyp == 52) THEN
        iorth_lay = 1
      ELSEIF (mtn==19.OR.mtn==15.OR.mtn==25) THEN
        iorth=1
      ELSE
        iorth=0
      ENDIF
C
      IF(iorth == 1) THEN
        mx  =mat(jft)
        DO i=jft,jlt
          e(i) = pm(20,mx)
          nu(i) = pm(21,mx)
          a1(i) = pm(24,mx)
          a2(i) = pm(25,mx)
C
          e11(i)  = pm(33,mx)
          e22(i)  = pm(34,mx)
          anu12   = pm(35,mx)
          anu21   = pm(36,mx)
C
          a12(i) = (1.-anu12*anu21)
          a1(i) = e11(i)/a12(i)
          a2(i) = e22(i)/a12(i)
          a12(i) = anu21*a1(i)
        ENDDO
      ENDIF
C-------- 2nd Step : Thermal stress computation -----------------------
 
        IF(iorth ==0)THEN
          mx  =mat(jft)
          DO i=jft,jlt
            p(i) =(a1(i)+a2(i))*eth(i)
            forth(i,1)=forth(i,1)+ p(i)
            forth(i,2)=forth(i,2)+ p(i)
          END DO
        ELSEIF(iorth == 1) THEN
          DO i=jft,jlt
             p1(i) = a1(i)*eth(i)+a12(i)*eth(i)
             p2(i) = a12(i)*eth(i)+a2(i)*eth(i)
             s1 = dir(i,1)*dir(i,1)*p1(i)
     .          + dir(i,2)*dir(i,2)*p2(i)
             s2 = dir(i,2)*dir(i,2)*p1(i)
     .          + dir(i,1)*dir(i,1)*p2(i)
             forth(i,1)=forth(i,1)+ s1
             forth(i,2)=forth(i,2)+ s2
          END DO   
        ENDIF 
C
C
      IF(iorth_lay > 0 ) THEN
         ethke(jft : jlt) = zero
        IF(mtn == 15 .OR. mtn == 25) THEN
          ipt_all = 0
          DO ilay=1,nlay
             nptt = elbuf_str%BUFLY(ilay)%NPTT
             j1 = 1+(ilay-1)*jlt     ! JMLY
             j3 = 1+(ilay-1)*jlt*2     ! jdir
            DO it=1,nptt
              lbuf => elbuf_str%BUFLY(ilay)%LBUF(ir,is,it) 
              ipt = ipt_all + it        ! count all NPTT through all layers
              j2 = 1+(ipt-1)*jlt        ! JPOS
              DO i=jft,jlt 
                jj = j2 - 1 + i
                mx = matly(j1+i-1)
C
                e11(i)  = pm(33,mx)
                e22(i)  = pm(34,mx)
                anu12   = pm(35,mx)
                anu21   = pm(36,mx)
C
                ifunc_alpha = ipm(219, mx) 
                fscal_alpha = pm(191, mx)
                alpha = finter(ifunc_alpha,tempel(i),npf,tf,df)
                alpha = alpha * fscal_alpha
                eth(i) = alpha*dtemp(i)
C                  
                a12(i) = (one - anu12*anu21)
                a1(i)  = e11(i)/a12(i)
                a2(i)  = e22(i)/a12(i)
                a12(i) =   anu21*a1(i)
C
                p1(i) = a1(i)*eth(i ) + a12(i)*eth(i)
                p2(i) = a12(i)*eth(i) + a2(i)*eth(i) 
C                
                lbuf%SIG(ii(1)+i)=lbuf%SIG(ii(1)+i) - p1(i) 
                lbuf%SIG(ii(2)+i)=lbuf%SIG(ii(2)+i) - p2(i)
C                
                for(i,1)=for(i,1) - thkly(jj)*p1(i) 
                for(i,2)=for(i,2) - thkly(jj)*p2(i)
C                
                wmc = posly(i,ipt)*thkly(jj)
                mom(i,1) = mom(i,1) - wmc*p1(i)
                mom(i,2) = mom(i,2) - wmc*p2(i)  
C       
                thklay = thkly(jj)*thk0(i)
                ethke(i) = ethke(i) + thklay*eth(i)             
              ENDDO
            ENDDO !!  DO IT=1,NPTT
             ipt_all = ipt_all + nptt
          ENDDO  ! DO ILAY=1,NLAY  
         ELSEIF(mtn > 26) THEN
           aa = zero
           ipt_all = 0
           zi2(jft:jlt) = zero
           sumalz(jft:jlt) = zero
           DO ilay=1,nlay
            nptt = elbuf_str%BUFLY(ilay)%NPTT
            j1 = 1+(ilay-1)*jlt     ! JMLY
            j3 = 1+(ilay-1)*jlt*2 
            DO it=1,nptt
              lbuf => elbuf_str%BUFLY(ilay)%LBUF(ir,is,it) 
              ipt = ipt_all + it        ! count all NPTT through all layers
              j2 = 1+(ipt-1)*jlt        ! JPOS
              DO i=jft,jlt
                jj = j2 - 1 + i
                mx = matly(j1+i-1)
                e(i)  = pm(20,mx)
                nu(i) = pm(21,mx)
                a1(i) = pm(24,mx)
                a2(i) = pm(25,mx)
                a1(i) = e(i)/ (one - nu(i)*nu(i))
                a2(i) = nu(i)*a1(i)
C                    
                ifunc_alpha = ipm(219, mx) 
                fscal_alpha = pm(191, mx)
                alpha = finter(ifunc_alpha,tempel(i),npf,tf,df)
                alpha = alpha * fscal_alpha
                eth(i) = alpha*dtemp(i)    
C            
                p(i) = a1(i)*eth(i) + a2(i)*eth(i)              
C                
                lbuf%SIG(ii(1)+i)=lbuf%SIG(ii(1)+i) - p(i) 
                lbuf%SIG(ii(2)+i)=lbuf%SIG(ii(2)+i) - p(i)
C forces  
                for(i,1)=for(i,1) - thkly(jj)*p(i) 
                for(i,2)=for(i,2) - thkly(jj)*p(i)
C                                                        
                wmc = posly(i,ipt)*thkly(jj)
                mom(i,1) = mom(i,1) - wmc*p(i)
                mom(i,2) = mom(i,2) - wmc*p(i)                
!!                ZI2(I) = ZI2(I) +  POSLY(I,IPT)**2
!!                SUMALZ(I) = SUMALZ(I) + POSLY(I,IPT)* ETH(I)
C      
                thklay = thkly(jj)*thk0(i)
                ethke(i) = ethke(i) + thklay*eth(i)
              ENDDO ! I 
            ENDDO !!  DO IT=1,NPTT
             ipt_all = ipt_all + nptt
          ENDDO  ! DO ILAY=1,NLAY
        ENDIF
          
         DO i=jft,jlt
C-------- 3rd Step : Energies computation -----------------------
!!             KXX = SUMALZ(I)/ ZI2(I)
!!             KYY = KXX  
!!    DEGFXTH(I) = -( MOM(I,1)*KXX + MOM(I,2)*KYY)*HALF*THK0(I)*THK0(I)*AREA(I) ! depending to kxx et kyy
             degmbth(i) = -(for(i,1)+for(i,2))*ethke(i)*half*area(i)
             eintth(i) = eintth(i) + degmbth(i)
             eint(i,1) = eint(i,1) + degmbth(i)
!!             EINT(I,2) = EINT(I,2) + DEGFXTH(I)
C------Thickness change due to thermal expansion-------           
             thk(i) = (thk(i)  + ethke(i))*off(i)
          ENDDO        
      ELSE
        IF(npt/=0) THEN
          DO ilay=1,nlay
            nptt = elbuf_str%BUFLY(ilay)%NPTT
            DO it=1,nptt
              lbuf => elbuf_str%BUFLY(ilay)%LBUF(ir,is,it) 
              IF(iorth ==0)THEN
                DO i=jft,jlt
                  lbuf%SIG(ii(1)+i)=lbuf%SIG(ii(1)+i)-p(i) 
                  lbuf%SIG(ii(2)+i)=lbuf%SIG(ii(2)+i)-p(i) 
                ENDDO
              ELSE
                DO i=jft,jlt
                  lbuf%SIG(ii(1)+i)=lbuf%SIG(ii(1)+i)-p1(i) 
                  lbuf%SIG(ii(2)+i)=lbuf%SIG(ii(2)+i)-p2(i) 
                ENDDO
              ENDIF
            ENDDO
          ENDDO     
         ENDIF
         IF(iorth ==0)THEN
           DO i=jft,jlt
              for(i,1)=for(i,1) - p(i) 
              for(i,2)=for(i,2) - p(i)  
           ENDDO
         ELSE
           DO i=jft,jlt
              for(i,1)=for(i,1) - p1(i) 
              for(i,2)=for(i,2) - p2(i) 
            ENDDO
         ENDIF
C
C-------- 3rd Step : Energies computation -----------------------
         DO i=jft,jlt
           degmbth(i) = -(for(i,1)+for(i,2))*eth(i)*half*thk0(i)*area(i)
         END DO
 
         DO i=jft,jlt
          eintth(i) = eintth(i) + degmbth(i)
          eint(i,1) = eint(i,1) + degmbth(i)
         ENDDO
C------Thickness change due to thermal expansion-------
          DO i=jft,jlt
            thk(i) = thk(i) *(1 + eth(i))*off(i)
          ENDDO          
        ENDIF ! IORTH_LAY > 0
C
      RETURN
      END