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

Functions/Subroutines

subroutine puff (iflag, nel, pmin, off, eint, mu, mu2, espe, dvol, df, vnew, pnew, dpdm, dpde, eos_struct)

Function/Subroutine Documentation

◆ puff()

subroutine puff_mod::puff ( 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) pnew,
dimension(nel), intent(inout) dpdm,
dimension(nel), intent(inout) dpde,
type(eos_param_), intent(in) eos_struct )

Definition at line 41 of file puff.F.

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
my_real
#define my_real
Definition cppsort.cpp:32
max
#define max(a, b)
Definition macros.h:21
constant_mod
Definition constant_mod.F:355