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 :: I,IFUNC,COUNT,ID,IFORM,IEOS
91 my_real :: rho0,rhor,e,nu,a0,a1,a2,amx,c0,c1,c2,c3,fac_y
92 my_real :: pmin,pext,xmumx,bunl,e0,g,pstar,psh,delta,det,young
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 ELSE
146 iform=2
147 ieos = 18
148 irec = irec+1
149 ENDIF
150
151 matparam%IEOS = ieos
152 ipm(4) = ieos
153
154 e0 = zero
155
156 IF(c1 == zero)c1=third*young/(one-two*anu)
157
158 IF(young <= zero)THEN
159 chain='YOUNG MODULUS MUST BE DEFINED '
160 chain2=''
161 CALL ancmsg(msgid=829, msgtype=msgerror, anmode=aninfo, i1=10, i2=mat_id, c1=
'ERROR', c2=titr, c3=chain, c4=chain2)
162 ENDIF
163
164 IF(anu <= zero)THEN
165 chain='POISSON RATIO MUST BE DEFINED '
166 chain2=''
167 CALL ancmsg(msgid=829,msgtype=msgerror,anmode=aninfo,i1=10,i2=mat_id,c1=
'ERROR',c2=titr,c3=chain,c4=chain2)
168 ENDIF
169
170 IF(a1 < zero .AND. a2 == zero)THEN
171 chain ='INVERTED YIELD SURFACE BECAUSE A1 IS NEGATIVE. '
172 chain2='CHECK A0,A1,A2 YIELD PARAMETERS '
173 CALL ancmsg(msgid=829,msgtype=msgwarning,anmode=aninfo,i1=10,i2=mat_id,c1=
'WARNING',c2=titr,c3=chain,c4=chain2)
174 ENDIF
175
176 IF(a2 < zero)THEN
177 chain ='UNEXPECTED YIELD SURFACE : A2 IS NEGATIVE '
178 chain2='CHECK A0,A1,A2 YIELD PARAMETERS '
179 CALL ancmsg(msgid=829,msgtype=msgwarning,anmode=aninfo,i1=10,i2=mat_id,c1=
'WARNING',c2=titr,c3=chain,c4=chain2)
180 ENDIF
181 IF(iform == 1 .AND. c1 <= zero) THEN
182 chain='TENSILE BULK MODULUS C1 IS LOWER OR EQUAL TO 0. '
183 chain2=''
184 CALL ancmsg(msgid=829,msgtype=msgerror,anmode=aninfo,i1=10,i2=mat_id,c1=
'ERROR',c2=titr,c3=chain,c4=chain2)
185 END IF
186
187 IF(a2==zero.AND.a1/=zero)THEN
188 pstar=-a0/a1
189 ELSEIF(a2 /= zero)THEN
190 delta = a1*a1-four*a0*a2
191
192 IF(delta >= zero)THEN
193 delta=sqrt(delta)
194 pstar = (-a1+delta)/two/a2
195
196 ELSE
197 pstar = -a1/two/a2
198 chain ='YIELD SURFACE J2(P)=A0+A1.P+A2.P^2 HAS NO INTERSECTION WITH '
199 chain2='PRESSURE AXIS. ASSUMING SURFACE CLOSURE AT P= '
200 WRITE(chain2(46:59),fmt=('(E13.6)')) pstar
201 CALL ancmsg( msgid=829,msgtype=msgwarning,anmode=aninfo,i1=10,i2=mat_id,c1=
'WARNING',c2=titr,c3=chain,c4=chain2)
202 ENDIF
203 ELSE
204
205 ENDIF
206
207 IF(amx == zero) amx = ep20
208 IF(pmin == zero) pmin =-infinity
209
210
211
212
213 IF(iform == 1)THEN
214 IF(xmumx == zero.AND.bunl /= zero)THEN
215 chain= 'MISSING MUMAX VALUE IS AUTOMATICALLY ESTIMATED FROM BUNL VALUE. '
216 chain2=''
217 CALL ancmsg(msgid=829,msgtype=msgwarning,anmode=aninfo,i1=10,i2=mat_id,c1=
'WARNING',c2=titr,c3=chain,c4=chain2)
218 IF(c3 == zero)THEN
219 IF(c2 == zero)THEN
220 xmumx=ep20
221 ELSE
222 xmumx=(bunl-c1)/(two*c2)
223 ENDIF
224 ELSE
225 det=sqrt(c2**2 + three*c3*(bunl-c1))
226 xmumx=(det-c2)/(three*c3)
227 ENDIF
228 ENDIF
229
230 IF(xmumx /= zero.AND.bunl == zero)THEN
231 chain= 'MISSING BUNL VALUE IS AUTOMATICALLY ESTIMATED FROM MUMAX '
232 chain2=''
233 CALL ancmsg(msgid=829,msgtype=msgwarning,anmode=aninfo,i1=
'WARNING',c2=titr,c3=chain,c4=chain2)
234 IF(c3 == zero)THEN
235 IF(c2 == zero)THEN
236 bunl = c1
237 ELSE
238 bunl = c1 + two*c2*xmumx
239 ENDIF
240 ELSE
241 bunl = c1 + two*c2*xmumx + three*c3*c3*xmumx**two
242 ENDIF
243 ENDIF
244 IF(bunl == zero) bunl = c1
245 ENDIF
246
247 g=young/(two*(one + anu))
248
249
250 pm(20)=young
251 pm(21)=anu
252 pm(22)=g
253 pm(23)=e0
254 pm(31)=c0
255 pm(32)=c1
256 pm(33)=c2
257 pm(34)=c3
258 pm(49)=c0
259 pm(35)=zero
260 pm(36)=zero
261 pm(37)=pmin
262 pm(38)=a0
263 pm(39)=a1
264 pm(40)=a2
265 pm(41)=amx
266 pm(43)=pext
267 pm(44)=pstar
268 pm(45)=bunl
269 pm(46)=xmumx
270 pm(47)=zero
271 pm(48)=iform
272 pm(79)=three100
273 pm(88)=psh
274 pm(104)=c0
275 pm(105) = two*g/(c1+four_over_3*g)
276
277
278 ipm(252)= 2
279
280
281
282 mtag%G_PLA = 1
283 mtag%G_EPSQ = 1
284 mtag%G_MU = 1
285
286 mtag%L_PLA = 1
287 mtag%L_EPSQ = 1
288 mtag%L_MU = 1
289
290 IF(count == 5) THEN
291
292 matparam%eos%nuparam = 3
293 matparam%eos%niparam = 1
294 matparam%eos%nfunc = 0
295 matparam%eos%ntable = 0
296 call matparam%eos%construct()
297 matparam%eos%uparam(1) = xmumx
298 matparam%eos%uparam(2) = zero
299 matparam%eos%uparam(3) = bunl
300 matparam%eos%iparam(1) = iform
301 ENDIF
302
303
305
306 ! EOS/Thermo keyword for pressure treatment in elements
307 CALL INIT_MAT_KEYWORD(MATPARAM,"hydro_eos")
308
309 CALL INIT_MAT_KEYWORD(MATPARAM,"elasto_plastic")
310
311 ! Properties compatibility
312 CALL INIT_MAT_KEYWORD(MATPARAM,"solid_isotropic")
313 CALL INIT_MAT_KEYWORD(MATPARAM,"sph")
314
315 !--- OUTPUT IN STARTER LISTING FILE ---!
316 WRITE(IOUT,1100) TRIM(TITR),MAT_ID,10
317 WRITE(IOUT,1000)
318 IF(IS_ENCRYPTED)THEN
319 WRITE(IOUT,'(5X,A,//)')'CONFIDENTIAL DATA'
320 ELSE
321 WRITE(IOUT,1200)RHO0,RHOR
322 WRITE(IOUT,1300)YOUNG,ANU,G
323 WRITE(IOUT,1400)A0,A1,A2,AMX
324 IF(IFORM == 1)THEN
325 WRITE(IOUT,1500)C0,C1,C2,C3,BUNL,XMUMX,PMIN,PSTAR
326 ELSE
327 WRITE(IOUT,1501)BUNL,XMUMX,PMIN,PSTAR
328 ENDIF
329 ENDIF
330
331 RETURN
332 1000 FORMAT(
333 & 5X,40H DRUCKER-PRAGER (LAW10) ,/,
334 & 5X,40H ---------------------- ,//)
335 1100 FORMAT(/
336 & 5X,A,/,
337 & 5X,'MATERIAL NUMBER . . . . . . . . . . . .=',I10/,
338 & 5X,'MATERIAL LAW. . . . . . . . . . . . . .=',I10/)
339 1200 FORMAT(
340 & 5X,'INITIAL DENSITY . . . . . . . . . . . .=',1PG20.13/,
341 & 5X,'REFERENCE DENSITY . . . . . . . . . . .=',1PG20.13/)
342 1300 FORMAT(
343 & 5X,40HYOUNG'S MODULUS . . . . . . . . . . . .=,E12.4/,
344 & 5X,40HPOISSON'S RATIO . . . . . . . . . . . .=,E12.4/,
345 & 5X,40HSHEAR MODULUS . . . . . . . . . . . . .=,E12.4//)
346 1400 FORMAT(
347 & 5X,40HYIELD COEFFICIENT A0. . . . . . . . . .=,E12.4/,
348 & 5X,40HYIELD COEFFICIENT A1. . . . . . . . . .=,E12.4/,
349 & 5X,40HYIELD COEFFICIENT A2. . . . . . . . . .=,E12.4/,
350 & 5X,40HA-MAX . . . . . . . . . . . . . . . . .=,E12.4//)
351 1500 FORMAT(
352 & 5X,40HC0. . . . . . . . . . . . . . . . . . .=,E12.4/,
353 & 5X,40HC1. . . . . . . . . . . . . . . . . . .=,E12.4/,
354 & 5X,40HC2. . . . . . . . . . . . . . . . . . .=,E12.4/,
355 & 5X,40HC3. . . . . . . . . . . . . . . . . . .=,E12.4/,
356 & 5X,40HUNLOADING BULK. . . . . . . . . . . . .=,E12.4/,
357 & 5X,40HMAX VOLUMIC COMPRESSION . . . . . . . .=,E12.4/,
358 & 5X,40HFRACTURE PRESSURE . . . . . . . . . . .=,E12.4/,
359 & 5X,40HYIELD SURFACE PRESSURE ROOT . . . . . .=,E12.4//)
360 1501 FORMAT(
361 & 5X,40HUNLOADING BULK. . . . . . . . . . . . .=,E12.4/,
362 & 5X,40HMAX VOLUMIC COMPRESSION . . . . . . . .=,E12.4/,
363 & 5X,40HFRACTURE PRESSURE . . . . . . . . . . .=,E12.4/,
364 & 5X,40HYIELD SURFACE PRESSURE ROOT . . . . . .=,E12.4//)
365
366 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)