OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
gruneisen.F
Go to the documentation of this file.
1Copyright> OpenRadioss
2Copyright> Copyright (C) 1986-2025 Altair Engineering Inc.
3Copyright>
4Copyright> This program is free software: you can redistribute it and/or modify
5Copyright> it under the terms of the GNU Affero General Public License as published by
6Copyright> the Free Software Foundation, either version 3 of the License, or
7Copyright> (at your option) any later version.
8Copyright>
9Copyright> This program is distributed in the hope that it will be useful,
10Copyright> but WITHOUT ANY WARRANTY; without even the implied warranty of
11Copyright> MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12Copyright> GNU Affero General Public License for more details.
13Copyright>
14Copyright> You should have received a copy of the GNU Affero General Public License
15Copyright> along with this program. If not, see <https://www.gnu.org/licenses/>.
16Copyright>
17Copyright>
18Copyright> Commercial Alternative: Altair Radioss Software
19Copyright>
20Copyright> As an alternative to this open-source version, Altair also offers Altair Radioss
21Copyright> software under a commercial license. Contact Altair to discuss further if the
22Copyright> commercial version may interest you: https://www.altair.com/radioss/.
23!||====================================================================
24!|| gruneisen_mod ../common_source/eos/gruneisen.F
25!||--- called by ------------------------------------------------------
26!|| eosmain ../common_source/eos/eosmain.F
27!|| eosmain51 ../engine/source/materials/mat/mat051/eosmain51.F90
28!||====================================================================
29 MODULE gruneisen_mod
30 IMPLICIT NONE
31 CONTAINS
32!||====================================================================
33!|| gruneisen ../common_source/eos/gruneisen.F
34!||--- called by ------------------------------------------------------
35!|| eosmain ../common_source/eos/eosmain.F
36!|| eosmain51 ../engine/source/materials/mat/mat051/eosmain51.F90
37!||--- uses -----------------------------------------------------
38!|| eos_param_mod ../common_source/modules/mat_elem/eos_param_mod.F90
39!|| precision_mod ../common_source/modules/precision_mod.F90
40!||====================================================================
41 SUBROUTINE gruneisen(
42 1 IFLAG , NEL ,PMIN ,OFF ,EINT ,MU ,MU2 ,
43 2 ESPE , DVOL ,DF ,VNEW ,RHO0,
44 3 PNEW , DPDM ,DPDE ,PSH ,
45 4 EOS_STRUCT)
46C-----------------------------------------------
47C M o d u l e s
48C-----------------------------------------------
49 USE eos_param_mod , ONLY : eos_param_
50 USE precision_mod, ONLY : wp
51C-----------------------------------------------
52C D e s c r i p t i o n
53C-----------------------------------------------
54!----------------------------------------------------------------------------
55!! \details STAGGERED SCHEME IS EXECUTED IN TWO PASSES IN EOSMAIN : IFLG=0 THEN IFLG=1
56!! \details COLLOCATED SCHEME IS DOING A SINGLE PASS : IFLG=2
57!! \details
58!! \details STAGGERED SCHEME
59!! \details EOSMAIN / IFLG = 0 : DERIVATIVE CALCULATION FOR SOUND SPEED ESTIMATION c[n+1] REQUIRED FOR PSEUDO-VISCOSITY (DPDE:partial derivative, DPDM:total derivative)
60!! \details MQVISCB : PSEUDO-VISCOSITY Q[n+1]
61!! \details MEINT : INTERNAL ENERGY INTEGRATION FOR E[n+1] : FIRST PART USING P[n], Q[n], and Q[n+1] CONTRIBUTIONS
62!! \details EOSMAIN / IFLG = 1 : UPDATE P[n+1], T[N+1]
63!! \details INTERNAL ENERGY INTEGRATION FOR E[n+1] : LAST PART USING P[n+1] CONTRIBUTION
64!! \details (second order integration dE = -P.dV where P = 0.5(P[n+1] + P[n]) )
65!! \details COLLOCATED SCHEME
66!! \details EOSMAIN / IFLG = 2 : SINGLE PASS FOR P[n+1] AND DERIVATIVES
67!----------------------------------------------------------------------------
68C-----------------------------------------------
69C I m p l i c i t T y p e s
70C-----------------------------------------------
71#include "implicit_f.inc"
72C-----------------------------------------------
73C C o m m o n B l o c k s
74C-----------------------------------------------
75C
76C-----------------------------------------------
77C D u m m y A r g u m e n t s
78C-----------------------------------------------
79 INTEGER,INTENT(IN) :: IFLAG !< flag for the operation (staggered scheme IFLAG=0,1 ; collocated scheme IFLAG=2)
80 INTEGER,INTENT(IN) :: NEL !< number of elements in the current group
81 my_real, INTENT(IN) :: PMIN !< minimum pressure allowed
82 my_real,INTENT(IN) :: OFF(NEL) !< OFF=1 if element is active, OFF=0 if element is inactive
83 my_real,INTENT(IN) :: mu(nel) !< volumetric strain : rho/rho0-1
84 my_real,INTENT(IN) :: mu2(nel) !< MU**2 IF MU > 0 , 0 otherwise
85 my_real,INTENT(IN) :: espe(nel) !< specific internal energy
86 my_real,INTENT(IN) :: dvol(nel) !< volume change
87 my_real,INTENT(IN) :: vnew(nel) !< current volume
88 my_real,INTENT(IN) :: rho0(nel) !< reference density
89 my_real,INTENT(IN) :: df(nel) !< relative volume (v/v0 = rho0/rho)
90 my_real,INTENT(INOUT) :: eint(nel) !< internal energy
91 my_real,INTENT(INOUT) :: psh(nel) !< pressure shift (for relative pressure modeling)
92 my_real,INTENT(INOUT) :: pnew(nel) !< current pressure
93 my_real,INTENT(INOUT) :: dpde(nel) !< partial derivative dP/dE where E=Eint/V0
94 my_real,INTENT(INOUT) :: dpdm(nel) !< total derivative : DP/Dmu = dPdmu + dPdE* P/(1+mu)**2
95 TYPE(eos_param_),INTENT(IN) :: EOS_STRUCT !< data structure for EoS parameters
96C-----------------------------------------------
97C L o c a l V a r i a b l e s
98C-----------------------------------------------
99 INTEGER I
100 real(kind=WP) :: AA, BB, DVV, FF, FG, FAC, XX, DFF, DFG, FAC1, PP
101 real(kind=WP) :: CC,S1,S2,S3,G0,GA
102C-----------------------------------------------
103C B o d y
104C-----------------------------------------------
105 cc = eos_struct%UPARAM(1)
106 s1 = eos_struct%UPARAM(2)
107 s2 = eos_struct%UPARAM(3)
108 s3 = eos_struct%UPARAM(4)
109 g0 = eos_struct%UPARAM(5)
110 ga = eos_struct%UPARAM(6)
111 psh(:) = eos_struct%PSH
112
113 IF(iflag == 0) THEN !--- SOUND SPEED DERIVATIVE (DPDE is partial derivative, DPDM is total derivative)
114 ! Derivative
115 DO i=1,nel
116 fac=one
117 fac1=one
118 IF(mu(i) > zero) THEN
119 xx= mu(i)/(one+mu(i))
120 ff=one+(one-half*g0)*mu(i)-half*ga*mu2(i)
121 fg=one-(s1-one+s2*xx+s3*xx*xx)*mu(i)
122 fac=ff/(fg*fg)
123 dff=one-half*g0-ga*mu(i)
124 dfg=one-s1+xx*(-two*s2+xx*(s2-three*s3)+two*s3*xx*xx)
125 fac1=fac*(one+mu(i)*(dff/ff-two*dfg/fg))
126 ENDIF
127 aa=fac*rho0(i)*cc*cc*mu(i)
128 bb=g0+ga*mu(i)
129 pp=max(aa+bb*espe(i),pmin)*off(i) ! total pressure need for sound speed calculation (=> no pressure shift)
130 dpdm(i)=fac1*rho0(i)*cc*cc+pp*df(i)*df(i)*bb+ga*espe(i)
131 dpde(i)=bb
132 pnew(i) = max(pp,pmin)*off(i)! P(mu[n+1],E[n])
133 pnew(i) = pnew(i) - psh(i)
134 ENDDO
135
136 ELSEIF(iflag == 1) THEN !--- P[n+1] and EINT[n+1]
137 ! Pressure P[n+1]
138 ! Internal Energy Eint[n+1] (2nd order integration requires P[n+1])
139 DO i=1,nel
140 fac=one
141 IF(mu(i) > zero) THEN
142 xx= mu(i)/(one+mu(i))
143 ff=one+(one-half*g0)*mu(i)-half*ga*mu2(i)
144 fg=one-(s1-one+s2*xx+s3*xx*xx)*mu(i)
145 fac=ff/(fg*fg)
146 ENDIF
147 aa=fac*rho0(i)*cc*cc*mu(i)
148 bb=g0+ga*mu(i)
149 dpde(i) = bb
150 dvv=half*dvol(i)*df(i) / max(em15,vnew(i))
151 pnew(i)=(aa+bb*espe(i))/(one+bb*dvv)
152 pnew(i)= max(pnew(i),pmin)*off(i) ! P(mu[n+1],E[n+1])
153 eint(i)=eint(i) - half*dvol(i)*pnew(i) !total pressure for energy integration
154 pnew(i) = pnew(i) - psh(i)
155 ENDDO
156
157 ELSEIF(iflag == 2) THEN !--- collocated scheme (law151)
158 ! pressure and sound speed derivatives
159 DO i=1,nel
160 IF (vnew(i) > zero) THEN
161 fac=one
162 fac1=one
163 IF(mu(i) > zero) THEN
164 xx= mu(i)/(one+mu(i))
165 ff=one+(one-half*g0)*mu(i)-half*ga*mu2(i)
166 fg=one-(s1-one+s2*xx+s3*xx*xx)*mu(i)
167 fac=ff/(fg*fg)
168 dff=one-half*g0-ga*mu(i)
169 dfg=one-s1+xx*(-two*s2+xx*(s2-three*s3)+two*s3*xx*xx)
170 fac1=fac*(one+mu(i)*(dff/ff-two*dfg/fg))
171 ENDIF
172 aa=fac*rho0(i)*cc*cc*mu(i)
173 bb=g0+ga*mu(i)
174 pnew(i)=max(aa+bb*espe(i),pmin)*off(i)
175 dpdm(i)=fac1*rho0(i)*cc*cc+pnew(i)*df(i)*df(i)*bb+ga*espe(i)
176 dpde(i)=bb
177 pnew(i)=pnew(i)-psh(i)
178 ENDIF
179 ENDDO
180
181 ENDIF
182
183 RETURN
184 END SUBROUTINE gruneisen
185 END MODULE gruneisen_mod
max
#define max(a, b)
Definition macros.h:21
gruneisen_mod
Definition gruneisen.F:29
gruneisen_mod::gruneisen
subroutine gruneisen(iflag, nel, pmin, off, eint, mu, mu2, espe, dvol, df, vnew, rho0, pnew, dpdm, dpde, psh, eos_struct)
Definition gruneisen.F:46