OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
tillotson.F File Reference
#include "implicit_f.inc"
#include "param_c.inc"
#include "com04_c.inc"
#include "com06_c.inc"
#include "com08_c.inc"
#include "vect01_c.inc"
#include "scr06_c.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine tillotson (iflag, nel, pm, off, eint, mu, mu2, espe, dvol, df, vnew, mat, psh, pnew, dpdm, dpde, vareos, nvareos)

Function/Subroutine Documentation

◆ tillotson()

subroutine tillotson ( integer, intent(in) iflag,
integer, intent(in) nel,
pm,
off,
eint,
mu,
mu2,
espe,
dvol,
df,
vnew,
integer, dimension(nel), intent(in) mat,
dimension(nel), intent(inout) psh,
pnew,
dpdm,
dpde,
dimension(nel,nvareos), intent(inout) vareos,
integer, intent(in) nvareos )

Definition at line 28 of file tillotson.F.

32C-----------------------------------------------
33C D e s c r i p t i o n
34C-----------------------------------------------
35!----------------------------------------------------------------------------
36!! \details STAGGERED SCHEME IS EXECUTED IN TWO PASSES IN EOSMAIN : IFLG=0 THEN IFLG=1
37!! \details COLLOCATED SCHEME IS DOING A SINGLE PASS : IFLG=2
38!! \details
39!! \details STAGGERED SCHEME
40!! \details EOSMAIN / IFLG = 0 : DERIVATIVE CALCULATION FOR SOUND SPEED ESTIMATION c[n+1] REQUIRED FOR PSEUDO-VISCOSITY (DPDE:partial derivative, DPDM:total derivative)
41!! \details MQVISCB : PSEUDO-VISCOSITY Q[n+1]
42!! \details MEINT : INTERNAL ENERGY INTEGRATION FOR E[n+1] : FIRST PART USING P[n], Q[n], and Q[n+1] CONTRIBUTIONS
43!! \details EOSMAIN / IFLG = 1 : UPDATE P[n+1], T[N+1]
44!! \details INTERNAL ENERGY INTEGRATION FOR E[n+1] : LAST PART USING P[n+1] CONTRIBUTION
45!! \details (second order integration dE = -P.dV where P = 0.5(P[n+1] + P[n]) )
46!! \details COLLOCATED SCHEME
47!! \details EOSMAIN / IFLG = 2 : SINGLE PASS FOR P[n+1] AND DERIVATIVES
48!----------------------------------------------------------------------------
49C
50C Solver for Tillotson Equation of State.
51C This equation of state depends of material state. It has several region in (P,v) diagram :
52C REGION=1 / (I) : compression
53C REGION=2 / (II) : cold expansion
54C REGION=3 / (III): transition (currently empty but can be implemented)
55C REGION=4 / (IV) : hot expansion
56C
57C in following source code regions are identified with following criteria :
58C !IF(MU(I) => ZERO)THEN
59C ! REGION = 1
60C !ELSEIF(MU(I) < ZERO)THEN
61C ! REGION = 2
62C ! IF( V>VSUBL_(I) .OR. (V<VSUBL .AND. ESPE(I)>=ESUBL_(I)) )THEN
63C ! REGION=4
64C ! ENDIF
65C !ENDIF
66C-----------------------------------------------
67C I m p l i c i t T y p e s
68C-----------------------------------------------
69#include "implicit_f.inc"
70C-----------------------------------------------
71C C o m m o n B l o c k s
72C-----------------------------------------------
73#include "param_c.inc"
74#include "com04_c.inc"
75#include "com06_c.inc"
76#include "com08_c.inc"
77#include "vect01_c.inc"
78#include "scr06_c.inc"
79C-----------------------------------------------
80C D u m m y A r g u m e n t s
81C-----------------------------------------------
82 INTEGER,INTENT(IN) :: MAT(NEL), IFLAG, NEL, NVAREOS
83 my_real pm(npropm,nummat),
84 . off(nel) , eint(nel), mu(nel) ,
85 . mu2(nel) , espe(nel), dvol(nel), df(nel) ,
86 . vnew(nel), pnew(nel), dpdm(nel), dpde(nel)
87 my_real,INTENT(INOUT) :: vareos(nel,nvareos)
88 my_real,INTENT(INOUT) :: psh(nel)
89C-----------------------------------------------
90C L o c a l V a r i a b l e s
91C-----------------------------------------------
92 INTEGER I, MX
93 my_real aa, bb, dvv, eta, enew, omega, xx, expa, expb,
94 . pp, facc1, facc2, facpb,
95 . c1(nel),c2(nel),ptia(nel),ptib(nel),ezero(nel),
96 . alpha(nel),beta(nel),esubl(nel),vsubl(nel),
97 . pc(nel), region
98C--------------------------------------------------------------------
99 IF(iflag == 0) THEN
100 DO i=1,nel
101 mx = mat(i)
102 c1(i) = pm(32,mx)
103 c2(i) = pm(33,mx)
104 ptia(i) = pm(34,mx)
105 ptib(i) = pm(35,mx)
106 pc(i) = pm(37,mx)
107 ezero(i)= pm(36,mx)
108 psh(i) = pm(88,mx)
109 esubl(i)= pm(160,mx)
110 vsubl(i)= pm(161,mx)
111 alpha(i)= pm(162,mx)
112 beta(i) = pm(163,mx)
113 ENDDO
114 DO i=1,nel
115 facc1=one
116 facc2=one
117 facpb=one
118 region=1
119 IF(mu(i)<zero) THEN
120 region=2
121 facc2=zero
122 IF(df(i)> vsubl(i) .OR. (df(i)<=vsubl(i).AND.espe(i)>=esubl(i))) THEN
123 region=4
124 xx=mu(i)/(one+mu(i))
125 expa=exp(-alpha(i)*xx*xx)
126 expb=exp(beta(i)*xx)
127 facc1=expa*expb
128 facpb=expa
129 ENDIF
130 ENDIF
131 eta=one+mu(i)
132 omega= one+espe(i)/(ezero(i)*eta**2)
133 aa=facc1*c1(i)*mu(i)+facc2*c2(i)*mu2(i)
134 bb=ptia(i)+facpb*ptib(i)/omega
135 pp=max(aa+bb*eta*espe(i),pc(i))*off(i)
136 dpdm(i)=facc1*c1(i)+two*facc2*c2(i)*mu(i)+bb*eta*pp*df(i)*df(i)
137 . +espe(i)*( bb+(two*espe(i)/eta-pp*df(i)*df(i))
138 . *ptib(i)*facpb/(ezero(i)*eta*omega**2) )
139 dpde(i)=bb*eta
140 vareos(i,1)=region
141 pnew(i) = max(pp,pc(i))*off(i)! P(mu[n+1],E[n])
142 ENDDO
143
144 ELSEIF(iflag == 1) THEN
145 DO i=1,nel
146 mx = mat(i)
147 c1(i) = pm(32,mx)
148 c2(i) = pm(33,mx)
149 ptia(i) = pm(34,mx)
150 ptib(i) = pm(35,mx)
151 pc(i) = pm(37,mx)
152 ezero(i)= pm(36,mx)
153 psh(i) = pm(88,mx)
154 esubl(i)= pm(160,mx)
155 vsubl(i)= pm(161,mx)
156 alpha(i)= pm(162,mx)
157 beta(i) = pm(163,mx)
158 ENDDO
159
160 DO i=1,nel
161 dvv=half*dvol(i)*df(i) / max(em15,vnew(i))
162 eta=one+mu(i)
163 omega= one+espe(i)/(ezero(i)*eta**2)
164 facc1=one
165 facc2=one
166 facpb=one
167 region=1
168 IF(mu(i)<zero) THEN
169 region=2
170 facc2=zero
171 IF(df(i)>vsubl(i).OR.(df(i)<=vsubl(i).AND.espe(i)>=esubl(i))) THEN
172 region=4
173 xx=mu(i)/(one+mu(i))
174 expa=exp(-alpha(i)*xx*xx)
175 expb=exp(beta(i)*xx)
176 facc1=expa*expb
177 facpb=expa
178 ENDIF
179 ENDIF
180 aa=facc1*c1(i)*mu(i)+facc2*c2(i)*mu2(i)
181 bb=(ptia(i)+facpb*ptib(i)/omega)*eta
182 pnew(i)=(aa +bb*espe(i))/(one+ bb*dvv)
183 enew=espe(i) - pnew(i)*dvv
184 !one iteration
185 omega= one+enew/(ezero(i)*eta**2)
186 bb=(ptia(i)+facpb*ptib(i)/omega)*eta
187 pnew(i)=(aa +bb*espe(i))/(one+ bb*dvv)
188 pnew(i)= max(pnew(i),pc(i))*off(i)! P(mu[n+1],E[n+1])
189 eint(i)= eint(i) - half*dvol(i)*pnew(i)
190 dpde(i) = bb
191 vareos(i,1)=region
192 ENDDO
193
194 ELSEIF(iflag == 2) THEN
195 DO i=1, nel
196 mx = mat(i)
197 c1(i) = pm(32,mx)
198 c2(i) = pm(33,mx)
199 ptia(i) = pm(34,mx)
200 ptib(i) = pm(35,mx)
201 pc(i) = pm(37,mx)
202 ezero(i)= pm(36,mx)
203 esubl(i)= pm(160,mx)
204 vsubl(i)= pm(161,mx)
205 alpha(i)= pm(162,mx)
206 beta(i) = pm(163,mx)
207 ENDDO
208 DO i=1, nel
209 region=1
210 IF (vnew(i) > zero) THEN
211 facc1=one
212 facc2=one
213 facpb=one
214 region=1
215 IF(mu(i)<zero) THEN
216 region=2
217 facc2=zero
218 IF(df(i) > vsubl(i) .OR. (df(i) <= vsubl(i) .AND. espe(i) >= esubl(i))) THEN
219 region=4
220 xx = mu(i)/(one+mu(i))
221 expa= exp(-alpha(i)*xx*xx)
222 expb= exp(beta(i)*xx)
223 facc1=expa*expb
224 facpb=expa
225 ENDIF
226 ENDIF
227 eta=one+mu(i)
228 omega= one+espe(i)/(ezero(i)*eta**2)
229 aa=facc1*c1(i)*mu(i)+facc2*c2(i)*mu2(i)
230 bb=ptia(i)+facpb*ptib(i)/omega
231 pp=max(aa+bb*eta*espe(i),pc(i))*off(i)
232 dpdm(i)=facc1*c1(i)+two*facc2*c2(i)*mu(i)+bb*eta*pp*df(i)*df(i)
233 . +espe(i)*( bb+(two*espe(i)/eta-pp*df(i)*df(i))
234 . *ptib(i)*facpb/(ezero(i)*eta*omega**2) )
235 dpde(i)=bb*eta
236 ENDIF
237 vareos(i,1)=region
238 ENDDO
239 ENDIF
240 RETURN
my_real
#define my_real
Definition cppsort.cpp:32
alpha
#define alpha
Definition eval.h:35
max
#define max(a, b)
Definition macros.h:21