OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
puff.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!|| puff ../common_source/eos/puff.F
25!||--- called by ------------------------------------------------------
26!|| eosmain ../common_source/eos/eosmain.F
27!||--- uses -----------------------------------------------------
28!|| constant_mod ../common_source/modules/constant_mod.F
29!||====================================================================
30 SUBROUTINE puff(IFLAG,NEL ,PM ,OFF ,EINT ,MU ,MU2,
31 2 ESPE ,DVOL ,DF ,VNEW ,MAT ,
32 3 PNEW ,DPDM ,DPDE )
33C-----------------------------------------------
34C D e s c r i p t i o n
35C-----------------------------------------------
36C This subroutine contains numerical solving
37C of PUFF EOS
38!----------------------------------------------------------------------------
39!! \details STAGGERED SCHEME IS EXECUTED IN TWO PASSES IN EOSMAIN : IFLG=0 THEN IFLG=1
40!! \details COLLOCATED SCHEME IS DOING A SINGLE PASS : IFLG=2
41!! \details
42!! \details STAGGERED SCHEME
43!! \details EOSMAIN / IFLG = 0 : DERIVATIVE CALCULATION FOR SOUND SPEED ESTIMATION c[n+1] REQUIRED FOR PSEUDO-VISCOSITY (DPDE:partial derivative, DPDM:total derivative)
44!! \details MQVISCB : PSEUDO-VISCOSITY Q[n+1]
45!! \details MEINT : INTERNAL ENERGY INTEGRATION FOR E[n+1] : FIRST PART USING P[n], Q[n], and Q[n+1] CONTRIBUTIONS
46!! \details EOSMAIN / IFLG = 1 : UPDATE P[n+1], T[N+1]
47!! \details INTERNAL ENERGY INTEGRATION FOR E[n+1] : LAST PART USING P[n+1] CONTRIBUTION
48!! \details (second order integration dE = -P.dV where P = 0.5(P[n+1] + P[n]) )
49!! \details COLLOCATED SCHEME
50!! \details EOSMAIN / IFLG = 2 : SINGLE PASS FOR P[n+1] AND DERIVATIVES
51!----------------------------------------------------------------------------
52C-----------------------------------------------
53C M o d u l e s
54C-----------------------------------------------
55 USE constant_mod , ONLY : zero, em15, half, one, three_half, two, three, three100
56C-----------------------------------------------
57C I m p l i c i t T y p e s
58C-----------------------------------------------
59 IMPLICIT NONE
60#include "my_real.inc"
61C-----------------------------------------------
62C C o m m o n B l o c k s
63C-----------------------------------------------
64#include "param_c.inc"
65#include "com04_c.inc"
66#include "com06_c.inc"
67#include "com08_c.inc"
68#include "vect01_c.inc"
69#include "scr06_c.inc"
70C-----------------------------------------------
71C D u m m y A r g u m e n t s
72C-----------------------------------------------
73 INTEGER MAT(NEL), IFLAG, NEL
74 my_real PM(NPROPM,NUMMAT),
75 . off(nel) , eint(nel), mu(nel) ,
76 . mu2(nel) , espe(nel), dvol(nel), df(nel) ,
77 . vnew(nel), pnew(nel), dpdm(nel), dpde(nel)
78C-----------------------------------------------
79C L o c a l V a r i a b l e s
80C-----------------------------------------------
81 INTEGER I, MX
82 my_real AA, BB, DVV, ETA, XX, GX, PRES, CC, EXPA, EE
83 my_real c1(nel),c2(nel),c3(nel),t1(nel),t2(nel), g0(nel),esubl(nel),hh(nel),pc(nel),psh(nel)
84C--------------------------------------------------------------------
85 IF(iflag == 0) THEN
86 DO i=1,nel
87 mx = mat(i)
88 c1(i) = pm(32,mx)
89 c2(i) = pm(33,mx)
90 c3(i) = pm(34,mx)
91 t1(i) = pm(35,mx)
92 t2(i) = pm(36,mx)
93 pc(i) = pm(37,mx)
94 esubl(i)= pm(160,mx)
95 g0(i) = pm(161,mx)
96 hh(i) = pm(162,mx)
97 psh(i) = pm(88,mx)
98 ENDDO
99 DO i=1,nel
100 xx =mu(i)/(one+mu(i))
101 IF(mu(i) >= zero) THEN
102 aa=(c1(i)+c3(i)*mu2(i))*mu(i)+c2(i)*mu2(i)
103 gx=one-g0(i)*half*xx
104 bb=g0(i)
105 pres=max(aa*gx+bb*espe(i),pc(i))*off(i)
106 dpdm(i)=(c1(i)+two*c2(i)*mu(i)+three*c3(i)*mu2(i))*gx+g0(i)*df(i)*df(i)*(pres-half*aa)
107 ELSEIF(espe(i) < esubl(i))THEN
108 aa=(t1(i)+t2(i)*mu(i))*mu(i)
109 gx=one-g0(i)*half*xx
110 bb=g0(i)
111 pres=max(aa*gx+bb*espe(i),pc(i))*off(i)
112 dpdm(i)=(t1(i)+two*t2(i)*mu(i))*gx+g0(i)*df(i)*df(i)*(pres-half*aa)
113 ELSE
114 eta=one+mu(i)
115 ee=sqrt(eta)
116 bb=(hh(i)+(g0(i)-hh(i))*ee)*eta
117 cc= c1(i)/(g0(i)*esubl(i))
118 expa=exp(cc*xx)
119 aa= bb*esubl(i)*(expa-one)
120 pres=max(aa+bb*espe(i),pc(i))*off(i)
121 dpdm(i)=bb*df(i)*df(i)*(pres+esubl(i)*expa*cc) +(espe(i)+esubl(i)*(expa-one))*(hh(i)+three_half*ee*(g0(i)-hh(i)))
122 ENDIF
123 dpde(i)=bb
124 pnew(i) = max(pres,pc(i))*off(i)! P(mu[n+1],E[n])
125 ENDDO
126
127 ELSEIF(iflag == 1) THEN
128 DO i=1,nel
129 mx = mat(i)
130 c1(i) = pm(32,mx)
131 c2(i) = pm(33,mx)
132 c3(i) = pm(34,mx)
133 t1(i) = pm(35,mx)
134 t2(i) = pm(36,mx)
135 pc(i) = pm(37,mx)
136 esubl(i)= pm(160,mx)
137 g0(i) = pm(161,mx)
138 hh(i) = pm(162,mx)
139 psh(i) = pm( 88,mx)
140 ENDDO
141 DO i=1,nel
142 dvv=half*dvol(i)*df(i) / max(em15,vnew(i))
143 xx =mu(i)/(one+mu(i))
144 IF(mu(i) >= zero) THEN
145 aa=(c1(i)+c3(i)*mu2(i))*mu(i)+c2(i)*mu2(i)
146 aa=aa*(one-g0(i)*half*xx)
147 bb=g0(i)
148 ELSEIF(espe(i) < esubl(i))THEN
149 aa=(t1(i)+t2(i)*mu(i))*mu(i)
150 aa=aa*(one-g0(i)*half*xx)
151 bb=g0(i)
152 ELSE
153 eta=one+mu(i)
154 bb=(hh(i)+(g0(i)-hh(i))*sqrt(eta))*eta
155 cc=c1(i)/(g0(i)*esubl(i))
156 expa=exp(cc*xx)
157 aa=bb*esubl(i)*(expa-one)
158 ENDIF
159 dpde(i)=bb
160 pnew(i)=(aa +bb*espe(i))/(one+ bb*dvv)
161 pnew(i)=max(pnew(i),pc(i))*off(i)! P(mu[n+1],E[n+1])
162 eint(i)=eint(i)-half*dvol(i)*pnew(i)
163 ENDDO
164
165 ELSEIF(iflag == 2) THEN
166 DO i=1, nel
167 mx = mat(i)
168 c1(i) = pm(32,mx)
169 c2(i) = pm(33,mx)
170 c3(i) = pm(34,mx)
171 t1(i) = pm(35,mx)
172 t2(i) = pm(36,mx)
173 pc(i) = pm(37,mx)
174 esubl(i)= pm(160,mx)
175 g0(i) = pm(161,mx)
176 hh(i) = pm(162,mx)
177 ENDDO
178 DO i=1, nel
179 IF (vnew(i) > zero) THEN
180 xx =mu(i)/(one+mu(i))
181 IF(mu(i) >= zero) THEN
182 aa=(c1(i)+c3(i)*mu2(i))*mu(i)+c2(i)*mu2(i)
183 gx=one-g0(i)*half*xx
184 bb=g0(i)
185 pres=max(aa*gx+bb*espe(i),pc(i))*off(i)
186 dpdm(i)=(c1(i)+two*c2(i)*mu(i)+three*c3(i)*mu2(i))*gx+g0(i)*df(i)*df(i)*(pres-half*aa)
187 ELSEIF(espe(i)<esubl(i))THEN
188 aa=(t1(i)+t2(i)*mu(i))*mu(i)
189 gx=one-g0(i)*half*xx
190 bb=g0(i)
191 pres=max(aa*gx+bb*espe(i),pc(i))*off(i)
192 dpdm(i)=(t1(i)+two*t2(i)*mu(i))*gx+g0(i)*df(i)*df(i)*(pres-half*aa)
193 ELSE
194 eta=one+mu(i)
195 ee=sqrt(eta)
196 bb=(hh(i)+(g0(i)-hh(i))*ee)*eta
197 cc= c1(i)/(g0(i)*esubl(i))
198 expa=exp(cc*xx)
199 aa= bb*esubl(i)*(expa-one)
200 pres=max(aa+bb*espe(i),pc(i))*off(i)
201 dpdm(i)=bb*df(i)*df(i)*(pres+esubl(i)*expa*cc) +
202 . (espe(i)+esubl(i)*(expa-one))*(hh(i)+three_half*ee*(g0(i)-hh(i)))
203 ENDIF
204 dpde(i)=bb
205 pnew(i) = pres
206 ENDIF
207 ENDDO
208 ENDIF
209C-----------------------------------------------
210 RETURN
211 END
max
#define max(a, b)
Definition macros.h:21
constant_mod
Definition constant_mod.F:294
puff
subroutine puff(iflag, nel, pm, off, eint, mu, mu2, espe, dvol, df, vnew, mat, pnew, dpdm, dpde)
Definition puff.F:33