42
43
44
49 USE matparam_def_mod
51 USE reader_old_mod , ONLY : irec
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66#include "implicit_f.inc"
67
68
69
70#include "scr17_c.inc"
71#include "units_c.inc"
72#include "param_c.inc"
73
74
75
76 TYPE (UNIT_TYPE_),INTENT(IN) ::UNITAB
77 INTEGER, INTENT(IN) :: MAT_ID
78 INTEGER, INTENT(INOUT) :: ISRATE
79 INTEGER, DIMENSION(NPROPMI) ,INTENT(INOUT) :: IPM
80 CHARACTER(LEN=NCHARTITLE) ,INTENT(IN) :: TITR
81 my_real,
DIMENSION(NPROPM) ,
INTENT(INOUT) :: pm
82 TYPE(SUBMODEL_DATA), DIMENSION(NSUBMOD),INTENT(IN) :: LSUBMODEL
83 TYPE(MLAW_TAG_), INTENT(INOUT) :: MTAG
84 TYPE(EOS_TAG_) , TARGET, DIMENSION(0:MAXEOS) ,INTENT(INOUT) :: EOS_TAG
85 TYPE(MATPARAM_STRUCT_),INTENT(INOUT) :: MATPARAM
86
87
88
89 LOGICAL :: IS_AVAILABLE,IS_ENCRYPTED
90 INTEGER :: COUNT,IFORM,IEOS
91 my_real :: rho0,rhor,a0,a1,a2,amx,c0,c1,c2,c3
92 my_real :: pmin,pext,xmumx,bunl,e0,g,pstar,psh,delta,det,young,anu
93 CHARACTER*64 :: chain,chain2
94
95
96
97 pstar = -infinity
98 c0 = zero
99 c1 = zero
100 c2 = zero
101 c3 = zero
102 pext = zero
103 psh = zero
104 pmin = zero
105 bunl = zero
106 xmumx = zero
107 israte = 0
108 is_encrypted = .false.
109 is_available = .false.
110
112
113 CALL hm_get_intv (
'Line_count' ,count ,is_available,lsubmodel)
114 IF(count==0)count=3
115 count=count+2
116
117 CALL hm_get_floatv(
'MAT_RHO' ,rho0 ,is_available, lsubmodel, unitab)
118 CALL hm_get_floatv(
'Refer_Rho' ,rhor ,is_available, lsubmodel, unitab)
119 CALL hm_get_floatv(
'MAT_E' ,young ,is_available, lsubmodel, unitab)
120 CALL hm_get_floatv(
'MAT_NU' ,anu ,is_available, lsubmodel, unitab)
121 CALL hm_get_floatv(
'MAT_A0' ,a0 ,is_available, lsubmodel, unitab)
122 CALL hm_get_floatv(
'MAT_A1' ,a1 ,is_available, lsubmodel, unitab)
123 CALL hm_get_floatv(
'MAT_A2' ,a2 ,is_available, lsubmodel, unitab)
124 CALL hm_get_floatv(
'MAT_AMAX' ,amx ,is_available, lsubmodel, unitab)
125 CALL hm_get_floatv(
'EOS_COM_C0' ,c0 ,is_available, lsubmodel, unitab)
126 CALL hm_get_floatv(
'EOS_COM_C1' ,c1 ,is_available, lsubmodel, unitab)
127 CALL hm_get_floatv(
'EOS_COM_C2' ,c2 ,is_available, lsubmodel, unitab)
128 CALL hm_get_floatv(
'EOS_COM_C3' ,c3 ,is_available, lsubmodel, unitab)
129 CALL hm_get_floatv(
'EOS_COM_B' ,bunl ,is_available, lsubmodel, unitab)
130 CALL hm_get_floatv(
'EOS_COM_Mue_max' ,xmumx ,is_available, lsubmodel, unitab)
131
132 CALL hm_get_floatv(
'MAT_PC' ,pmin ,is_available, lsubmodel, unitab)
133 CALL hm_get_floatv(
'PEXT' ,pext ,is_available, lsubmodel, unitab)
134
135
136 IF (rhor == zero) rhor=rho0
137 pm(1) = rhor
138 pm(89)= rho0
139
140 IF(count == 5) THEN
141 iform= 1
142 ieos = 13
143 eos_tag(ieos)%L_MU = 1
144 matparam%IEOS = ieos
145
146 matparam%EOS%NUPARAM = 7
147 matparam%EOS%NIPARAM = 1
148 matparam%EOS%NFUNC = 0
149 matparam%EOS%NTABLE = 0
150 CALL matparam%EOS%CONSTRUCT()
151 matparam%EOS%uparam(1) = xmumx
152 matparam%EOS%uparam(2) = zero
153 matparam%EOS%uparam(3) = bunl
154 matparam%EOS%uparam(4) = c0
155 matparam%EOS%uparam(5) = c1
156 matparam%EOS%uparam(6) = c2
157 matparam%EOS%uparam(7) = c3
158 matparam%EOS%iparam(1) = 1
159 matparam%EOS%psh = psh
160 ELSE
161 iform=2
162 ieos = 18
163 irec = irec+1
164 ENDIF
165
166 matparam%IEOS = ieos
167 ipm(4) = ieos
168
169 e0 = zero
170
171 IF(c1 == zero)c1=third*young/(one-two*anu)
172
173 IF(young <= zero)THEN
174 chain='YOUNG MODULUS MUST BE DEFINED '
175 chain2=''
176 CALL ancmsg(msgid=829, msgtype=msgerror, anmode=aninfo, i1=10, i2=mat_id, c1=
'ERROR', c2=titr, c3=chain, c4=chain2)
177 ENDIF
178
179 IF(anu <= zero)THEN
180 chain='POISSON RATIO MUST BE DEFINED '
181 chain2=''
182 CALL ancmsg(msgid=829,msgtype=msgerror,anmode=aninfo,i1=10,i2=mat_id,c1=
'ERROR',c2=titr,c3=chain,c4=chain2
183 ENDIF
184
185 IF(a1 < zero .AND. a2 == zero)THEN
186 chain ='INVERTED YIELD SURFACE BECAUSE A1 IS NEGATIVE. '
187 chain2='CHECK A0,A1,A2 YIELD PARAMETERS '
188 CALL ancmsg(msgid=829,msgtype=msgwarning,anmode=aninfo,i1=10,i2=mat_id,c1=
'WARNING',c2=titr,c3=chain,c4=chain2)
189 ENDIF
190
191 IF(a2 < zero)THEN
192 chain ='UNEXPECTED YIELD SURFACE : A2 IS NEGATIVE '
193 chain2='CHECK A0,A1,A2 YIELD PARAMETERS '
194 CALL ancmsg(msgid=829,msgtype=msgwarning,anmode=aninfo,i1=10,i2=mat_id,c1=
'WARNING',c2=titr,c3=chain,c4=chain2)
195 ENDIF
196 IF(iform == 1 .AND. c1 <= zero) THEN
197 chain='TENSILE BULK MODULUS C1 IS LOWER OR EQUAL TO 0. '
198 chain2=''
199 CALL ancmsg(msgid=829,msgtype=msgerror,anmode=aninfo,i1=10,i2=mat_id,c1=
'ERROR',c2=titr,c3=chain,c4=chain2)
200 END IF
201
202 IF(a2==zero.AND.a1/=zero)THEN
203 pstar=-a0/a1
204 ELSEIF(a2 /= zero)THEN
205 delta = a1*a1-four*a0*a2
206
207 IF(delta >= zero)THEN
208 delta=sqrt(delta)
209 pstar = (-a1+delta)/two/a2
210
211 ELSE
212 pstar = -a1/two/a2
213 chain ='YIELD SURFACE J2(P)=A0+A1.P+A2.P^2 HAS NO INTERSECTION WITH '
214 chain2='PRESSURE AXIS. ASSUMING SURFACE CLOSURE AT P= '
215 WRITE(chain2(46:59),fmt=('(E13.6)')) pstar
216 CALL ancmsg( msgid=829,msgtype=msgwarning,anmode=aninfo,i1=10,i2=mat_id,c1=
'WARNING',c2=titr,c3=chain,c4=chain2)
217 ENDIF
218 ELSE
219
220 ENDIF
221
222 IF(amx == zero) amx = ep20
223 IF(pmin == zero) pmin =-infinity
224
225
226
227
228 IF(iform == 1)THEN
229 IF(xmumx == zero.AND.bunl /= zero)THEN
230 chain= 'MISSING MUMAX VALUE IS AUTOMATICALLY ESTIMATED FROM BUNL VALUE. '
231 chain2=''
232 CALL ancmsg(msgid=829,msgtype=msgwarning,anmode=aninfo,i1=10,i2=mat_id,c1=
'WARNING',c2=titr,c3=chain,c4=chain2)
233 IF(c3 == zero)THEN
234 IF(c2 == zero)THEN
235 xmumx=ep20
236 ELSE
237 xmumx=(bunl-c1)/(two*c2)
238 ENDIF
239 ELSE
240 det=sqrt(c2**2 + three*c3*(bunl-c1))
241 xmumx=(det-c2)/(three*c3)
242 ENDIF
243 ENDIF
244
245 IF(xmumx /= zero.AND.bunl == zero)THEN
246 chain= 'MISSING BUNL VALUE IS AUTOMATICALLY ESTIMATED FROM MUMAX '
247 chain2=''
248 CALL ancmsg(msgid=829,msgtype=msgwarning,anmode=aninfo,i1=10,i2=mat_id,c1=
'WARNING',c2=titr,c3=chain,c4=chain2)
249 IF(c3 == zero)THEN
250 IF(c2 == zero)THEN
251 bunl = c1
252 ELSE
253 bunl = c1 + two*c2*xmumx
254 ENDIF
255 ELSE
256 bunl = c1 + two*c2*xmumx + three*c3*c3*xmumx**two
257 ENDIF
258 ENDIF
259 IF(bunl == zero) bunl = c1
260 ENDIF
261
262 g=young/(two*(one + anu))
263
264
265 pm(20)=young
266 pm(21)=anu
267 pm(22)=g
268 pm(23)=e0
269 pm(31)=c0
270 pm(32)=c1
271 pm(33)=c2
272 pm(34)=c3
273 pm(49)=c0
274 pm(35)=zero
275 pm(36)=zero
276 pm(37)=pmin
277 pm(38)=a0
278 pm(39)=a1
279 pm(40)=a2
280 pm(41)=amx
281 pm(43)=pext
282 pm(44)=pstar
283 pm(45)=bunl
284 pm(46)=xmumx
285 pm(47)=zero
286 pm(48)=iform
287 pm(79)=three100
288 pm(88)=psh
289 pm(104)=c0
290 pm(105) = two*g/(c1+four_over_3*g)
291
292
293 ipm(252)= 2
294
295
296
297 mtag%G_PLA = 1
298 mtag%G_EPSQ = 1
299 mtag%G_MU = 1
300
301 mtag%L_PLA = 1
302 mtag%L_EPSQ = 1
303 mtag%L_MU = 1
304
305
307
308 ! EOS/Thermo keyword for pressure treatment in elements
309 CALL INIT_MAT_KEYWORD(MATPARAM,"hydro_eos")
310
311 CALL INIT_MAT_KEYWORD(MATPARAM,"elasto_plastic")
312
313 ! Properties compatibility
314 CALL INIT_MAT_KEYWORD(MATPARAM,"solid_isotropic")
315 CALL INIT_MAT_KEYWORD(MATPARAM,"sph")
316
317 !--- OUTPUT IN STARTER LISTING FILE ---!
318 WRITE(IOUT,1100) TRIM(TITR),MAT_ID,10
319 WRITE(IOUT,1000)
320 IF(IS_ENCRYPTED)THEN
321 WRITE(IOUT,'(5X,A,//)')'CONFIDENTIAL DATA'
322 ELSE
323 WRITE(IOUT,1200)RHO0,RHOR
324 WRITE(IOUT,1300)YOUNG,ANU,G
325 WRITE(IOUT,1400)A0,A1,A2,AMX
326 IF(IFORM == 1)THEN
327 WRITE(IOUT,1500)C0,C1,C2,C3,BUNL,XMUMX,PMIN,PSTAR
328 ELSE
329 WRITE(IOUT,1501)BUNL,XMUMX,PMIN,PSTAR
330 ENDIF
331 ENDIF
332
333 RETURN
334 1000 FORMAT(
335 & 5X,40H DRUCKER-PRAGER (LAW10) ,/,
336 & 5X,40H ---------------------- ,//)
337 1100 FORMAT(/
338 & 5X,A,/,
339 & 5X,'MATERIAL NUMBER . . . . . . . . . . . .=',I10/,
340 & 5X,'MATERIAL LAW. . . . . . . . . . . . . .=',I10/)
341 1200 FORMAT(
342 & 5X,'INITIAL DENSITY . . . . . . . . . . . .=',1PG20.13/,
343 & 5X,'REFERENCE DENSITY . . . . . . . . . . .=',1PG20.13/)
344 1300 FORMAT(
345 & 5X,40HYOUNG'S MODULUS . . . . . . . . . . . .=,E12.4/,
346 & 5X,40HPOISSON'S RATIO . . . . . . . . . . . .=,E12.4/,
347 & 5X,40HSHEAR MODULUS . . . . . . . . . . . . .=,E12.4//)
348 1400 FORMAT(
349 & 5X,40HYIELD COEFFICIENT A0. . . . . . . . . .=,E12.4/,
350 & 5X,40HYIELD COEFFICIENT A1. . . . . . . . . .=,E12.4/,
351 & 5X,40HYIELD COEFFICIENT A2. . . . . . . . . .=,E12.4/,
352 & 5X,40HA-MAX . . . . . . . . . . . . . . . . .=,E12.4//)
353 1500 FORMAT(
354 & 5X,40HC0. . . . . . . . . . . . . . . . . . .=,E12.4/,
355 & 5X,40HC1. . . . . . . . . . . . . . . . . . .=,E12.4/,
356 & 5X,40HC2. . . . . . . . . . . . . . . . . . .=,E12.4/,
357 & 5X,40HC3. . . . . . . . . . . . . . . . . . .=,E12.4/,
358 & 5X,40HUNLOADING BULK. . . . . . . . . . . . .=,E12.4/,
359 & 5X,40HMAX VOLUMIC COMPRESSION . . . . . . . .=,E12.4/,
360 & 5X,40HFRACTURE PRESSURE . . . . . . . . . . .=,E12.4/,
361 & 5X,40HYIELD SURFACE PRESSURE ROOT . . . . . .=,E12.4//)
362 1501 FORMAT(
363 & 5X,40HUNLOADING BULK. . . . . . . . . . . . .=,E12.4/,
364 & 5X,40HMAX VOLUMIC COMPRESSION . . . . . . . .=,E12.4/,
365 & 5X,40HFRACTURE PRESSURE . . . . . . . . . . .=,E12.4/,
366 & 5X,40HYIELD SURFACE PRESSURE ROOT . . . . . .=,E12.4//)
367
368 RETURN
subroutine hm_get_floatv(name, rval, is_available, lsubmodel, unitab)
subroutine hm_get_intv(name, ival, is_available, lsubmodel)
subroutine hm_option_is_encrypted(is_encrypted)
subroutine init_mat_keyword(matparam, keyword)
integer, parameter nchartitle
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)