OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
idealgas_vt_mod Module Reference

Functions/Subroutines

subroutine idealgas_vt (iflag, nel, pmin, off, eint, mu, mu2, espe, dvol, df, vnew, psh, rho0, rho, pnew, dpdm, dpde, theta, pold, eos_struct)

Function/Subroutine Documentation

◆ idealgas_vt()

subroutine idealgas_vt_mod::idealgas_vt ( integer, intent(in) iflag,
integer, intent(in) nel,
intent(in) pmin,
dimension(nel), intent(in) off,
dimension(nel), intent(inout) eint,
dimension(nel), intent(in) mu,
dimension(nel), intent(in) mu2,
dimension(nel), intent(in) espe,
dimension(nel), intent(in) dvol,
dimension(nel), intent(in) df,
dimension(nel), intent(in) vnew,
dimension(nel), intent(inout) psh,
dimension(nel), intent(in) rho0,
dimension(nel), intent(in) rho,
dimension(nel), intent(inout) pnew,
dimension(nel), intent(inout) dpdm,
dimension(nel), intent(inout) dpde,
dimension(nel), intent(inout) theta,
dimension(nel), intent(in) pold,
type(eos_param_), intent(in) eos_struct )

Definition at line 41 of file idealgas_vt.F.

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
my_real
#define my_real
Definition cppsort.cpp:32
constant_mod
Definition constant_mod.F:355