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