-
Notifications
You must be signed in to change notification settings - Fork 0
/
umat41c_tmp.f
474 lines (472 loc) · 16.4 KB
/
umat41c_tmp.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c yield surface a and flow rules
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
subroutine murcia_a_flow(
& fcpr, hr, k3, k4, k5, sigma, tau, pp, pn, r, ma1)
double precision fcpr, hr, k3, k4, k5
double precision sigma, tau, pp, pn, r
double precision ma1
double precision rhosigma
rhosigma = max(0.0D+00,-sigma) / fcpr
ma1 = k3 * max(0.0D+00,1.0D+00-rhosigma)*exp(-k4*(pp+pn)/hr)
& - k5*rhosigma*max(0.0D+00,r)/hr
return
end
subroutine murcia_a_yield(
& fcpr, sr, c0, mua0, k1, k2, sigma, tau, pp, pn, fa)
double precision fcpr, sr, c0, mua0, k1, k2
double precision sigma, tau, pp, pn
double precision fa
double precision rhop, c, mua
rhop = (pp+pn) / sr
c = c0 * max(0.0D+00,1.0D+00-rhop)
mua = mua0 * exp(-k2*rhop)
fa = abs(sigma/fcpr)**k1 - (c/fcpr)**k1 + mua*tau/fcpr
return
end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c yield surface b and flow rules
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
subroutine murcia_angle(li, lt, alpha0, pp, pn, s, alpha)
c input arguments
double precision li, lt, alpha0
double precision pp, pn, s
c output arguments
double precision alpha
c function-scale parameters
double precision spp, spn
spp = s + pp
spn = s - pn
if (spp.le.-li) then
alpha = 0.0
elseif (spp.le.-li+lt) then
alpha = alpha0 * (spp-li) / lt
elseif (spp.le.-lt) then
alpha = -alpha0
elseif (spp.le.0.0) then
alpha = alpha0 * spp / lt
elseif (spn.le.0.0) then
alpha = 0.0
elseif (spn.le.lt) then
alpha = alpha0 * spn / lt
elseif (spn.le.li-lt) then
alpha = alpha0
elseif (spn.le.li) then
alpha = alpha0 * (li-spn) / lt
else
alpha = 0.0
endif
return
end
subroutine murcia_b_flow(li, lt, alpha0, pp, pn, s, mb1)
c input arguments
double precision li, lt, alpha0
double precision pp, pn, s
c output arguments
double precision mb1
c function-scale parameters
double precision alpha
call murcia_angle(li, lt, alpha0, pp, pn, s, alpha)
mb1 = -tan(alpha)
return
end
subroutine murcia_b_yield(
& li, lt, alpha0, mub, sigma, tau, pp, pn, s, fb)
c input arguments
double precision li, lt, alpha0, mub
double precision sigma, tau, pp, pn, s
c output arguments
double precision fb
c function-scale parameters
double precision alpha, sina, cosa
call murcia_angle(li, lt, alpha0, pp, pn, s, alpha)
sina = sin(alpha)
cosa = cos(alpha)
fb = abs(tau*cosa+sigma*sina) + mub*(sigma*cosa-tau*sina)
return
end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c return to a
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
subroutine murcia_a_f2(n, x, fvec, iflag, cm, lcm, ext, lext)
integer n, iflag, lcm, lext
double precision x(n), fvec(n), cm(lcm), ext(lext)
double precision ma1t, ma1t1, ddns
call murcia_a_flow(
& cm(3), cm(5), cm(16), cm(17), cm(18), ext(3), ext(4), ext(5),
& ext(6), ext(7), ma1t)
call murcia_a_flow(
& cm(3), cm(5), cm(16), cm(17), cm(18), x(1), ext(10), ext(11),
& ext(11), ext(12), x(2), ma1t1)
ddns = ext(18) * ((1.0D+00-cm(19))*ma1t+cm(19)*ma1t1)
fvec(1) = ext(1) - cm(12)*ddns - ext(9)
fvec(2) = ext(7) + ddns - ext(13)
return
end
subroutine murcia_a_f1(n, x, fvec, iflag, cm, lcm, ext, lext)
integer n, iflag, lcm, lext
double precision x(n), fvec(n), cm(lcm), ext(lext)
external murcia_a_f2
integer n2, lwa2
parameter (n2 = 2, lwa2 = (n2*(3*n2+13)) / 2)
double precision x2(n2), fvec2(n2), info2, wa2(lwa2)
double precision tol2/0.00001D+00/
ext(10) = ext(2) - x(1)*cm(13)*ext(15)
ext(11) = ext(5) + x(1)*max(0.0D+00,ext(15))
ext(12) = ext(6) + x(1)*max(0.0D+00,-ext(15))
ext(18) = x(1)
x2(1) = ext(3)
x2(2) = ext(7)
call hybrdext1(
& murcia_a_f2, n2, x2, fvec2, tol2, info2, wa2, lwa2, cm, lcm,
& ext, lext)
ext(9) = x2(1)
ext(13) = x2(2)
call murcia_a_yield(
& cm(3), cm(4), cm(9), cm(10), cm(14), cm(15), ext(9), ext(10),
& ext(11), ext(12), fvec(1))
return
end
subroutine murcia_a_return(
& cm, lcm, ext, lext, sigmat1, taut1, ppt1, pnt1, rt1, st1)
integer lcm, lext
double precision cm(lcm), ext(lext)
external murcia_a_f1
integer n, lwa
parameter (n = 1, lwa = (n*(3*n+13)) / 2)
double precision x(n), fvec(n), info, wa(lwa)
double precision tol/0.00001D+00/
ext(15) = sign(1.0D+00, ext(2))
x(1) = 0.0D+00
call hybrdext1(
& murcia_a_f1, n, x, fvec, tol, info, wa, lwa, cm, lcm,
& ext, lext)
ext(14) = ext(8)
return
end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c return to b
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
subroutine murcia_b_f1(n, x, fvec, iflag, cm, lcm, ext, lext)
integer n, iflag, lcm, lext
double precision x(n), fvec(n), cm(lcm), ext(lext)
double precision mb1t, mb1t1
ext(10) = ext(2) - x(1)*cm(13)*ext(16)
ext(14) = ext(8) + x(1)*ext(16)
call murcia_b_flow(
& cm(6), cm(7), cm(8), ext(5), ext(6), ext(8), mb1t)
call murcia_b_flow(
& cm(6), cm(7), cm(8), ext(5), ext(6), ext(14), mb1t1)
ext(9) =
& ext(1) - x(1)*cm(12)*((1.0D+00-cm(19))*mb1t+cm(19)*mb1t1)
& *ext(16)
call murcia_b_yield(
& cm(6), cm(7), cm(8), cm(11), ext(9), ext(10), ext(5), ext(6),
& ext(14), fvec(1))
return
end
subroutine murcia_b_return(
& cm, lcm, ext, lext, sigmat1, taut1, ppt1, pnt1, rt1, st1)
integer lcm, lext
double precision cm(lcm), ext(lext)
external murcia_a_f1
integer n, lwa
parameter (n = 1, lwa = (n*(3*n+13)) / 2)
double precision x(n), fvec(n), info, wa(lwa)
double precision tol/0.00001D+00/
if (ext(17).le.0.0D+00) then
ext(9) = 0.0D+00
ext(10) = 0.0D+00
ext(11) = ext(5)
ext(12) = ext(6)
ext(13) = ext(7)
ext(14) = ext(8) + ext(2)/cm(13)
return
endif
x(1) = 0.0D+00
call hybrdext1(
& murcia_b_f1, n, x, fvec, tol, info, wa, lwa, cm, lcm,
& ext, lext)
ext(11) = ext(5)
ext(12) = ext(6)
ext(13) = ext(7)
return
end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c return to the intersection a and b
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
subroutine murcia_ab_f2(n, x, fvec, iflag, cm, lcm, ext, lext)
integer n, iflag, lcm, lext
double precision x(n), fvec(n), cm(lcm), ext(lext)
double precision ma1t, ma1t1, ddns, mb1t, mb1t1
call murcia_a_flow(
& cm(3), cm(5), cm(16), cm(17), cm(18), ext(3), ext(4), ext(5),
& ext(6), ext(7), ma1t)
call murcia_a_flow(
& cm(3), cm(5), cm(16), cm(17), cm(18), x(1), ext(10), ext(11),
& ext(11), ext(12), x(2), ma1t1)
ddns = ext(18) * ((1.0D+00-cm(19))*ma1t+cm(19)*ma1t1)
call murcia_b_flow(
& cm(6), cm(7), cm(8), ext(5), ext(6), ext(8), mb1t)
call murcia_b_flow(
& cm(6), cm(7), cm(8), ext(5), ext(6), ext(14), mb1t1)
fvec(1) = ext(1) - cm(12)*ddns - ext(9)
& - ext(19)*cm(12)*((1.0D+00-cm(19))*mb1t+cm(19)*mb1t1)*ext(16)
fvec(2) = ext(7) + ddns - ext(13)
return
end
subroutine murcia_ab_f1(n, x, fvec, iflag, cm, lcm, ext, lext)
integer n, iflag, lcm, lext
double precision x(n), fvec(n), cm(lcm), ext(lext)
external murcia_ab_f2
integer n2, lwa2
parameter (n2 = 2, lwa2 = (n2*(3*n2+13)) / 2)
double precision x2(n2), fvec2(n2), info2, wa2(lwa2)
double precision tol2/0.00001D+00/
ext(10) = ext(2) - x(1)*cm(13)*ext(15) - x(2)*cm(11)*ext(16)
ext(11) = ext(5) + x(1)*max(0.0D+00,ext(15))
ext(12) = ext(6) + x(1)*max(0.0D+00,-ext(15))
ext(14) = ext(8) + x(2)*ext(16)
ext(18) = x(1)
ext(19) = x(2)
x2(1) = ext(3)
x2(2) = ext(7)
call hybrdext1(
& murcia_ab_f2, n2, x2, fvec2, tol2, info2, wa2, lwa2, cm, lcm,
& ext, lext)
ext(9) = x2(1)
ext(13) = x2(2)
call murcia_a_yield(
& cm(3), cm(4), cm(9), cm(10), cm(14), cm(15), ext(9), ext(10),
& ext(11), ext(12), fvec(1))
call murcia_b_yield(
& cm(6), cm(7), cm(8), cm(11), ext(9), ext(10), ext(5), ext(6),
& ext(14), fvec(2))
return
end
subroutine murcia_ab_return(cm, lcm, ext, lext)
integer lcm, lext
double precision cm(lcm), ext(lext)
external murcia_ab_f1
integer n, lwa
parameter (n = 2, lwa = (n*(3*n+13)) / 2)
double precision x(n), fvec(n), info, wa(lwa)
double precision tol/0.00001D+00/
x(1) = 0.0D+00
x(2) = 0.0D+00
call hybrdext1(
& murcia_ab_f1, n, x, fvec, tol, info, wa, lwa, cm, lcm,
& ext, lext)
return
end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c return to the intersection a and b
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
subroutine murcia_return(cm, lcm, ext, lext)
integer lcm, lext
double precision cm(lcm), ext(lext)
double precision alphat
double precision fa, fb, swa
double precision sigmat1a, taut1a, ppt1a, pnt1a, rt1a, st1a
double precision sigmat1b, taut1b, ppt1b, pnt1b, rt1b, st1b
double precision fab, fba
call murcia_a_yield(
& cm(3), cm(4), cm(9), cm(10), cm(14), cm(15), ext(1), ext(2),
& ext(11), ext(12), fa)
call murcia_b_yield(
& cm(6), cm(7), cm(8), cm(11), ext(9), ext(10), ext(1), ext(2),
& ext(14), fb)
if (fa.le.0.0D+00.and.fb.le.0.0D+00) then
ext(9) = ext(3)
ext(10) = ext(4)
ext(11) = ext(5)
ext(12) = ext(6)
ext(13) = ext(7)
ext(14) = ext(8)
return
endif
ext(15) = sign(1.0D+00, ext(2))
call murcia_angle(
& cm(6), cm(7), cm(8), ext(5), ext(6), ext(8), alphat)
ext(16) = sign(1.0D+00, ext(2)*cos(alphat)+ext(1)*sin(alphat))
swa =
& abs(ext(2))/cm(13)*cm(16)*exp(-cm(17)*(ext(5)+ext(6))/cm(5))
& - ext(1)/cm(12)
ext(17) = ext(2)/cm(13)*tan(alphat) - ext(1)/cm(12)
if (fa.gt.0.0D+00.and.swa.gt.0.0D+00.and.fb.gt.0.0D+00) then
call murcia_a_return(cm, lcm, ext, lext)
sigmat1a = ext(9)
taut1a = ext(10)
ppt1a = ext(11)
pnt1a = ext(12)
rt1a = ext(13)
st1a = ext(14)
call murcia_b_return(cm, lcm, ext, lext)
sigmat1b = ext(9)
taut1b = ext(10)
ppt1b = ext(11)
pnt1b = ext(12)
rt1b = ext(13)
st1b = ext(14)
call murcia_a_yield(
& cm(3), cm(4), cm(9), cm(10), cm(14), cm(15), sigmat1b, taut1b,
& ppt1b, pnt1b, fab)
call murcia_b_yield(
& cm(6), cm(7), cm(8), cm(11), sigmat1a, taut1a, ppt1a, pnt1a,
& st1a, fba)
if (fab.gt.0.0D+00) then
if (fba.gt.0.0D+00) then
call murcia_ab_return(cm, lcm, ext, lext)
else
ext(9) = sigmat1a
ext(10) = taut1a
ext(11) = ppt1a
ext(12) = pnt1a
ext(13) = rt1a
ext(14) = st1a
endif
else
if (fba.gt.0.0D+00) then
ext(9) = sigmat1b
ext(10) = taut1b
ext(11) = ppt1b
ext(12) = pnt1b
ext(13) = rt1b
ext(14) = st1b
else
call murcia_ab_return(cm, lcm, ext, lext)
endif
endif
return
endif
if (fa.gt.0.0D+00.and.swa.gt.0.0D+00) then
call murcia_a_return(cm, lcm, ext, lext)
call murcia_b_yield(
& cm(6), cm(7), cm(8), cm(11), ext(9), ext(10), ext(5), ext(6),
& ext(14), fba)
if (fba.gt.0.0D+00) then
call murcia_ab_return(cm, lcm, ext, lext)
endif
return
endif
call murcia_b_return(cm, lcm, ext, lext)
call murcia_a_yield(
& cm(3), cm(4), cm(9), cm(10), cm(14), cm(15), ext(9), ext(10),
& ext(11), ext(12), fab)
if (fab.gt.0.0D+00) then
call murcia_ab_return(cm, lcm, ext, lext)
endif
return
end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c usermat definition
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
subroutine umat41c(
& idpart, cm, lft, llt, fc, dx, dxdt, aux, ek, ifail, dt1siz, crv,
& nhxbwp,cma)
c
c******************************************************************
c| Livermore Software Technology Corporation (LSTC) |
c| ------------------------------------------------------------ |
c| Copyright 1987-2008 Livermore Software Tech. Corp |
c| All rights reserved |
c******************************************************************
c
c vectorized cohesive material user model example
c
c*** variables
c idpart ---- part ID
c cm -------- material constants
c lft,llt --- start and end of block
c fc -------- components of the cohesive force
c dx -------- components of the displacement
c dxdt ------ components of the velocity
c aux ------- history storage
c ek -------- max. stiffness/area for time step calculation
c ifail ----- =.false. not failed
c =.true. failed
c dt1siz ---- time step size
c crv ------- curve array
c nhxbwp ---- internal element id array, lqfinv(nhxbwp(i),2)
c gives external element id
c cma=additional memory for material data defined by LMCA at
c 6th field of 2nd crad of *DATA_USER_DEFINED
c
c*** dx, dxdt, and fc are in the local coordinate system:
c components 1 and 2 are in the plane of the cohesive surface
c component 3 is normal to the plane
c
c*** cm storage convention
c (1) =0 density is per area
c =1 density is per volume
c (2) number of integration points for element deletion
c =0 no deletion
c (3:48) material model constants
c
include 'nlqparm'
logical ifail
dimension cm(*), fc(nlq,*), dx(nlq,*), dxdt(nlq,*), aux(nlq,*),
& ek(*), ifail(*), dt1siz(*), crv(lq1,2,*), nhxbwp(*), cma(*)
integer lext, lcm
parameter (lext = 19, lcm = 19)
double precision ext(lext)
c
! et=cm(3)
! en=cm(4)
! eki=max(et,en)
! fcfail=cm(5)
c
do i=lft,llt
ext(1) = aux(i,1) + cm(12)*dx(i,3)
ext(2) = aux(i,2) + cm(13)*dx(i,1)
ext(3) = aux(i,1)
ext(4) = aux(i,2)
ext(5) = aux(i,3)
ext(6) = aux(i,4)
ext(7) = aux(i,5)
ext(8) = aux(i,6)
call murcia_return(cm, lcm, ext, lext)
fc(i,3) = ext(9)
fc(i,1) = ext(10)
aux(i,1) = ext(9)
aux(i,2) = ext(10)
aux(i,3) = ext(11)
aux(i,4) = ext(12)
aux(i,5) = ext(13)
c fc(i,1)=et*dx(i,1)
c fc(i,2)=et*dx(i,2)
c fc(i,3)=en*dx(i,3)
c ek(i)=eki
c ifail(i)=fc(i,3).gt.fcfail
enddo
c
return
end