adiff3.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/.
!||====================================================================
!||    adiff3                 ../engine/source/ale/ale3d/adiff3.F
!||--- called by ------------------------------------------------------
!||    aeturb                 ../engine/source/ale/turbulence/aeturb.F
!||    akturb                 ../engine/source/ale/turbulence/akturb.F
!||    atherm                 ../engine/source/ale/atherm.f
!||--- uses       -----------------------------------------------------
!||    ale_connectivity_mod   ../common_source/modules/ale/ale_connectivity_mod.F
!||====================================================================
      SUBROUTINE adiff3(PHIN,PHI,GRAD,ALPHA,ALE_CONNECT,VOL,TEMP,RHOCP,NEL)
C-----------------------------------------------
C   D e s c r i p t i o n
C-----------------------------------------------
C Solving Heat equation
C    alpha = k/rhocp is thermal diffusivity
c    coeff is k (factor simplification since originally only EINT/V was updated)
c       PHIN is EINT/V : updated at the end of the subroutine
C       TEMP : updated at the end of the subroutine.
C              Since temperature is calculated incrementally, both energy and temperature must be consistently updated.
C-----------------------------------------------
C   M o d u l e s
C-----------------------------------------------
      USE ale_connectivity_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   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER,INTENT(IN) :: NEL
      my_real,INTENT(INOUT) :: temp(nel)
      my_real,INTENT(INOUT) :: phin(nel) ! EINT/V
      my_real,INTENT(IN) :: rhocp
      my_real,INTENT(IN) :: phi(*), grad(6,nel), vol(nel)
      my_real,INTENT(IN) :: alpha(*)
      TYPE(t_ale_connectivity), INTENT(IN) :: ALE_CONNECT
C-----------------------------------------------
C   C o m m o n   B l o c k s
C-----------------------------------------------
#include      "com08_c.inc"
#include      "vect01_c.inc"
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I, IE, IV1, IV2, IV3, IV4, IV5, IV6, IAD2
      my_real dphi(mvsiz)
      my_real aa(0:6)     !< thermal diffusivity (updated from ALPHA for ghost cells)
      my_real aa_face(6)  ! k = alpha * rhocp
C-----------------------------------------------
C   S o u r c e   L i n e s
C-----------------------------------------------
      DO i=1,nel
        ie =nft+i
        iad2 = ale_connect%ee_connect%iad_connect(ie)
        iv1 = ale_connect%ee_connect%connected(iad2 + 1 - 1)
        iv2 = ale_connect%ee_connect%connected(iad2 + 2 - 1)
        iv3 = ale_connect%ee_connect%connected(iad2 + 3 - 1)
        iv4 = ale_connect%ee_connect%connected(iad2 + 4 - 1)
        iv5 = ale_connect%ee_connect%connected(iad2 + 5 - 1)
        iv6 = ale_connect%ee_connect%connected(iad2 + 6 - 1)
 
        ! adjacent cells (IV=0 => ghost cell with same value as IE)
        IF(iv1 <= 0)iv1=ie
        IF(iv2 <= 0)iv2=ie
        IF(iv3 <= 0)iv3=ie
        IF(iv4 <= 0)iv4=ie
        IF(iv5 <= 0)iv5=ie
        IF(iv6 <= 0)iv6=ie
 
        !thermal diffusivity ( *rho0cp)
        aa(0) = alpha(ie)
        aa(1) = alpha(iv1)
        aa(2) = alpha(iv2)
        aa(3) = alpha(iv3)
        aa(4) = alpha(iv4)
        aa(5) = alpha(iv5)
        aa(6) = alpha(iv6)
        !ghost cells
        IF(aa(1) == zero) aa(1)=aa(0)
        IF(aa(2) == zero) aa(2)=aa(0)
        IF(aa(3) == zero) aa(3)=aa(0)
        IF(aa(4) == zero) aa(4)=aa(0)
        IF(aa(5) == zero) aa(5)=aa(0)
        IF(aa(6) == zero) aa(6)=aa(0)
        !harmonic interpolation
        aa_face(1) = (aa(0)*aa(1)) / max(em20,(aa(0)+aa(1)))
        aa_face(2) = (aa(0)*aa(2)) / max(em20,(aa(0)+aa(2)))
        aa_face(3) = (aa(0)*aa(3)) / max(em20,(aa(0)+aa(3)))
        aa_face(4) = (aa(0)*aa(4)) / max(em20,(aa(0)+aa(4)))
        aa_face(5) = (aa(0)*aa(5)) / max(em20,(aa(0)+aa(5)))
        aa_face(6) = (aa(0)*aa(6)) / max(em20,(aa(0)+aa(6)))
C-----------------------------------------------------------
        ! time evolution
        dphi(i) = aa_face(1) * (phi(iv1)-phi(ie))*grad(1,i)
     2          + aa_face(2) * (phi(iv2)-phi(ie))*grad(2,i)
     3          + aa_face(3) * (phi(iv3)-phi(ie))*grad(3,i)
     4          + aa_face(4) * (phi(iv4)-phi(ie))*grad(4,i)
     5          + aa_face(5) * (phi(iv5)-phi(ie))*grad(5,i)
     6          + aa_face(6) * (phi(iv6)-phi(ie))*grad(6,i)
      enddo!next I
C-----------------------------------------------------------
      ! time integration for Eint/V
      ! %EINT is here Eint / V (J/m3)
      ! => DPHI = is finally m.cp.dT = dT *rhoCp   (ALPHA is k instead k/rhocp )
      DO i=1,nel
        dphi(i) = two*dphi(i)*dt1/max(vol(i),em20)
      ENDDO
C-----------------------------------------------------------
      ! Eint/V updated
      DO i=1,nel
        phin(i) = phin(i)+dphi(i)
      ENDDO
C-----------------------------------------------------------
      ! temperature updated
      IF(rhocp > zero)THEN
        DO i=1,nel
          temp(i) = temp(i) + dphi(i)/rhocp
        ENDDO
      ENDIF
C-----------------------------------------------------------
      RETURN
      END