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_mod ../common_source/eos/puff.F
25!||--- called by ------------------------------------------------------
26!|| eosmain ../common_source/eos/eosmain.F
27!|| eosmain51 ../engine/source/materials/mat/mat051/eosmain51.F90
28!||====================================================================
29 MODULE puff_mod
30 IMPLICIT NONE
31 CONTAINS
32!||====================================================================
33!|| puff ../common_source/eos/puff.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 puff(
42 1 IFLAG ,NEL ,PMIN ,OFF ,EINT ,MU ,MU2,
43 2 ESPE ,DVOL ,DF ,VNEW ,
44 3 PNEW ,DPDM ,DPDE ,EOS_STRUCT)
45C-----------------------------------------------
46C D e s c r i p t i o n
47C-----------------------------------------------
48C This subroutine contains numerical solving
49C of PUFF EOS
50!----------------------------------------------------------------------------
51!! \details STAGGERED SCHEME IS EXECUTED IN TWO PASSES IN EOSMAIN : IFLG=0 THEN IFLG=1
52!! \details COLLOCATED SCHEME IS DOING A SINGLE PASS : IFLG=2
53!! \details
54!! \details STAGGERED SCHEME
55!! \details EOSMAIN / IFLG = 0 : DERIVATIVE CALCULATION FOR SOUND SPEED ESTIMATION c[n+1] REQUIRED FOR PSEUDO-VISCOSITY (DPDE:partial derivative, DPDM:total derivative)
56!! \details MQVISCB : PSEUDO-VISCOSITY Q[n+1]
57!! \details MEINT : INTERNAL ENERGY INTEGRATION FOR E[n+1] : FIRST PART USING P[n], Q[n], and Q[n+1] CONTRIBUTIONS
58!! \details EOSMAIN / IFLG = 1 : UPDATE P[n+1], T[N+1]
59!! \details INTERNAL ENERGY INTEGRATION FOR E[n+1] : LAST PART USING P[n+1] CONTRIBUTION
60!! \details (second order integration dE = -P.dV where P = 0.5(P[n+1] + P[n]) )
61!! \details COLLOCATED SCHEME
62!! \details EOSMAIN / IFLG = 2 : SINGLE PASS FOR P[n+1] AND DERIVATIVES
63!----------------------------------------------------------------------------
64C-----------------------------------------------
65C M o d u l e s
66C-----------------------------------------------
67 USE constant_mod , ONLY : zero, em15, half, one, three_half, two, three, three100
68 USE eos_param_mod , ONLY : eos_param_
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 C o m m o n B l o c k s
76C-----------------------------------------------
77C
78C-----------------------------------------------
79C D u m m y A r g u m e n t s
80C-----------------------------------------------
81 INTEGER,INTENT(IN) :: IFLAG, NEL
82 my_real,INTENT(IN) :: PMIN, OFF(NEL), MU(NEL),DVOL(NEL),MU2(NEL),DF(NEL),VNEW(NEL),ESPE(NEL)
83 my_real,INTENT(INOUT) :: DPDM(NEL), DPDE(NEL), PNEW(NEL)
84 my_real,INTENT(INOUT) :: eint(nel)
85 TYPE(eos_param_) , INTENT(IN) :: EOS_STRUCT
86C-----------------------------------------------
87C L o c a l V a r i a b l e s
88C-----------------------------------------------
89 INTEGER I
90 my_real AA, BB, DVV, ETA, XX, GX, PRES, CC, EXPA, EE
91 my_real C1,C2,C3,T1,T2, G0,ESUBL,HH,PSH(NEL)
92C--------------------------------------------------------------------
93
94 c1 = eos_struct%UPARAM(1)
95 c2 = eos_struct%UPARAM(2)
96 c3 = eos_struct%UPARAM(3)
97 t1 = eos_struct%UPARAM(4)
98 t2 = eos_struct%UPARAM(5)
99 esubl = eos_struct%UPARAM(6)
100 g0 = eos_struct%UPARAM(7)
101 hh = eos_struct%UPARAM(8)
102
103 psh(1:nel) = eos_struct%PSH
104
105
106 IF(iflag == 0) THEN
107 DO i=1,nel
108 xx =mu(i)/(one+mu(i))
109 IF(mu(i) >= zero) THEN
110 aa=(c1+c3*mu2(i))*mu(i)+c2*mu2(i)
111 gx=one-g0*half*xx
112 bb=g0
113 pres=max(aa*gx+bb*espe(i),pmin)*off(i)
114 dpdm(i)=(c1+two*c2*mu(i)+three*c3*mu2(i))*gx+g0*df(i)*df(i)*(pres-half*aa)
115 ELSEIF(espe(i) < esubl)THEN
116 aa=(t1+t2*mu(i))*mu(i)
117 gx=one-g0*half*xx
118 bb=g0
119 pres=max(aa*gx+bb*espe(i),pmin)*off(i)
120 dpdm(i)=(t1+two*t2*mu(i))*gx+g0*df(i)*df(i)*(pres-half*aa)
121 ELSE
122 eta=one+mu(i)
123 ee=sqrt(eta)
124 bb=(hh+(g0-hh)*ee)*eta
125 cc= c1/(g0*esubl)
126 expa=exp(cc*xx)
127 aa= bb*esubl*(expa-one)
128 pres=max(aa+bb*espe(i),pmin)*off(i)
129 dpdm(i)=bb*df(i)*df(i)*(pres+esubl*expa*cc) +(espe(i)+esubl*(expa-one))*(hh+three_half*ee*(g0-hh))
130 ENDIF
131 dpde(i)=bb
132 pnew(i) = max(pres,pmin)*off(i)! P(mu[n+1],E[n])
133 ENDDO
134
135 ELSEIF(iflag == 1) THEN
136 DO i=1,nel
137 dvv=half*dvol(i)*df(i) / max(em15,vnew(i))
138 xx =mu(i)/(one+mu(i))
139 IF(mu(i) >= zero) THEN
140 aa=(c1+c3*mu2(i))*mu(i)+c2*mu2(i)
141 aa=aa*(one-g0*half*xx)
142 bb=g0
143 ELSEIF(espe(i) < esubl)THEN
144 aa=(t1+t2*mu(i))*mu(i)
145 aa=aa*(one-g0*half*xx)
146 bb=g0
147 ELSE
148 eta=one+mu(i)
149 bb=(hh+(g0-hh)*sqrt(eta))*eta
150 cc=c1/(g0*esubl)
151 expa=exp(cc*xx)
152 aa=bb*esubl*(expa-one)
153 ENDIF
154 dpde(i)=bb
155 pnew(i)=(aa +bb*espe(i))/(one+ bb*dvv)
156 pnew(i)=max(pnew(i),pmin)*off(i)! P(mu[n+1],E[n+1])
157 eint(i)=eint(i)-half*dvol(i)*pnew(i)
158 ENDDO
159
160 ELSEIF(iflag == 2) THEN
161 DO i=1, nel
162 IF (vnew(i) > zero) THEN
163 xx =mu(i)/(one+mu(i))
164 IF(mu(i) >= zero) THEN
165 aa=(c1+c3*mu2(i))*mu(i)+c2*mu2(i)
166 gx=one-g0*half*xx
167 bb=g0
168 pres=max(aa*gx+bb*espe(i),pmin)*off(i)
169 dpdm(i)=(c1+two*c2*mu(i)+three*c3*mu2(i))*gx+g0*df(i)*df(i)*(pres-half*aa)
170 ELSEIF(espe(i)<esubl)THEN
171 aa=(t1+t2*mu(i))*mu(i)
172 gx=one-g0*half*xx
173 bb=g0
174 pres=max(aa*gx+bb*espe(i),pmin)*off(i)
175 dpdm(i)=(t1+two*t2*mu(i))*gx+g0*df(i)*df(i)*(pres-half*aa)
176 ELSE
177 eta=one+mu(i)
178 ee=sqrt(eta)
179 bb=(hh+(g0-hh)*ee)*eta
180 cc= c1/(g0*esubl)
181 expa=exp(cc*xx)
182 aa= bb*esubl*(expa-one)
183 pres=max(aa+bb*espe(i),pmin)*off(i)
184 dpdm(i)=bb*df(i)*df(i)*(pres+esubl*expa*cc) + (espe(i)+esubl*(expa-one))*(hh+three_half*ee*(g0-hh))
185 ENDIF
186 dpde(i)=bb
187 pnew(i) = pres
188 ENDIF
189 ENDDO
190 ENDIF
191C-----------------------------------------------
192 RETURN
193 END SUBROUTINE puff
194 END MODULE puff_mod
max
#define max(a, b)
Definition macros.h:21
constant_mod
Definition constant_mod.F:355
puff_mod
Definition puff.F:29
puff_mod::puff
subroutine puff(iflag, nel, pmin, off, eint, mu, mu2, espe, dvol, df, vnew, pnew, dpdm, dpde, eos_struct)
Definition puff.F:45