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

Functions/Subroutines

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

Function/Subroutine Documentation

◆ tillotson()

subroutine tillotson_mod::tillotson ( integer, intent(in) iflag,
integer, intent(in) nel,
real(kind=wp) pmin,
real(kind=wp), dimension(nel) off,
real(kind=wp), dimension(nel) eint,
real(kind=wp), dimension(nel) mu,
real(kind=wp), dimension(nel) mu2,
real(kind=wp), dimension(nel) espe,
real(kind=wp), dimension(nel) dvol,
real(kind=wp), dimension(nel) df,
real(kind=wp), dimension(nel) vnew,
real(kind=wp), dimension(nel), intent(inout) psh,
real(kind=wp), dimension(nel) pnew,
real(kind=wp), dimension(nel) dpdm,
real(kind=wp), dimension(nel) dpde,
real(kind=wp), dimension(nel,nvareos), intent(inout) vareos,
integer, intent(in) nvareos,
type(eos_param_), intent(in) eos_struct )

Definition at line 42 of file tillotson.F.

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
alpha
#define alpha
Definition eval.h:35
max
#define max(a, b)
Definition macros.h:21
constant_mod
Definition constant_mod.F:355