OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
lecm51__check_initial_state.F File Reference
#include "implicit_f.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine lecm51__check_initial_state (av, r0, c0, c1, c2, c3, c4, c5, e0, pm, rho0, rhor, iexp, pext, iflg, a1, a2, pla, id, titr, ssp1, ssp2, ssp3, ssp4, lc1, lc2, lc3, lc4, p0a)

Function/Subroutine Documentation

◆ lecm51__check_initial_state()

subroutine lecm51__check_initial_state ( av,
r0,
c0,
c1,
c2,
c3,
c4,
c5,
e0,
pm,
rho0,
rhor,
integer iexp,
pext,
integer iflg,
a1,
a2,
integer, dimension(4) pla,
integer id,
character(len=nchartitle) titr,
ssp1,
ssp2,
ssp3,
ssp4,
lc1,
lc2,
lc3,
lc4,
p0a )

Definition at line 32 of file lecm51__check_initial_state.F.

40 USE message_mod
41 USE names_and_titles_mod , ONLY : nchartitle
42C-----------------------------------------------
43C D e s c r i p t i o n
44C-----------------------------------------------
45C Check initial state of the different material
46C
47C 1 INIVOL FRAC in [0 1]
48C 2 AV1+AV2+AV3+AV4 = 1 (tol 0.1%)
49C 3 Negative or null density
50C 5 Unbalanced initial Pressure
51C 4 Global Density not consistent with material
52C densities
53C-----------------------------------------------
54C I m p l i c i t T y p e s
55C-----------------------------------------------
56#include "implicit_f.inc"
57C-----------------------------------------------
58C D u m m y A r g u m e n t s
59C-----------------------------------------------
60 my_real
61 . av(4),r0(4), c0(4), c1(4),c2(4), c3(4),
62 . a1(4),a2(4),
63 . c4(4),c5(4), e0(4), pm(4), pext, rho0, rhor,
64 . ssp(4),lc(4),
65 . ssp1,ssp2,ssp3,ssp4,
66 . lc1,lc2,lc3,lc4,
67 . p0a(4)
68 INTEGER IEXP, IFLG,PLA(4)
69 INTEGER ID
70 CHARACTER(LEN=NCHARTITLE)::TITR
71C-----------------------------------------------
72C L o c a l V a r i a b l e s
73C-----------------------------------------------
74 INTEGER i,j, TEST1,TEST2,TEST3,TEST4,TEST5,
75 . TEST6, IMAX,IPLA(4)
76
77 my_real sum, p0b, err
78
79 character*8 chain
80 character*47 chain3
81 character*64 chain1
82C===============================================|
83 ipla(1)=pla(1)
84 ipla(2)=pla(2)
85 ipla(3)=pla(3)
86 ipla(4)=pla(4)
87
88 imax=iexp+3
89 !IMAX number of phases to be checked 3 or 4
90
91 ssp(1:4) = (/ssp1,ssp2,ssp3,ssp4/)
92 lc(1:4) = (/ lc1, lc2, lc3, lc4/)
93
94 !0: warning or error not detected
95 !1: test failed, message in listing file.
96 test1=0 !AV in [0,1]
97 test2=0 !Sum(AV(i))=1
98 test3=0 !Rho not defined
99 test4=0 !Unbalanced Pressure t=0
100 test5=0 !Global density = mean(rho_i)
101 test6=0 !C1,C4,C5 >=0, C1>0 phase 4
102
103 IF(iflg==6)THEN
104 sum=0
105 DO i=1,imax
106 sum=sum+av(i)*r0(i)
107 END DO
108 IF(sum/=zero)THEN
109 rho0=sum
110 rhor=rho0
111 ELSE
112 rho0=one !will be computed in engine, not need for sound speed. This formulation has SSP=EM20
113 rhor=one
114 ENDIF
115 RETURN
116 ENDIF
117 !TEST1 (quel que soit iflg)
118 ! Initial Volume fraction checked : must be between 0 and 1
119 ! Otherwise : error
120 DO i=1,imax
121 IF ((av(i)<zero).OR.(av(i)>one)) THEN
122 chain='SUBMAT-0'
123 write(chain(8:8),'(i1)')i
124 chain1='INITIAL VOLUMETRIC FRACTION MUST BE BETWEEN 0 AND 1 FOR '//chain(1:8)
125 CALL ancmsg(msgid=99,msgtype=msgerror,anmode=aninfo,i1=id,c1=titr,c2=chain1)
126 test1=1
127 !exit
128 END IF
129 END DO
130
131
132 !TEST2 (quel que soit iflg)
133 ! Initial Volume fraction : Sum must be 1
134 ! Otherwise : warning
135 ! if TEST1/=0 it is useless to do this test)
136 sum=zero
137 IF (test1==0) THEN
138 DO i=1,imax
139 sum=sum+av(i) ! AV(I) in [0,1], SUM=0 <=> AV(I)=0 , i=1..4
140 END DO
141 IF(sum==zero)THEN
142 av(1)=1
143 CALL ancmsg(msgid=98,msgtype=msgwarning,anmode=aninfo_blind_1,i1=id,c1=titr,
144 . c2='INITIAL VOLUMETRIC FRACTIONS ARE NULL , SUBMAT-1 SET TO 100%')
145 ELSEIF (abs(1-sum)>em03) THEN
146 CALL ancmsg(msgid=99,msgtype=msgerror,anmode=aninfo_blind_1,i1=id,c1=titr,
147 . c2='SUM OF INITIAL VOLUMETRIC FRACTION IS NOT EQUAL TO 1')
148 test2=1
149 END IF
150 END IF
151
152
153 !TEST3(whatever is iflg)
154 !Negative Density or Missing density material
155 ! No negative density
156 ! If AV>0 THEN R0 must be >0
157 DO i=1,imax
158 IF ( ((r0(i)<zero).OR.((av(i)>zero).AND.(r0(i)<=zero))) ) THEN
159 chain='SUBMAT-0'
160 write(chain(8:8),'(i1.1)')i
161 !print *,"chain(1:8) =", chain(1:8)
162 chain3='NULL OR NEGATIVE INITIAL DENSITIES FOR '//chain
163 CALL ancmsg(msgid=99,msgtype=msgerror,anmode=aninfo,i1=id,c1=titr,c2=chain3)
164 test3=1
165 !exit
166 END IF
167 END DO
168
169
170 !TEST4 (whatever is iflg)
171 !Unbalanced Pressure, material are checked two by two (Cnp couples)
172 !IF (TEST1+TEST2+TEST3==0) THEN !If no error
173 sum=zero
174 IF(imax==4)p0a(4) = c0(4)
175 DO i=1,imax-1
176 if (r0(i)==zero)cycle !if empty submaterial
177 p0a(i)=c0(i)+c4(i)*e0(i)
178 DO j=i,imax
179 IF (r0(j)==zero)cycle !if empty submaterial
180 !comparing Pi(t=0) with Pj(t=0)
181 p0b=c0(j)+c4(j)*e0(j)
182 err=abs(p0a(i)-p0b)
183 IF (err<=em20) err=zero !avoid epsilon issue if P0a ou P0b is zero
184 IF (err/=zero) THEN ! JWL under pressure may detonate.
185 IF (err/max(abs(p0a(i)),abs(p0b))>em06) THEN
186 IF (j==4) THEN
187 test4=1
188 CALL ancmsg(msgid=98,msgtype=msgwarning,anmode=aninfo_blind_1,
189 . i1=id,c1=titr,c2='PRESSURES ARE UNBALANCED WITH JWL MATERIAL')
190 exit
191 ELSE
192 test4=1
193 CALL ancmsg(msgid=98,msgtype=msgwarning,anmode=aninfo_blind_1,
194 . i1=id,c1=titr,c2='INITIAL PRESSURES ARE UNBALANCED')
195 EXIT
196 END IF
197 END IF
198 END IF
199 END DO !next J
200 IF(test4==1)EXIT
201 END DO !next I
202 !END IF
203
204
205 !TEST5 (whatever is iflg)
206 !Check the consistency between global density and materials densities.
207 !IF(TEST1+TEST2+TEST3==0)THEN !If no error
208 sum=0
209 DO i=1,imax
210 sum=sum+av(i)*r0(i)
211 END DO
212 test5=1
213 rho0=sum
214 rhor=rho0
215 !END IF
216
217
218 !TEST6 (whatever is iflg)
219 ! check if c1>=0, c4>=0, c5>=0
220 ! Otherwise : error
221 DO i=1,3
222 IF ((c1(i)<zero)) THEN
223 chain='SUBMAT-0'
224 write(chain(8:8),'(i1.1)')i
225 chain1='POLYNOMIAL COEFFICIENT C1 MUST BE POSITIVE OR NULL FOR '//chain
226 CALL ancmsg(msgid=99,msgtype=msgerror,anmode=aninfo,i1=id,c1=titr,c2=chain1)
227 test6=1
228 !exit
229 END IF
230
231 IF ((c2(i)<zero)) THEN
232 chain='submat-0'
233 write(chain(8:8),'(i1.1)')I
234 chain1='polynomial coefficient c2 is negative for '//chain
235 CALL ANCMSG(MSGID=98,MSGTYPE=MSGWARNING,ANMODE=ANINFO_BLIND_1,I1=ID,C1=TITR,C2=chain1)
236 TEST6=1
237 !exit
238 END IF
239
240 IF ((C3(I)<ZERO)) THEN
241 chain='submat-0'
242 write(chain(8:8),'(i1.1)')I
243 chain1='polynomial coefficient c3 is negative for '//chain
244 CALL ANCMSG(MSGID=98,MSGTYPE=MSGWARNING,ANMODE=ANINFO_BLIND_1,I1=ID,C1=TITR,C2=chain1)
245 TEST6=1
246 !exit
247 END IF
248
249 IF ((C4(I)<ZERO)) THEN
250 chain='submat-0'
251 write(chain(8:8),'(i1.1)')I
252 chain1='polynomial coefficient c4 must be positive or null for '//chain
253 CALL ANCMSG(MSGID=99,MSGTYPE=MSGERROR,ANMODE=ANINFO_BLIND_1,I1=ID,C1=TITR,C2=chain1)
254 TEST6=1
255 !exit
256 END IF
257
258 IF ((C5(I)<ZERO)) THEN
259 chain='submat-0'
260 write(chain(8:8),'(i1.1)')I
261 chain1='polynomial coefficient c5 must be positive or null for '//chain
262 CALL ANCMSG(MSGID=99,MSGTYPE=MSGERROR,ANMODE=ANINFO_BLIND_1,I1=ID,C1=TITR,C2=chain1)
263 TEST6=1
264 !exit
265 END IF
266 END DO !next I
267
268
269 !check pressure cut-off pressure for perfect gas
270 IF(IFLG<=1)THEN
271 DO I=1,3
272 IF(IPLA(I)>0)CYCLE
273 P0a(I) = PEXT+C0(I)+C4(I)*E0(I)
274.AND..AND. IF(((PEXT+C0(I))==C1(I))(C2(I)==ZERO)
275.AND..AND. . (C3(I)==ZERO)(C4(I)==C5(I))
276.AND. . (C4(I)>ZERO)(P0a(I)>ZERO))THEN
277 !we have a perfect gas
278.AND. IF(PM(I)/=ZERO PEXT/=ZERO)THEN
279 IF(abs(PM(I)+PEXT)>EM06)THEN
280 !PM(I)=-PEXT
281 chain='submat-0'
282 write(chain(8:8),'(i1.1)')I
283 chain1=
284 . chain//' is a perfect gas:minimum pressure should be -pext '
285 CALL ANCMSG(MSGID=98,MSGTYPE=MSGWARNING,ANMODE=ANINFO_BLIND_1,I1=ID,C1=TITR,C2=chain1)
286 END IF
287 ENDIF
288 END IF
289 END DO !next I
290 END IF !(IFLG<=1)
291
292 !check Drucker Prager Yield Surface
293 IF(IFLG<=1)THEN
294 DO I=1,3
295 chain='submat-0'
296 write(chain(8:8),'(i1.1)')I
297 !
298.AND. IF(A1(I) < ZERO A2(I) == ZERO)THEN
299 chain1 = chain//': inverted yield surface. check a1 sign. '
300 CALL ANCMSG(MSGID=829,MSGTYPE=MSGWARNING,ANMODE=ANINFO,I1=51,I2=ID,C1='warning',C2=TITR,C3=chain)
301 ENDIF
302 IF(A2(I) < ZERO)THEN
303 chain1 = chain//': untipycal yield surface. check a2 sign. '
304 CALL ANCMSG(MSGID=829,MSGTYPE=MSGWARNING,ANMODE=ANINFO,I1=51,I2=ID,C1='warning',C2=TITR,C3=chain)
305 ENDIF
306 !
307 END DO !next I
308 END IF !(IFLG<=1)
309
310 RETURN
my_real
#define my_real
Definition cppsort.cpp:32
max
#define max(a, b)
Definition macros.h:21
id
initmumps id
Definition multiplerhs_example.m:4
for
for(i8=*sizetab-1;i8 >=0;i8--)
Definition mumps_common.c:91
message_mod
Definition message_mod.F:1249
names_and_titles_mod
Definition names_and_titles_mod.F:997
names_and_titles_mod::nchartitle
integer, parameter nchartitle
Definition names_and_titles_mod.F:1003
ancmsg
subroutine ancmsg(msgid, msgtype, anmode, i1, i2, i3, i4, i5, i6, i7, i8, i9, i10, i11, i12, i13, i14, i15, i16, i17, i18, i19, i20, r1, r2, r3, r4, r5, r6, r7, r8, r9, c1, c2, c3, c4, c5, c6, c7, c8, c9, prmode)
Definition message.F:889