OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps105.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!|| sigeps105 ../engine/source/materials/mat/mat105/sigeps105.F
25!||--- called by ------------------------------------------------------
26!|| mulaw ../engine/source/materials/mat_share/mulaw.F90
27!||--- calls -----------------------------------------------------
28!|| finter ../engine/source/tools/curve/finter.F
29!||--- uses -----------------------------------------------------
30!|| elbufdef_mod ../common_source/modules/mat_elem/elbufdef_mod.F90
31!|| i22bufbric_mod ../common_source/modules/interfaces/cut-cell-search_mod.F
32!|| i22tri_mod ../common_source/modules/interfaces/cut-cell-search_mod.F
33!||====================================================================
34 SUBROUTINE sigeps105 (
35 1 NEL ,NUPARAM,NUVAR ,NFUNC ,IFUNC ,TBURN ,
36 2 NPF ,TF ,TIME ,TIMESTEP,UPARAM ,BFRAC ,
37 3 RHO0 ,RHO ,VOLUME ,EINT ,SIGY ,DELTAX ,
38 4 EPSPXX ,EPSPYY ,EPSPZZ ,EPSPXY ,EPSPYZ ,EPSPZX ,
39 5 DEPSXX ,DEPSYY ,DEPSZZ ,DEPSXY ,DEPSYZ ,DEPSZX ,
40 6 EPSXX ,EPSYY ,EPSZZ ,EPSXY ,EPSYZ ,EPSZX ,
41 7 SIGOXX ,SIGOYY ,SIGOZZ ,SIGOXY ,SIGOYZ ,SIGOZX ,
42 8 SIGNXX ,SIGNYY ,SIGNZZ ,SIGNXY ,SIGNYZ ,SIGNZX ,
43 9 SIGVXX ,SIGVYY ,SIGVZZ ,SIGVXY ,SIGVYZ ,SIGVZX ,
44 A SOUNDSP,VISCMAX,UVAR ,OFF ,NFT ,V ,
45 B W ,X ,IX ,NIX ,JTHE ,
46 C GEO ,PID ,ILAY ,NG ,ELBUF_TAB,PM ,
47 D IPARG ,BUFVOIS ,IPM ,BUFMAT ,STIFN ,
48 E VD2 ,VDX ,VDY ,VDZ ,MAT ,VOLN ,
49 F QNEW ,DDVOL ,QOLD ,PSH)
50C-----------------------------------------------
51C D e s c r i p t i o n
52C-----------------------------------------------
53C This subroutine computes Powder Burn EOS and is based on (Atwood,Friis,Moxnes) publication.
54C ELEMENT BUFFER
55C %BFRAC : is ignition fraction (depends on flamme position regarding the cell)
56C UVAR(6) : is burn fraction (depends on F(t) function and also on gas pressure and grains shape)
57C EOS
58C GAS : Exponantial EOS
59C POWDER : Linear EOS
60C MATERIAL BUFFER
61C UPARAM(1:NUPARAM) : material parameters
62C UVAR (1:NUVAR) : cell buffer specific for this material law : Pgas, etc...
63C
64C An Alternative implementation is to use /EOS option : see ./common_source/eos/powder_burn.F
65C
66C-----------------------------------------------
67C M o d u l e s
68C-----------------------------------------------
69 USE elbufdef_mod
70 USE i22bufbric_mod
71 USE i22tri_mod
72C-----------------------------------------------
73C I m p l i c i t T y p e s
74C-----------------------------------------------
75#include "implicit_f.inc"
76#include "comlock.inc"
77C-----------------------------------------------
78C C o m m o n B l o c k s
79C-----------------------------------------------
80#include "param_c.inc"
81#include "com01_c.inc"
82#include "com08_c.inc"
83#include "mvsiz_p.inc"
84C-----------------------------------------------
85C I N P U T A r g u m e n t s
86C-----------------------------------------------
87 INTEGER NEL, NUPARAM, NUVAR,NFT,NIX,JTHE,
88 . IX(NIX,*), PID(*), ILAY, NG,PM(NPROPM,*),IPARG(*),MAT(*),
89 . IPM(NPROPMI,*)
90 my_real
91 . TIME ,TIMESTEP ,UPARAM(NUPARAM),
92 . SIGY(NEL) ,VK(NEL) ,
93 . RHO(NEL) ,RHO0(NEL) ,VOLUME(NEL),
94 . EINT(NEL) ,BUFVOIS(*) ,QNEW(NEL) ,
95 . DDVOL(NEL) ,QOLD(NEL) ,
96 . EPSPXX(NEL) ,EPSPYY(NEL),EPSPZZ(NEL),
97 . EPSPXY(NEL) ,EPSPYZ(NEL),EPSPZX(NEL),
98 . DEPSXX(NEL) ,DEPSYY(NEL),DEPSZZ(NEL),
99 . DEPSXY(NEL) ,DEPSYZ(NEL),DEPSZX(NEL),
100 . EPSXX(NEL) ,EPSYY(NEL) ,EPSZZ(NEL) ,
101 . EPSXY(NEL) ,EPSYZ(NEL) ,EPSZX(NEL) ,
102 . sigoxx(nel) ,sigoyy(nel),sigozz(nel),
103 . sigoxy(nel) ,sigoyz(nel),sigozx(nel),
104 . v(3,*) ,w(3,*) ,x(3,*) ,
105 . geo(npropg,*),bufmat(*) ,
106 . stifn(mvsiz) ,vd2(*) ,vdx(*) ,
107 . vdy(*) ,vdz(*) ,ssp ,
108 . voln(*)
109 my_real,intent(inout) :: psh(nel)
110 TYPE (ELBUF_STRUCT_), TARGET, DIMENSION(NGROUP) :: ELBUF_TAB
111C-----------------------------------------------
112C O U T P U T A r g u m e n t s
113C-----------------------------------------------
114 my_real
115 . SIGNXX(NEL),SIGNYY(NEL),SIGNZZ(NEL),
116 . SIGNXY(NEL),SIGNYZ(NEL),SIGNZX(NEL),
117 . SIGVXX(NEL),SIGVYY(NEL),SIGVZZ(NEL),
118 . SIGVXY(NEL),SIGVYZ(NEL),SIGVZX(NEL),
119 . SOUNDSP(NEL),VISCMAX(NEL)
120C-----------------------------------------------
121C I N P U T O U T P U T A r g u m e n t s
122C-----------------------------------------------
123 my_real UVAR(NEL,NUVAR), OFF(NEL), TBURN(NEL), BFRAC(NEL), DELTAX(NEL)
124C-----------------------------------------------
125C VARIABLES FOR FUNCTION INTERPOLATION
126C-----------------------------------------------
127 INTEGER NPF(*), NFUNC, IFUNC(NFUNC)
128 my_real TF(*)
129C-----------------------------------------------
130C E x t e r n a l F u n c t i o n s
131C-----------------------------------------------
132 my_real,EXTERNAL :: finter
133 ! EXTERNAL FINTER
134 ! Y = FINTER(IFUNC(J),X,NPF,TF,DYDX)
135 ! Y : y = f(x)
136 ! X : x
137 ! DYDX : f'(x) = dy/dx
138 ! IFUNC(J): FUNCTION INDEX
139 ! J : FIRST(J=1), SECOND(J=2) .. FUNCTION USED FOR THIS LAW
140 ! NPF,TF : FUNCTION PARAMETER
141C-----------------------------------------------
142C L o c a l V a r i a b l e s
143C-----------------------------------------------
144 my_real :: p,qa,qb,qal,qbl, pold,vold,einc,pnew
145 my_real :: bulk, d, eg, gr, c, alpha,beta, c1, c2,fscale_g,fscale_b,fscale_p,fscale_rho
146 my_real :: p0,xl, dd, volo, espe, e0, total_bfrac, delta_bf, brate
147 my_real :: pg,pg_,ps, mu_g, mu_p, mu, rho_g,rho_g_, rho_s, rhog, rho_si
148 INTEGER :: funcb, funcg
149 INTEGER :: IBID, IBFRAC, I , K
150 my_real :: MASS, DF, TB, DIFF, G_RHOC2, P_RHOC2, RHOC2, SSP_G, SSP_S, MASS_G, MASS_S, VS
151 my_real :: VOL_G, VOL_S, fac,compac,INVrhog
152 my_real :: m0,rho_s0,F,Y_old,Y,FUNC_DOT,FUNC,V0,TMP
153
154 TYPE(G_BUFEL_) ,POINTER :: GBUF
155 TYPE(L_BUFEL_) ,POINTER :: LBUF
156 TYPE(BUF_LAY_) ,POINTER :: BUFLY
157C-----------------------------------------------
158C S o u r c e L i n e s
159C-----------------------------------------------
160 gbuf => elbuf_tab(ng)%GBUF
161 lbuf => elbuf_tab(ng)%BUFLY(ilay)%LBUF(1,1,1)
162 bufly => elbuf_tab(ng)%BUFLY(ilay)
163C-----------------------------------------------
164C S t a n d a r d J W L E O S
165C-----------------------------------------------
166 bulk = uparam(01)
167 p0 = uparam(02)
168 psh(:) = uparam(03)
169 e0 = uparam(04)
170 d = uparam(05)
171 eg = uparam(06)
172 gr = uparam(07)
173 c = uparam(08)
174 alpha = uparam(09)
175 c1 = uparam(10)
176 c2 = uparam(11)
177 fscale_g = uparam(12)
178 fscale_b = uparam(13)
179 fscale_p = uparam(14)
180 fscale_rho = uparam(15)
181 funcb = ifunc(1)
182 funcg = ifunc(2)
183
184 compac=0.93
185
186 IF(dt1 == zero)THEN
187 DO i=1,nel
188 eint(i) = eg*rho0(i)*volume(i)
189 uvar(i,1) = p0-psh(i) !PS
190 uvar(i,2) = p0-psh(i) !PG
191 uvar(i,3) = rho0(i)/compac !RHO_S
192 uvar(i,4) = zero !RHO_G
193 uvar(i,5) = p0-psh(i) !POLD
194 uvar(i,6) = zero !F(t_old)
195 uvar(i,7) = rho0(i)*volume(i) !Mass0
196 ENDDO
197 signxx(1:nel) = - (p0-psh(i))
198 signyy(1:nel) = - (p0-psh(i))
199 signzz(1:nel) = - (p0-psh(i))
200 signxy(1:nel) = zero
201 signyz(1:nel) = zero
202 signzx(1:nel) = zero
203 sigvxx(1:nel) = zero
204 sigvyy(1:nel) = zero
205 sigvzz(1:nel) = zero
206 sigvxy(1:nel) = zero
207 sigvyz(1:nel) = zero
208 sigvzx(1:nel) = zero
209 DO i=1,nel
210 soundsp(i)= sqrt(bulk/rho0(i))
211 ENDDO
212 ENDIF
213
214 DO i=1,nel
215 !--------------------------------!
216 ! INIT. !
217 !--------------------------------!
218 mass = rho(i)*volume(i)
219 xl = deltax(i)
220 espe = eint(i)/mass
221 ps = uvar(i,1)
222 pg = uvar(i,2)
223 rho_s = uvar(i,3)
224 rho_g = uvar(i,4)
225 pold = uvar(i,5)
226
227 !--------------------------------!
228 ! IGNITION FRACTION !
229 !--------------------------------!
230 !--> valid for const flamme velocity
231 IF(bfrac(i) < one) THEN
232 tb = - tburn(i)
233 bfrac(i) = zero
234 IF(time > tb) bfrac(i) = c1*(time-tb)*two_third/xl
235 IF(bfrac(i) < em04) THEN
236 bfrac(i) = zero
237 ELSEIF(bfrac(i) > one) THEN
238 bfrac(i) = one
239 ENDIF
240 ENDIF
241
242 !------------------------------------------------!
243 ! BURN FRACTION !
244 ! BFRAC : ignition fraction (flamme evolution) !
245 ! UVAR(6): burn fraction (grain burn evolution) !
246 !------------------------------------------------!
247 IF(bfrac(i) <= zero)THEN
248 total_bfrac = zero
249 uvar(i,6) = zero
250 ELSEIF(uvar(i,6) /= one)THEN ! F<1
251 pg_=pg
252 !PG_=PG*(ONE-MASS_S/VOLN(I)/RHO_S) !effective gaz pressure Pg*(1-RHO_S_/rhop)
253 brate = fscale_b*finter(funcb,(pg_)/fscale_p,npf,tf,diff)
254 delta_bf = gr*exp(c*log((one-alpha*uvar(i,6))))*brate*timestep !F_dot * DT
255 uvar(i,6) = uvar(i,6) + delta_bf !F
256 uvar(i,6) = max(min(one,uvar(i,6)),zero)
257 total_bfrac = bfrac(i)*uvar(i,6)
258 ELSE
259 ! avoid LOG(0)
260 total_bfrac = bfrac(i)*uvar(i,6)
261 ENDIF
262
263 !--------------------------------!
264 ! EOS SOLVING !
265 !--------------------------------!
266 mass_s=(one-uvar(i,6))*uvar(i,7) !(1-F)*M0
267 mass_g=uvar(i,7)-mass_s !Mass_Gas=F(t)*Mass
268
269 rho_s=(one-uvar(i,6))*rho(i) !MASS_S/VOLN(I)
270 rho_g=zero
271
272 !in case of compressible solid (linear EOS)
273 rho_si=rho0(i)/compac*(uvar(i,2)/bulk+one)
274 !in case of constant solid density
275 !rho_si=rho0(I)/compac
276
277 IF(mass_g > zero)rho_g=mass_g/(voln(i)-mass_s/ max(em10,rho_si))
278 pg = p0+rho_g*eg*exp(rho_g/d)
279 ps = p0+bulk*(rho_si/rho0(i)-one)
280 ssp_s = sqrt(bulk/rho0(i))
281 IF(rho_g > zero)THEN
282 ssp_g = rho0(i)/d*pg + exp(rho_g/d)* pg / (rho_g/rho0(i))**2
283 ssp_g = sqrt(one/rho_g * ssp_g)
284 ELSE
285 ssp_g = zero
286 ENDIF
287
288 !--------------------------------!
289 ! Global Pressure & SoundSpeed !
290 !--------------------------------!
291 ssp = total_bfrac*ssp_g + (one-total_bfrac)*ssp_s
292 pnew = total_bfrac*pg + (one-total_bfrac)*ps
293 pnew = (pnew-psh(i))*off(i)
294 !--------------------------------!
295 ! Returning values !
296 !--------------------------------!
297 signxx(i) = -pnew
298 signyy(i) = -pnew
299 signzz(i) = -pnew
300 signxy(i) = zero
301 signyz(i) = zero
302 signzx(i) = zero
303 sigvxx(i) = zero
304 sigvyy(i) = zero
305 sigvzz(i) = zero
306 sigvxy(i) = zero
307 sigvyz(i) = zero
308 sigvzx(i) = zero
309 soundsp(i) = ssp
310
311 uvar(i,1)=ps
312 uvar(i,2)=pg
313 uvar(i,3)=rho_s
314 uvar(i,4)=rho_g
315 uvar(i,5)=pnew
316 !UVAR(I,6) -> F(t)
317 !UVAR(I,7) -> Mass0
318
319 ! --- verification output
320 if(time == zero)then
321 open(unit=1,file='F.txt' )
322 open(unit=2,file='Pg.txt' )
323 open(unit=3,file='rhoS.txt' )
324 open(unit=4,file='Ifrac.txt' )
325 else
326 open(unit=1,file='F.txt' , position="append")
327 open(unit=2,file='Pg.txt' , position="append")
328 open(unit=3,file='rhoS.txt', position="append")
329 open(unit=4,file='Ifrac.txt',position="append")
330 endif
331
332 WRITE(1,*)time,uvar(i,6)
333 WRITE(2,*)time,pg
334 WRITE(3,*)time,rho_s
335 WRITE(4,*)bfrac(i)
336
337 CLOSE(1)
338 CLOSE(2)
339 CLOSE(3)
340
341 enddo! next I
342
343C-----------------------------------------------
344 RETURN
345C-----------------------------------------------
346 END SUBROUTINE sigeps105
347C-----------------------------------------------
alpha
#define alpha
Definition eval.h:35
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
i22bufbric_mod
Definition cut-cell-search_mod.F:95
i22tri_mod
Definition cut-cell-search_mod.F:212
sigeps105
subroutine sigeps105(nel, nuparam, nuvar, nfunc, ifunc, tburn, npf, tf, time, timestep, uparam, bfrac, rho0, rho, volume, eint, sigy, deltax, epspxx, epspyy, epspzz, epspxy, epspyz, epspzx, depsxx, depsyy, depszz, depsxy, depsyz, depszx, epsxx, epsyy, epszz, epsxy, epsyz, epszx, sigoxx, sigoyy, sigozz, sigoxy, sigoyz, sigozx, signxx, signyy, signzz, signxy, signyz, signzx, sigvxx, sigvyy, sigvzz, sigvxy, sigvyz, sigvzx, soundsp, viscmax, uvar, off, nft, v, w, x, ix, nix, jthe, geo, pid, ilay, ng, elbuf_tab, pm, iparg, bufvois, ipm, bufmat, stifn, vd2, vdx, vdy, vdz, mat, voln, qnew, ddvol, qold, psh)
Definition sigeps105.F:50