OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
idealgas_vt.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!|| idealgas_vt_mod ../common_source/eos/idealgas_vt.F
25!||--- called by ------------------------------------------------------
26!|| eosmain ../common_source/eos/eosmain.F
27!|| eosmain51 ../engine/source/materials/mat/mat051/eosmain51.F90
28!||====================================================================
29 MODULE idealgas_vt_mod
30 IMPLICIT NONE
31 CONTAINS
32!||====================================================================
33!|| idealgas_vt ../common_source/eos/idealgas_vt.F
34!||--- called by ------------------------------------------------------
35!|| eosmain ../common_source/eos/eosmain.F
36!|| eosmain51 ../engine/source/materials/mat/mat051/eosmain51.F90
37!||--- uses -----------------------------------------------------
38!|| constant_mod ../common_source/modules/constant_mod.F
39!|| eos_param_mod ../common_source/modules/mat_elem/eos_param_mod.F90
40!||====================================================================
41 SUBROUTINE idealgas_vt(
42 1 IFLAG ,NEL ,PMIN ,OFF ,EINT ,MU ,MU2 ,
43 2 ESPE ,DVOL ,DF ,VNEW ,PSH ,RHO0, RHO,
44 3 PNEW ,DPDM ,DPDE ,THETA,POLD ,EOS_STRUCT)
45C-----------------------------------------------
46C M o d u l e s
47C-----------------------------------------------
48 USE constant_mod , ONLY : zero, em05, half, one, two, third, fourth, one_fifth
49 USE eos_param_mod , ONLY : eos_param_
50C-----------------------------------------------
51C D e s c r i p t i o n
52C-----------------------------------------------
53C This subroutine contains numerical solving
54C of IDEAL GAS EOS with variable Cp(T) parameter
55!----------------------------------------------------------------------------
56!! \details STAGGERED SCHEME IS EXECUTED IN TWO PASSES IN EOSMAIN : IFLG=0 THEN IFLG=1
57!! \details COLLOCATED SCHEME IS DOING A SINGLE PASS : IFLG=2
58!! \details
59!! \details STAGGERED SCHEME
60!! \details EOSMAIN / IFLG = 0 : DERIVATIVE CALCULATION FOR SOUND SPEED ESTIMATION c[n+1] REQUIRED FOR PSEUDO-VISCOSITY (DPDE:partial derivative, DPDM:total derivative)
61!! \details MQVISCB : PSEUDO-VISCOSITY Q[n+1]
62!! \details MEINT : INTERNAL ENERGY INTEGRATION FOR E[n+1] : FIRST PART USING P[n], Q[n], and Q[n+1] CONTRIBUTIONS
63!! \details EOSMAIN / IFLG = 1 : UPDATE P[n+1], T[N+1]
64!! \details INTERNAL ENERGY INTEGRATION FOR E[n+1] : LAST PART USING P[n+1] CONTRIBUTION
65!! \details (second order integration dE = -P.dV where P = 0.5(P[n+1] + P[n]) )
66!! \details COLLOCATED SCHEME
67!! \details EOSMAIN / IFLG = 2 : SINGLE PASS FOR P[n+1] AND DERIVATIVES
68!----------------------------------------------------------------------------
69C-----------------------------------------------
70C I m p l i c i t T y p e s
71C-----------------------------------------------
72 IMPLICIT NONE
73#include "my_real.inc"
74C-----------------------------------------------
75C D u m m y A r g u m e n t s
76C-----------------------------------------------
77 INTEGER,INTENT(IN) :: IFLAG, NEL
78 my_real,INTENT(IN) :: PMIN,OFF(NEL),MU(NEL)
79 my_real,INTENT(IN) :: MU2(NEL) ,ESPE(NEL) ,DVOL(NEL) ,DF(NEL), VNEW(NEL)
80 my_real, INTENT(INOUT) :: psh(nel),pnew(nel),dpdm(nel),dpde(nel),eint(nel),theta(nel)
81 my_real,INTENT(IN) :: pold(nel), rho0(nel), rho(nel)
82 TYPE(eos_param_),INTENT(IN) :: EOS_STRUCT
83C-----------------------------------------------
84C L o c a l V a r i a b l e s
85C-----------------------------------------------
86 INTEGER I
87 my_real :: P0,GAMMA
88 my_real :: r_GAS,A0,A1,A2,A3,A4,CP(NEL),CV
89 my_real :: temp, fun, dfun, tol, error, incr
90 INTEGER :: ITER, MAX_ITER
91C-----------------------------------------------
92C S o u r c e L i n e s
93C-----------------------------------------------
94 max_iter = 50
95 tol = em05
96
97 cp(1:nel)=zero
98 cv = eos_struct%CV
99 psh = eos_struct%PSH
100 a0 = eos_struct%UPARAM(1)
101 a1 = eos_struct%UPARAM(2)
102 a2 = eos_struct%UPARAM(3)
103 a3 = eos_struct%UPARAM(4)
104 a4 = eos_struct%UPARAM(5)
105 r_gas = eos_struct%UPARAM(6)
106 p0 = eos_struct%UPARAM(7)
107
108 IF(iflag == 0) THEN
109 !SOLVE TEMPERATURE
110 DO i=1,nel
111 ! Init newton
112 temp = pold(i)/rho(i)/r_gas
113 iter = 0
114 error = huge(zero)
115 DO WHILE(error > tol .AND. iter < max_iter)
116 ! f(x) = 0 <=> int(cv(T), T) - eint = 0
117 ! <=> int(cp(T) - r_gas, T) - eint = 0
118 fun = a0 * temp + half * a1 * temp**2 + third * a2 * temp**3 +
119 . fourth * a3 * temp**4 + one_fifth * a4 * temp**5 - r_gas * temp - espe(i) / rho0(i)
120 IF (abs(fun) < tol) EXIT
121 dfun = a0 + a1 * temp + a2 * temp**2 + a3 * temp**3 + a4 * temp**4 - r_gas
122 incr = - fun / dfun
123 temp = temp + incr
124 error = abs(incr / temp)
125 iter = iter + 1
126 ENDDO
127 ! Store
128 theta(i) = temp
129 cp(i) = a0 + a1 * temp + a2 * temp**2 + a3 * temp**3 + a4 * temp**4
130 ENDDO
131 DO i=1,nel
132 cv = cp(i) - r_gas
133 gamma = cp(i)/cv
134 dpdm(i) = rho0(i)*gamma*r_gas*theta(i) !total derivative
135 dpde(i) = gamma*(one+mu(i))
136 pnew(i) = rho0(i)*(one+mu(i))*r_gas*theta(i) ! P(mu[n+1],E[n])
137 pnew(i) = pnew(i)-psh(i)
138 ENDDO
139
140C-----------------------------------------------
141 ELSEIF(iflag == 1) THEN
142 !SOLVE TEMPERATURE
143 DO i=1,nel
144 ! Init newton
145 temp = pold(i)/rho(i)/r_gas
146 iter = 0
147 error = huge(zero)
148 DO WHILE(error > tol .AND. iter < max_iter)
149 ! f(x) = 0 <=> int(cv(T), T) - eint = 0
150 ! <=> int(cp(T) - r_gas, T) - eint = 0
151 fun = a0 * temp + half * a1 * temp**2 + third * a2 * temp**3 +
152 . fourth * a3 * temp**4 + one_fifth * a4 * temp**5 - r_gas * temp - espe(i) / rho0(i)
153 IF (abs(fun) < tol) EXIT
154 dfun = a0 + a1 * temp + a2 * temp**2 + a3 * temp**3 + a4 * temp**4 - r_gas
155 incr = - fun / dfun
156 temp = temp + incr
157 error = abs(incr / temp)
158 iter = iter + 1
159 ENDDO
160 ! Store
161 theta(i) = temp
162 cp(i) = a0 + a1 * temp + a2 * temp**2 + a3 * temp**3 + a4 * temp**4
163 ENDDO
164 DO i=1,nel
165 pnew(i) = rho0(i)*(one+mu(i))*r_gas*theta(i) ! P(mu[n+1],E[n+1])
166 eint(i) = eint(i) - half*dvol(i)*pnew(i)
167 pnew(i) = pnew(i)-psh(i)
168 cv = cp(i) - r_gas
169 gamma = cp(i)/cv
170 dpde(i) = gamma*(one+mu(i))
171 ENDDO
172
173 ELSEIF (iflag == 2) THEN
174 !SOLVE TEMPERATURE
175 DO i=1,nel
176 ! Init newton
177 temp = pold(i)/rho(i)/r_gas
178 iter = 0
179 error = huge(zero)
180 DO WHILE(error > tol .AND. iter < max_iter)
181 ! f(x) = 0 <=> int(cv(T), T) - eint = 0
182 ! <=> int(cp(T) - r_gas, T) - eint = 0
183 fun = a0 * temp + half * a1 * temp**2 + third * a2 * temp**3 +
184 . fourth * a3 * temp**4 + one_fifth * a4 * temp**5 - r_gas * temp - espe(i) / rho0(i)
185 IF (abs(fun) < tol) EXIT
186 dfun = a0 + a1 * temp + a2 * temp**2 + a3 * temp**3 + a4 * temp**4 - r_gas
187 incr = - fun / dfun
188 temp = temp + incr
189 error = abs(incr / temp)
190 iter = iter + 1
191 ENDDO
192 ! Store
193 theta(i) = temp
194 cp(i) = a0 + a1 * temp + a2 * temp**2 + a3 * temp**3 + a4 * temp**4
195 ENDDO
196 DO i=1, nel
197 IF (vnew(i) > zero) THEN
198 cv = cp(i) - r_gas
199 gamma = cp(i)/cv
200 dpdm(i) = rho0(i)*gamma*r_gas*theta(i) !total derivative
201 dpde(i) = gamma*(one+mu(i)) !partial derivative
202 pnew(i) = rho0(i)*(one+mu(i))*r_gas*theta(i)
203 ENDIF
204 ENDDO
205
206 ENDIF
207C-----------------------------------------------
208 RETURN
209 END SUBROUTINE idealgas_vt
210 END MODULE idealgas_vt_mod
constant_mod
Definition constant_mod.F:355
idealgas_vt_mod
Definition idealgas_vt.F:29
idealgas_vt_mod::idealgas_vt
subroutine idealgas_vt(iflag, nel, pmin, off, eint, mu, mu2, espe, dvol, df, vnew, psh, rho0, rho, pnew, dpdm, dpde, theta, pold, eos_struct)
Definition idealgas_vt.F:45