OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
tillotson.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!|| tillotson ../common_source/eos/tillotson.F
25!||--- called by ------------------------------------------------------
26!|| eosmain ../common_source/eos/eosmain.F
27!||====================================================================
28 SUBROUTINE tillotson(
29 1 IFLAG,NEL ,PM ,OFF ,EINT ,MU , MU2,
30 2 ESPE ,DVOL ,DF ,VNEW ,MAT ,PSH ,
31 3 PNEW ,DPDM ,DPDE ,VAREOS, NVAREOS)
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
241 END
alpha
#define alpha
Definition eval.h:35
max
#define max(a, b)
Definition macros.h:21
tillotson
subroutine tillotson(iflag, nel, pm, off, eint, mu, mu2, espe, dvol, df, vnew, mat, psh, pnew, dpdm, dpde, vareos, nvareos)
Definition tillotson.F:32