33
34
35
36
37
38#include "implicit_f.inc"
39
40
41
42#include "units_c.inc"
43#include "comlock.inc"
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64 INTEGER NEL,NUVAR,ISMSTR
65 INTEGER NGL(NEL)
67 my_real uparam(*),epsxx(nel) ,epsyy(nel),epsxy(nel),dfmax(nel)
68
69
70
71 my_real uvar(nel,nuvar), off(nel),offl(nel),
72 . signxx(nel),signyy(nel),signxy(nel)
73
74
75
76
77
78 INTEGER I,J,LENG,NINDX,INDX(NEL),TYPE_MAX,F_FLAG,STRDEF,STRFLAG
79 my_real e11,sig_a,sig_b,sig_11,eps_a,eps_b,
80 . c_min,c_max,ema,damage
81 DOUBLE PRECISION :: A0(2),A1(2),A2(2),C0,C1,C2,C3,C4,C5,C6,C7,C8,C9,
82 . X1,X2,X3,Y1,Y2,Y3,Z1,Z2,Z3,F,FF,D,DD,D2,DP,E,G
83 DOUBLE PRECISION, PARAMETER :: PI8 = 0.3926990817d0
84 DOUBLE PRECISION, PARAMETER :: PI38 = 1.178097245d0
85 DOUBLE PRECISION, PARAMETER :: SPI8 = 0.3826834324d0
86 DOUBLE PRECISION, PARAMETER :: SPI38 = 0.9238795325d0
87
88
89
90
91
92
93
94
95
96 type_max = int(uparam(1))
97 c_min = uparam(2)
98 c_max = uparam(3)
99 ema = uparam(4)
100 ff = uparam(5)
101 f_flag = int(uparam(6))
102 strdef = int(uparam(7))
103
104 nindx = 0
105
106
107
108 strflag = 0
109 IF (strdef == 2) THEN
110 IF (ismstr == 10 .or. ismstr == 12) THEN
111 strflag = 1
112 ELSE IF (ismstr == 0 .or. ismstr == 2 .or. ismstr == 4) THEN
113 strflag = 2
114 END IF
115 ELSE IF (strdef == 3) THEN
116 IF (ismstr == 1 .or. ismstr == 3 .or. ismstr == 11) THEN
117 strflag = 3
118 ELSE IF (ismstr == 10 .or. ismstr == 12) THEN
119 strflag = 4
120 END IF
121 END IF
122
123
124
125
126
127 DO i=1,nel
128 IF (uvar(i,1) < 1.0d0 .AND. off(i) == one)THEN
129
130 IF (type_max == 1) THEN
131
132
133 sig_a = (signxx(i) + signyy(i))/two
134 sig_b = sqrt(((signxx(i)-signyy(i))/two)**2+signxy(i)**2)
135 sig_11 = sig_a + sig_b
136 e11 = sig_11
137
138 ELSE
139
140
141 eps_a = (epsxx(i) + epsyy(i)) * half
142 eps_b = sqrt(((epsxx(i)-epsyy(i))/two)**2+(half*epsxy(i))**2)
143 e11 = eps_a + eps_b
144 IF (strflag == 1) THEN
145 e11 = sqrt(e11 + one) - one
146 ELSE IF (strflag == 2) THEN
147 e11 = exp(e11) - one
148 ELSE IF (strflag == 3) THEN
149 e11 = log(e11 + one)
150 ELSE IF (strflag == 4) THEN
151 e11 = log(sqrt(e11+one))
152 END IF
153 ENDIF
154
155
156
157
158 IF (ema == one .and. ff /= zero .and. f_flag > 1) THEN
159
160
161
162 f =
min(ff,zep4 /
max(timestep,em20))
163 d = tan(pi*f*
max(timestep,em20))
164 dd = d*d
165 d2 = two*d
166 dp = one + dd
167 e = d2*spi8
168 g = e + dp
169 g = one/g
170
171 c0 = dd * g
172 c1 = two* c0
173 c2 = c0
174 c3 = two * g - c1
175 c4 = (e - dp) * g
176
177 e = d2*spi38
178 g = e + dp
179 g = one/g
180
181 c5 = dd * g
182 c6 = two * c5
183 c7 = c5
184 c8 = two * g - c6
185 c9 = (e - dp) * g
186
187
188
189
190
191 a0(1) = uvar(i,3)*uvar(i,9)
192 a0(2) = uvar(i,4)*uvar(i,9)
193 a1(1) = uvar(i,5)*uvar(i,10)
194 a1(2) = uvar(i,6)*uvar(i,10)
195 a2(1) = uvar(i,7)*uvar(i,11)
196 a2(2) = uvar(i,8)*uvar(i,11)
197
198 x1 = a0(2)
199 x2 = a0(1)
200
201 x3 = e11
202 y1 = a1(2)
203 y2 = a1(1)
204 y3 = c0 * x3
205 y3 = y3 + c1 * x2
206 y3 = y3 + c2 * x1
207 y3 = y3 + c3 * y2
208 y3 = y3 + c4 * y1
209 z1 = a2(2)
210 z2 = a2(1)
211 z3 = c5 * y3
212 z3 = z3 + c6 * y2
213 z3 = z3 + c7 * y1
214 z3 = z3 + c8 * z2
215 z3 = z3 + c9 * z1
216
217 a0(2) = x2
218 a0(1) = x3
219 a1(2) = y2
220 a1(1) = y3
221 a2(2) = z2
222 a2(1) = z3
223
224 IF ((x3 /= zero).AND.(x2 /= zero)) THEN
225 uvar(i,3) = a0(1)/x2
226 uvar(i,4) = a0(2)/x2
227 uvar(i,5) = a1(1)/(c0*x3)
228 uvar(i,6) = a1(2)/(c0*x3)
229 uvar(i,7) = a2(1)/(c0*y3)
230 uvar(i,8) = a2(2)/(c0*y3)
231 uvar(i,9) = x2
232 uvar(i,10) = c0*x3
233 uvar(i,11) = c0*y3
234 ELSE
235 uvar(i,3) = a0(1)
236 uvar(i,4) = a0(2)
237 uvar(i,5) = a1(1)
238 uvar(i,6) = a1(2)
239 uvar(i,7) = a2(1)
240 uvar(i,8) = a2(2)
241 uvar(i,9) = one
242 uvar(i,10) = one
243 uvar(i,11) = one
244 ENDIF
245
246 e11
247
248 ELSE
249
250
251
252
253 e11 = ema * e11 + ( one - ema ) * uvar(i,2)
254 uvar(i,2) = e11
255 ENDIF
256
257
258
259
260
261
262
263
264
265
266
267 damage =
max(zero ,
min(one ,(e11-c_min)/
max(em6,(c_max-c_min)) ))
268 uvar(i,1) =
max(uvar(i,1),damage)
269 dfmax(i) = uvar(i,1)
270
271 IF (uvar(i,1) >= one) THEN
272 nindx = nindx+1
273 indx(nindx) = i
274 ENDIF
275
276 ENDIF
277
278
279
280
281
282 ENDDO
283
284 DO j=1,nindx
285 i = indx(j)
286#include "lockon.inc"
287 WRITE(iout, 1000) ngl(i),time
288 WRITE(istdo,1100) ngl(i),time
289#include "lockoff.inc"
290
291 ENDDO
292
293 1000 FORMAT(1x,'SHELL ELEMENT NUMBER (VISUAL) el#',i10,
294 . ' LIMIT REACHED AT TIME :',1pe12.4)
295 1100 FORMAT(1x,'SHELL ELEMENT NUMBER (VISUAL) el#',i10,
296 . ' LIMIT REACHED AT TIME :',1pe12.4)
297
298 RETURN