38
39
40
41
42
43
44
45#include "implicit_f.inc"
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#include "units_c.inc"
78#include "scr17_c.inc"
79#include "comlock.inc"
80
81
82 INTEGER NEL, NUPARAM, NUVAR,NGL(NEL)
83
84 my_real time,timestep,uparam(nuparam),
85 . dpla(nel),epsp(nel),pla(nel),
86 . epspxx(nel),epspyy(nel),epspzz(nel),
87 . epspxy(nel),epspyz(nel),epspzx(nel),
88 . epsxx(nel) ,epsyy(nel) ,epszz(nel) ,
89 . epsxy(nel) ,epsyz(nel) ,epszx(nel) ,
90 . tdele(nel) ,dfmax(nel)
91
92
93
94 my_real uvar(nel,nuvar), off(nel),
95 . signxx(nel),signyy(nel),signzz(nel),
96 . signxy(nel),signyz(nel),signzx(nel)
97
98
99
100 INTEGER I,J,LENG,NINDX,INDX(NEL),NCYC
101 my_real c0,ema,e_hyd,e_11,e_22,e_33,e_12,e_23,
102 . e_13,eeq,d_eeq
104 . q,r,r_inter,phi,lockroft
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126 c0 = uparam(1)
127 ema = uparam(2)
128
129 lockroft = zero
130 e_hyd = zero
131 eeq = zero
132 d_eeq = zero
133 r_inter = zero
134 phi = zero
135
136 leng = 0
137
138
139
140 nindx = 0
141
142
143
144 DO i=1,nel
145 IF(off(i)<one .AND. off(i) > zero) THEN
146 signxx(i) = signxx(i) * off(i)
147 signyy(i) = signyy(i) * off(i)
148 signzz(i) = signzz(i) * off(i)
149 signxy(i) = signxy(i) * off(i)
150 signyz(i) = signyz(i) * off(i)
151 signzx(i) = signzx(i) * off(i)
152 off(i) = off(i) * four_over_5
153 ENDIF
154 IF(off(i)<em01) off(i) = zero
155 END DO
156
157
158
159 DO i=1,nel
160 IF(off(i)==one) THEN
161
162
163
164 IF(c0 < zero)THEN
165 e_hyd = third * (epsxx(i) + epsyy(i) + epszz(i))
166 e_11 = epsxx(i) - e_hyd
167 e_22 = epsyy(i) - e_hyd
168 e_33 = epszz(i) - e_hyd
169 e_12 = half*epsxy(i)
170 e_23 = half*epsyz(i)
171 e_13 = half*epszx(i)
172
173 eeq = e_11**2 + e_22**2 + e_33**2
174 eeq = eeq + two * (e_12**2) + two * (e_23**2) + two * (e_13
175 eeq = 0.8164965809 * sqrt(eeq)
176
177 d_eeq = eeq - uvar(i,1)
178 IF (d_eeq <= zero) d_eeq = zero
179 uvar(i,1) = eeq
180 ELSE
181 d_eeq = dpla(i)
182 uvar(i,1) = uvar(i,1) + dpla(i)
183 ENDIF
184
185
186
187 i1 = signxx(i)+signyy(i)+signzz(i)
188 i2 = signxx(i)*signyy(i)+signyy(i)*signzz(i)+signzz(i)*signxx(i)-
189 . signxy(i)*signxy(i)-signzx(i)*signzx(i)-signyz(i)*signyz(i
190 i3 = signxx(i)*signyy(i)*signzz(i)-signxx(i)*signyz(i)*signyz(i)-
191 . signyy(i)*signzx(i)*signzx(i)-signzz(i)*signxy(i)*signxy(i)+
192 . two*signxy(i)*signzx(i)*signyz(i)
193 q = (three*i2 - i1*i1)/nine
194 r = (two*i1*i1*i1-nine*i1*i2+twenty7*i3)/cinquante4
195
196 r_inter =
min(r/sqrt(
max(em20,(-q**3))),one)
197 phi = acos(
max(r_inter,-one))
198
199 s11 = two*sqrt(-q)*cos(phi/three)+third*i1
200 s22 = two*sqrt(-q)*cos((phi+two*pi)/three)+third*i1
201 s33 = two*sqrt(-q)*cos((phi+four*pi)/three)+third*i1
202
203 IF (s11<s22) THEN
204 r_inter = s11
205 s11 = s22
206 s22 = r_inter
207 ENDIF
208 IF (s22<s33)THEN
209 r_inter = s22
210 s22 = s33
211 s33 = r_inter
212 ENDIF
213 IF (s11<s22)THEN
214 r_inter = s11
215 s11 = s22
216 s22 = r_inter
217 ENDIF
218
219 IF(s11 > zero) THEN
220 s11 = s11 * ema + (one-ema)* uvar(i,3)
221 uvar(i,3) = s11
222
223 uvar(i,2) = uvar(i,2) +
max(s11,zero) * d_eeq
224 ENDIF
225
226 dfmax(i) = uvar(i,2) /
max(abs(c0),em20)
227 dfmax(i) =
min(one,dfmax(i))
228
229 IF (uvar(i,2) >= abs(c0)) THEN
230 tdele(i) = time
231 nindx = nindx+1
232 indx(nindx) = i
233
234 off(i) = four_over_5
235
236 ENDIF
237
238
239
240 ENDIF
241
242 ENDDO
243
244 DO j=1,nindx
245 i = indx(j)
246#include "lockon.inc"
247 WRITE(iout, 1000) ngl(i),time
248 WRITE(istdo,1100) ngl(i),time
249#include "lockoff.inc"
250
251 ENDDO
252
253 1000 FORMAT(1x,'DELETE SOLID ELEMENT NUMBER (COCKROFT-LATHAM) el#',i10,
254 . ' AT TIME :',1pe12.4)
255 1100 FORMAT(1x,'DELETE SOLID ELEMENT NUMBER (COCKROFT-LATHAM) el#',i10,
256 . ' AT TIME :',1pe12.4)
257
258 RETURN