35
36
37
38#include "implicit_f.inc"
39#include "impl1_c.inc"
40#include "comlock.inc"
41
42
43
44#include "mvsiz_p.inc"
45
46
47
48#include "param_c.inc"
49#include "com08_c.inc"
50#include "units_c.inc"
51#include "scr17_c.inc"
52
53
54
55 INTEGER JFT, JLT,IMAT,IPLA,NEL
56
58 . pm(npropm,*), off(*), dir(nel,2), epseq(*),
59 . ezz(*),dt1c(*),dpla1(*),epspd(*),hardm(*),
60 . signxx(nel),signyy(nel),signxy(nel),signyz(nel),signzx(nel),
61 . depsxx(mvsiz),depsyy(mvsiz),depsxy(mvsiz),depsyz(mvsiz),
62 . depszx(mvsiz)
63
64
65
66 INTEGER ICC
67 INTEGER I,NMAX,J, NINDX, INDEX(MVSIZ),NGL(MVSIZ),
68 . INDX(MVSIZ),N
69
71 . e,ca, ce, cn, nu,
72 .
ymax,cm,epdr,epchk(mvsiz),
73 . a11,a22,a1122,a12,
74 . seq,yld,scale, epsp,s11,s22,s12,g3(mvsiz),nu1,
75 . d11, d22, d12, dpla, dtinv, p,a,b,c,a1s1,a2s2,q,umr,
76 . s_11(mvsiz),s_22(mvsiz),s_12(mvsiz),s1,s2,s3,seqh(mvsiz),
77 . h(mvsiz),dpla_i(mvsiz),dpla_j(mvsiz),etse(mvsiz),nu2(mvsiz),
78 . nu3(mvsiz),nu4(mvsiz),nu5(mvsiz),nnu1(mvsiz),q12(mvsiz),
79 . q21(mvsiz),jq(mvsiz),jq2(mvsiz),a_1,a_2,a_3,
80 . st11(mvsiz),st22(mvsiz),st12(mvsiz),axx(mvsiz),ayy(mvsiz),
81 . a_xy(mvsiz),axy(mvsiz),b_1(mvsiz),b_2(mvsiz),b_3(mvsiz),
82 . yldp(mvsiz),dr,p_1,p_2,p_3,pp1,pp2,pp3,f,df,yld_i,a1,
83 . a2,epsmax,sige(mvsiz,5)
84 DATA nmax/3/
85
86 e = pm(20,imat)
87 nu = pm(21,imat)
88 ca = pm(38,imat)
89 ce = pm(39,imat)
90 cn = pm(40,imat)
92 cm = pm(43,imat)
93 epdr = pm(44,imat)
94 a11 = pm(45,imat)
95 a22 = pm(46,imat)
96 a1122= pm(47,imat)
97 a12 = pm(48,imat)
98 icc = nint(pm(49,imat))
99 a1 = pm(24,imat)
100 a2 = pm(25,imat)
101 epsmax = pm(41,imat)
102
103 IF (ipla == 0) THEN
104
105
106
107
108 DO i=jft,jlt
109 d11 = dir(i,1)*dir(i,1)
110 d22 = dir(i,2)*dir(i,2)
111 d12 = dir(i,1)*dir(i,2)
112 s11 = d11*signxx(i) + d22*signyy(i) + two*d12*signxy(i)
113 s22 = d22*signxx(i) + d11*signyy(i) - two*d12*signxy(i)
114 s12 = d12*(signyy(i) - signxx(i)) + ( d11 - d22 )*signxy(i)
115
116
117
118 epsp =
max(abs(depsxx(i)),abs(depsyy(i)),half*abs(depsxy(i)))
119 dtinv= dt1c(i)/
max(dt1c(i)**2,em20)
120 epsp = epsp*dtinv
121 epsp =
max(epsp ,epdr)
122 epspd(i) = epsp
123 nnu1(i) = nu / (one - nu)
124 nu5(i) = one-nnu1(i)
125
126
127
128 yld = ca*(ce+epseq(i))**cn*epsp**cm
130
131
132
133 seq = sqrt(a11*s11*s11 + a22*s22*s22
134 . -a1122*s11*s22 + a12*s12*s12)
135 scale =
min(one,yld /
max(seq ,em20))
136 signxx(i)=signxx(i)*scale
137 signyy(i)=signyy(i)*scale
138 signxy(i)=signxy(i)*scale
139 dpla = off(i)*
max(zero,(seq -yld )/e)
140 ezz(i) = - nu5(i)*dpla*half*(signxx(i)+signyy(i))/
max(em20,yld)
141 epseq(i) = epseq(i) + dpla
142 dpla1(i) = dpla
143 ENDDO
144 ELSEIF (ipla == 2) THEN
145
146
147
148
149 DO i=jft,jlt
150 d11 = dir(i,1)*dir(i,1)
151 d22 = dir
152 d12 = dir(i,1)*dir(i,2)
153 s11 = d11*signxx(i) + d22*signyy(i) + two*d12*signxy(i)
154 s22 = d22*signxx(i) + d11*signyy(i) - two*d12*signxy(i)
155 s12 = d12*(signyy(i) - signxx(i)) + ( d11 - d22 )*signxy(i)
156
157
158
159 epsp =
max(abs(depsxx(i)),abs(depsyy(i)),half*abs(depsxy(i)))
160 dtinv= dt1c(i)/
max(dt1c(i)**2,em20)
161 epsp = epsp*dtinv
162 epsp =
max(epsp ,epdr)
163 epspd(i) = epsp
164
165
166
167 yld =ca*(ce+epseq(i))**cn*epsp**cm
169
170
171
172 nu1 = nu/(one - nu)
173 nu5(i) = one-nu1
174 g3(i) = three_half*e/(one + nu)
175 p = -(s11+s22) * third
176 q = (one -nu1)*p
177 s11 = s11+q
178 s22 = s22+q
179 a1s1 = a11*s11
180 a2s2 = a22*s22
181 a = a1s1*s11 + a2s2*s22 - a1122*s11*s22 + a12*s12*s12
182 b = -q*(a1s1 + a2s2 - half*a1122*(s11+s22))
183 c = (a11+a22-a1122)*q*q
184 seq = sqrt(a+b+b+c)
185 c = c - yld*yld
186 scale =
min(one,(-b + sqrt(
max(zero,b*b-a*c)))/
max(a ,em20))
187 umr = one - scale
188 q = q*umr
189 signxx(i) = signxx(i)*scale - q
190 signyy(i) = signyy(i)*scale - q
191 signxy(i) = signxy(i)*scale
192 dpla = off(i)*seq*umr/
max(em20,g3(i))
193 ezz(i) = -nu5(i)*dpla*half*(signxx(i)+signyy(i)) / yld
194 epseq(i) = epseq(i) + dpla
195 dpla1(i) = dpla
196 ENDDO
197
198 ELSEIF (ipla == 1) THEN
199
200
201
202 DO i=jft,jlt
203 d11 = dir(i,1)*dir(i,1)
204 d22 = dir(i,2)*dir(i,2)
205 d12 = dir(i,1)*dir(i,2)
206 s_11(i) = d11*signxx(i) + d22*signyy(i) + two*d12*signxy(i)
207 s_22(i) = d22*signxx(i) + d11*signyy(i) - two*d12*signxy(i)
208 s_12(i) = d12*(signyy(i) - signxx(i)) + ( d11 - d22 )*signxy(i)
209 nnu1(i) = nu / (one - nu)
210
211
212
213 epsp =
max(abs(depsxx(i)),abs(depsyy(i)),half*abs(depsxy(i)))
214 dtinv= dt1c(i)/
max(dt1c(i)**2,em20)
215 epsp = epsp*dtinv
216 epsp =
max(epsp ,epdr)
217 epspd(i) = epsp
218
219
220
221 s1=a11*s_11(i)*s_11(i)
222 s2=a22*s_22(i)*s_22(i)
223 s3=a1122*s_11(i)*s_22(i)
224 axy(i)=a12*s_12(i)*s_12(i)
225 seqh(i)=sqrt(s1+s2-s3+axy(i))
226 ezz(i) = -(depsxx(i)+depsyy(i))*nnu1(i)
227 g3(i) = three_half*e/(one+nu)
228 yldp(i) = ca*(ce+epseq(i))**cn*epsp**cm
230 IF (yldp(i) >=
ymax)
THEN
231 h(i)=zero
232 ELSE
233 h(i)=ca*cn*(ce+epseq(i))**(cn-one)*epsp**cm
234 ENDIF
235 ENDDO
236
237
238
239 DO i=jft,jlt
240 hardm(i) = h(i)
241 ENDDO
242
243
244
245 nindx=0
246 DO i=jft,jlt
247 IF (seqh(i) > yldp(i) .AND. off(i) == one) THEN
248 nindx=nindx+1
249 index(nindx)=i
250 ENDIF
251 ENDDO
252 IF (nindx > 0) THEN
253
254
255
256#include "vectorize.inc"
257 DO j=1,nindx
258 i=index(j)
259 dpla_j(i)= (seqh(i)-yldp(i))/(g3(i)+h(i))
260 etse(i)= h(i)/(h(i)+e)
261 nu2(i) = one-nu*nu
262 nu3(i) = nu*half
263 nu4(i) = half*(one-nu)
264 nu5(i) = one-nnu1(i)
265 s1=a11*nu*two-a1122
266 s2=a22*nu*two-a1122
267 s12=a1122-nu*(a11+a22)
268 s3=sqrt(nu2(i)*(a11-a22)*(a11-a22)+s12*s12)
269 IF (abs(s1) < em20) THEN
270 q12(i)=zero
271 ELSE
272 q12(i)=-(a11-a22+s3)/s1
273 ENDIF
274 IF (abs(s2) < em20) THEN
275 q21(i)=zero
276 ELSE
277 q21(i)=(a11-a22+s3)/s2
278 ENDIF
279 jq(i)=one/(one-q12(i)*q21(i))
280 jq2(i)=jq(i)*jq(i)
281 a=a11*q12(i)
282 b=a22*q21(i)
283 a_1=(a11+a1122*q21(i)+b*q21(i))*jq2(i)
284 a_2=(a22+a1122*q12(i)+a*q12(i))*jq2(i)
285 a_3=(a+b)*jq2(i)*two+a1122*(jq2(i)*two-jq(i))
286 st11(i)=s_11(i)+s_22(i)*q12(i)
287 st22(i)=q21(i)*s_11(i)+s_22(i)
288 axx(i)=a_1*st11(i)*st11(i)
289 ayy(i)=a_2*st22(i)*st22(i)
290 a_xy(i)=a_3*st11(i)*st22(i)
291 a=a1122*nu3(i)
292 b=s3*jq(i)
293 b_1(i)=a22-a-b
294 b_2(i)=a11-a+b
295 b_3(i)=a12*nu4(i)
296 ENDDO
297
298 DO n=1,nmax
299#include "vectorize.inc"
300 DO j=1,nindx
301 i=index(j)
302 dpla_i(i)=dpla_j(i)
303 IF (dpla_i(i) > zero) THEN
304 yld_i =
min(yldp(i)+h(i)*dpla_i(i),
ymax)
305 dr =a1*dpla_i(i)/yld_i
306 p_1=one/(one+b_1(i)*dr)
307 pp1=p_1*p_1
308 p_2=one/(one+b_2(i)*dr)
309 pp2=p_2*p_2
310 p_3=one/(one+b_3(i)*dr)
311 pp3=p_3*p_3
312 f = axx(i)*pp1+ayy(i)*pp2-a_xy(i)*p_1*p_2+axy(i)*pp3
313 . -yld_i*yld_i
314 df = -((axx(i)*p_1-a_xy(i)*p_2*half)*pp1*b_1(i)+
315 . (ayy(i)*p_2-a_xy(i)*p_1*half)*pp2*b_2(i)+
316 . axy(i)*pp3*p_3*b_3(i))*(a1-dr*h(i))/yld_i
317 . -h(i)*yld_i
318 dpla_j(i)=
max(zero,dpla_i(i)-f*half/df)
319 ELSE
320 dpla_j(i)=zero
321 ENDIF
322 ENDDO
323 ENDDO
324
325
326
327#include "vectorize.inc"
328 DO j=1,nindx
329 i=index(j)
330 epseq(i) = epseq(i) + dpla_j(i)
331 yld_i =yldp(i)+h(i)*dpla_j(i)
333 dr =a1*dpla_j(i)/yld_i
334 p_1=one/(one+b_1(i)*dr)
335 p_2=one/(one+b_2(i)*dr)
336 p_3=one/(one+b_3(i)*dr)
337 s1=st11(i)*p_1
338 s2=st22(i)*p_2
339 s_11(i)=jq(i)*(s1-s2*q12(i))
340 s_22(i)=jq(i)*(s2-s1*q21(i))
341 s_12(i)=s_12(i)*p_3
342 s1=a11*s_11(i)+a22*s_22(i)
343 . -a1122*(s_11(i)+s_22(i))*half
344 ezz(i) = - nu5(i)*dpla_j(i)*s1/yld_i
345 s1 =a1122*half
346 p_1=a11*s_11(i)-s1*s_22(i)
347 p_2=a22*s_22(i)-s1*s_11(i)
348 p_3=a12*s_12(i)
349 dpla1(i) = dpla_j(i)
350 ENDDO
351
352 nindx=0
353 DO i=jft,jlt
354 IF (epseq(i) > epsmax .AND. off(i) == one) THEN
355 nindx=nindx+1
356 indx(nindx)=i
357 off(i)=zero
358 ENDIF
359 ENDDO
360
361 DO i=jft,jlt
362 sige(i,1) = s_11(i)
363 sige(i,2) = s_22(i)
364 sige(i,3) = s_12(i)
365 sige(i,4) = zero
366 sige(i,5) = zero
367 ENDDO
368
369 CALL urotov(jft,jlt,sige,dir,nel)
370
371 DO i=jft,jlt
372 signxx(i)= sige(i,1)
373 signyy(i)= sige(i,2)
374 signxy(i)= sige(i,3)
375 ENDDO
376 ENDIF
377
378 IF (nindx > 0) THEN
379 IF (imconv == 1) THEN
380 DO i = 1, nindx
381#include "lockon.inc"
382 WRITE(iout, 1000) ngl(indx(i))
383 WRITE(istdo,1100) ngl(indx(i)),tt
384#include "lockoff.inc"
385 ENDDO
386 ENDIF
387 ENDIF
388 ENDIF
389
390 1000 FORMAT(1x,'-- RUPTURE OF SHELL ELEMENT NUMBER ',i10)
391 1100 FORMAT(1x,'-- RUPTURE OF SHELL ELEMENT :',i10,' AT TIME :',g11.4)
392
393 RETURN
subroutine ymax(idn, fac, npc, pld, stiffmin, stiffmax, stiffini, stiffavg)
subroutine urotov(jft, jlt, sig, dir, nel)
1.15.0