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_mod ../common_source/eos/tillotson.F
25!||--- called by ------------------------------------------------------
26!|| eosmain ../common_source/eos/eosmain.F
27!|| eosmain51 ../engine/source/materials/mat/mat051/eosmain51.F90
28!||====================================================================
29 MODULE tillotson_mod
30 IMPLICIT NONE
31 CONTAINS
32!||====================================================================
33!|| tillotson ../common_source/eos/tillotson.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!|| precision_mod ../common_source/modules/precision_mod.F90
41!||====================================================================
42 SUBROUTINE tillotson(
43 1 IFLAG,NEL ,PMIN ,OFF ,EINT ,MU , MU2,
44 2 ESPE ,DVOL ,DF ,VNEW ,PSH ,
45 3 PNEW ,DPDM ,DPDE ,VAREOS, NVAREOS, EOS_STRUCT)
46C-----------------------------------------------
47C M o d u l e s
48C-----------------------------------------------
49 USE constant_mod , ONLY : zero, em15, half, one, two
50 USE eos_param_mod , ONLY : eos_param_
51 use precision_mod , only : wp
52 IMPLICIT NONE
53C-----------------------------------------------
54C D e s c r i p t i o n
55C-----------------------------------------------
56!----------------------------------------------------------------------------
57!! \details STAGGERED SCHEME IS EXECUTED IN TWO PASSES IN EOSMAIN : IFLG=0 THEN IFLG=1
58!! \details COLLOCATED SCHEME IS DOING A SINGLE PASS : IFLG=2
59!! \details
60!! \details STAGGERED SCHEME
61!! \details EOSMAIN / IFLG = 0 : DERIVATIVE CALCULATION FOR SOUND SPEED ESTIMATION c[n+1] REQUIRED FOR PSEUDO-VISCOSITY (DPDE:partial derivative, DPDM:total derivative)
62!! \details MQVISCB : PSEUDO-VISCOSITY Q[n+1]
63!! \details MEINT : INTERNAL ENERGY INTEGRATION FOR E[n+1] : FIRST PART USING P[n], Q[n], and Q[n+1] CONTRIBUTIONS
64!! \details EOSMAIN / IFLG = 1 : UPDATE P[n+1], T[N+1]
65!! \details INTERNAL ENERGY INTEGRATION FOR E[n+1] : LAST PART USING P[n+1] CONTRIBUTION
66!! \details (second order integration dE = -P.dV where P = 0.5(P[n+1] + P[n]) )
67!! \details COLLOCATED SCHEME
68!! \details EOSMAIN / IFLG = 2 : SINGLE PASS FOR P[n+1] AND DERIVATIVES
69!----------------------------------------------------------------------------
70C
71C Solver for Tillotson Equation of State.
72C This equation of state depends of material state. It has several region in (P,v) diagram :
73C REGION=1 / (I) : compression
74C REGION=2 / (II) : cold expansion
75C REGION=3 / (III): transition (currently empty but can be implemented)
76C REGION=4 / (IV) : hot expansion
77C
78C in following source code regions are identified with following criteria :
79C !IF(MU(I) => ZERO)THEN
80C ! REGION = 1
81C !ELSEIF(MU(I) < ZERO)THEN
82C ! REGION = 2
83C ! IF( V>VSUBL_(I) .OR. (V<VSUBL .AND. ESPE(I)>=ESUBL_(I)) )THEN
84C ! REGION=4
85C ! ENDIF
86C !ENDIF
87C-----------------------------------------------
88C D u m m y A r g u m e n t s
89C-----------------------------------------------
90 INTEGER,INTENT(IN) :: IFLAG, NEL, NVAREOS
91 real(kind=WP) :: PMIN,
92 . OFF(NEL) , EINT(NEL), MU(NEL) ,
93 . mu2(nel) , espe(nel), dvol(nel), df(nel),
94 . vnew(nel), pnew(nel), dpdm(nel), dpde(nel)
95 real(kind=wp),INTENT(INOUT) :: vareos(nel,nvareos)
96 real(kind=wp),INTENT(INOUT) :: psh(nel)
97 TYPE(eos_param_),INTENT(IN) :: eos_struct
98C-----------------------------------------------
99C L o c a l V a r i a b l e s
100C-----------------------------------------------
101 INTEGER :: i
102 real(kind=WP) :: aa, bb, dvv, eta, enew, omega, xx, expa, expb,
103 . pp, facc1, facc2, facpb,
104 . c1,c2,ptia,ptib,ezero,
105 . alpha,beta,esubl,vsubl,
106 . region
107C--------------------------------------------------------------------
108 c1 = eos_struct%UPARAM(1)
109 c2 = eos_struct%UPARAM(2)
110 ptia = eos_struct%UPARAM(3)
111 ptib = eos_struct%UPARAM(4)
112 ezero = eos_struct%UPARAM(5)
113 esubl = eos_struct%UPARAM(6)
114 vsubl = eos_struct%UPARAM(7)
115 alpha = eos_struct%UPARAM(8)
116 beta = eos_struct%UPARAM(9)
117 psh(:) = eos_struct%PSH
118
119 IF(iflag == 0) THEN
120 DO i=1,nel
121 facc1=one
122 facc2=one
123 facpb=one
124 region=1
125 IF(mu(i) < zero) THEN
126 region=2
127 facc2=zero
128 IF(df(i) > vsubl .OR. (df(i) <= vsubl .AND. espe(i) >= esubl)) THEN
129 region=4
130 xx=mu(i)/(one+mu(i))
131 expa=exp(-alpha*xx*xx)
132 expb=exp(beta*xx)
133 facc1=expa*expb
134 facpb=expa
135 ENDIF
136 ENDIF
137 eta=one+mu(i)
138 omega= one+espe(i)/(ezero*eta**2)
139 aa=facc1*c1*mu(i)+facc2*c2*mu2(i)
140 bb=ptia+facpb*ptib/omega
141 pp=max(aa+bb*eta*espe(i),pmin)*off(i)
142 dpdm(i)=facc1*c1+two*facc2*c2*mu(i)+bb*eta*pp*df(i)*df(i)
143 . +espe(i)*( bb+(two*espe(i)/eta-pp*df(i)*df(i))
144 . *ptib*facpb/(ezero*eta*omega**2) )
145 dpde(i)=bb*eta
146 vareos(i,1)=region
147 pnew(i) = max(pp,pmin)*off(i)! P(mu[n+1],E[n])
148 ENDDO
149
150 ELSEIF(iflag == 1) THEN
151
152 DO i=1,nel
153 dvv=half*dvol(i)*df(i) / max(em15,vnew(i))
154 eta=one+mu(i)
155 omega= one+espe(i)/(ezero*eta**2)
156 facc1=one
157 facc2=one
158 facpb=one
159 region=1
160 IF(mu(i) < zero) THEN
161 region=2
162 facc2=zero
163 IF(df(i) > vsubl .OR. (df(i) <= vsubl .AND. espe(i) >= esubl)) THEN
164 region=4
165 xx=mu(i)/(one+mu(i))
166 expa=exp(-alpha*xx*xx)
167 expb=exp(beta*xx)
168 facc1=expa*expb
169 facpb=expa
170 ENDIF
171 ENDIF
172 aa=facc1*c1*mu(i)+facc2*c2*mu2(i)
173 bb=(ptia+facpb*ptib/omega)*eta
174 pnew(i)=(aa +bb*espe(i))/(one+ bb*dvv)
175 enew=espe(i) - pnew(i)*dvv
176 !one iteration
177 omega= one+enew/(ezero*eta**2)
178 bb=(ptia+facpb*ptib/omega)*eta
179 pnew(i)=(aa +bb*espe(i))/(one+ bb*dvv)
180 pnew(i)= max(pnew(i),pmin)*off(i)! P(mu[n+1],E[n+1])
181 eint(i)= eint(i) - half*dvol(i)*pnew(i)
182 dpde(i) = bb
183 vareos(i,1)=region
184 ENDDO
185
186 ELSEIF(iflag == 2) THEN
187
188 DO i=1, nel
189 region=1
190 IF (vnew(i) > zero) THEN
191 facc1=one
192 facc2=one
193 facpb=one
194 region=1
195 IF(mu(i)<zero) THEN
196 region=2
197 facc2=zero
198 IF(df(i) > vsubl .OR. (df(i) <= vsubl .AND. espe(i) >= esubl)) THEN
199 region=4
200 xx = mu(i)/(one+mu(i))
201 expa= exp(-alpha*xx*xx)
202 expb= exp(beta*xx)
203 facc1=expa*expb
204 facpb=expa
205 ENDIF
206 ENDIF
207 eta=one+mu(i)
208 omega= one+espe(i)/(ezero*eta**2)
209 aa=facc1*c1*mu(i)+facc2*c2*mu2(i)
210 bb=ptia+facpb*ptib/omega
211 pp=max(aa+bb*eta*espe(i),pmin)*off(i)
212 dpdm(i)=facc1*c1+two*facc2*c2*mu(i)+bb*eta*pp*df(i)*df(i)
213 . +espe(i)*( bb+(two*espe(i)/eta-pp*df(i)*df(i))
214 . *ptib*facpb/(ezero*eta*omega**2) )
215 dpde(i)=bb*eta
216 pnew(i) = pp
217 ENDIF
218 vareos(i,1)=region
219 ENDDO
220 ENDIF
221 RETURN
222 END SUBROUTINE tillotson
223 END MODULE tillotson_mod
alpha
#define alpha
Definition eval.h:35
max
#define max(a, b)
Definition macros.h:21
constant_mod
Definition constant_mod.F:355
tillotson_mod
Definition tillotson.F:29
tillotson_mod::tillotson
subroutine tillotson(iflag, nel, pmin, off, eint, mu, mu2, espe, dvol, df, vnew, psh, pnew, dpdm, dpde, vareos, nvareos, eos_struct)
Definition tillotson.F:46