36
37
38
39#include "implicit_f.inc"
40
41
42
43#include "units_c.inc"
44#include "comlock.inc"
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 INTEGER ,INTENT(IN) :: NEL,NUPARAM,NUVAR,ILAY,IPT,NPTOT,IPTT
87 INTEGER ,DIMENSION(NEL), INTENT(IN) :: NGL
88 my_real ,
DIMENSION(NUPARAM),
INTENT(IN) :: uparam
89 my_real ,
DIMENSION(NEL),
INTENT(IN) :: dpla,pla
90
91
92
93
94
95 INTEGER ,DIMENSION(NEL), INTENT(INOUT) :: FOFF,FLD_IDX
97 my_real ,
DIMENSION(NEL),
INTENT(IN) :: signxx,signyy,signxy,signyz,signzx
98 my_real ,
DIMENSION(NEL),
INTENT(INOUT) :: dfmax, uel,off,dam
99 my_real ,
DIMENSION(NEL),
INTENT(OUT) :: tdel
100 my_real,
DIMENSION(NEL,NUVAR),
INTENT(INOUT) :: uvar
101
102
103
104
105
106 INTEGER :: I,J,K,N,NINDX,IFAIL_SH
107 INTEGER, DIMENSION(NEL) :: INDX , CONDITION1, CONDITION2,CONDITION3
108
109 my_real afrac,bfrac,cfrac,nfrac,inst,p_thickfail
111 my_real lode,dam1,dam2,f1,f2,f3,epsf
113
114 g1 = zero
115 g2 = zero
116 epsf = zero
117 dsse = zero
118
119 afrac = uparam(1)
120 bfrac = uparam(2)
121 cfrac = uparam(3)
122 nfrac = uparam(4)
123 inst = uparam(5)
124 ifail_sh = uparam(7)
125 condition1(1:nel) = 0
126 condition2(1:nel) = 0
127 condition3(1:nel) = 0
128 nindx = 0
129
130
131
132 DO i =1,nel
133
134 dam1 = dfmax(i)
135 dam2 = uvar(i,2)
136 IF (fld_idx(i) == 0) fld_idx(i) = 1
137
138 IF (dpla(i) > zero .and. off(i) == one .and. dam1 < one) THEN
139 sigm = (signxx(i)+ signyy(i))*third
140 sigvm = sqrt((signxx(i)**2)+(signyy(i)**2)-(signxx(i)*signyy(i))+(3*signxy(i)**2))
141 eta = sigm /
max(em10,sigvm)
142
143 xi = -27.0*0.5*eta*(eta**2 - third)
144 IF (xi < -one) xi =-one
145 IF (xi > one) xi = one
146
147 lode = one - two*acos(xi)/pi
148
149 f1 = two/three*cos((one -lode)*pi/six)
150 f2 = two/three*cos((three+lode)*pi/six)
151 f3 =-two/three*cos((one +lode)*pi/six)
152
153 ghc= (half*((f1-f2)**afrac
154 . +(f2-f3)**afrac
155 . +(f1-f3)**afrac))**(one/afrac)
156 . +cfrac*(two*eta+f1+f3)
157
158 IF (eta < -third) THEN
159 epsf = 100.00
160 ghc = -one
161 ELSEIF(eta < two_third)THEN
162 ghc= (half*(
max(f1-f2,zero)**afrac
163 . +
max(f2-f3,zero)**afrac
164 . +
max(f1-f3,zero)**afrac))**(one/afrac)
165 . +cfrac*(two*eta+f1+f3)
166 epsf = bfrac*((one+cfrac)/ghc)**(1.0/nfrac)
167 IF (eta > third .and. dam2 < one) THEN
168 g1 = three/two*eta+sqrt(third+em20-0.75*(eta**2))
169 g2 = three/two*eta-sqrt(third+em20-0.75*(eta**2))
170 g3 = (half*((
max(zero,(g1-g2)))**inst+g1**inst+
172 dsse = bfrac * g3**(-100)
173 dam2 = dam2 + dpla(i)/
max(em6,dsse)
174 uvar(i,2) = dam2
175 ENDIF
176 ELSE
177 eta = two_third
178 f1 = one/three
179 f2 = one/three
180 f3 = -two/three
181
182 ghc= (half*(
max(f1-f2,zero)**afrac
183 . +
max(f2-f3,zero)**afrac
184 . +
max(f1-f3,zero)**afrac))**(one/afrac)
185 . +cfrac*(two*eta+f1+f3)
186
187 epsf = bfrac*(one+cfrac)**(one/nfrac)
188 . * ((half*(
max(f1-f2,zero)**afrac
189 . +
max(f2-f3,zero)**afrac
190 . +
max(f1-f3,zero)**afrac))**(one/afrac)
191 . + cfrac*(two*eta+f1+f3))**(-one/nfrac)
192
193
194 ENDIF
195
196 dam1 = dam1 + dpla(i)/
max(em6,epsf)
197 dfmax(i) =
min(one,dam1)
198 IF ((dam1 >= one).AND.(ifail_sh < 3)) THEN
199 nindx = nindx + 1
200 indx(nindx) = i
201 foff(i) = 0
202 uel(i) = zero
203 tdel(i) = time
204 condition1(nindx) = ipt
205 ENDIF
206
207 IF (dam2 >= one .AND. uvar(i,1) == zero .AND. foff(i) == 1 .AND. ifail_sh < 3)THEN
208 uvar(i,1) = one
209 uel(i) = uel(i) + one
210 IF (int(uel(i)+1e-6) == nptot) THEN
211 nindx = nindx + 1
212 indx(nindx) = i
213 tdel(i)= time
214 off(i) = zero
215 condition2(nindx) = 1
216 condition3(nindx) = 1
217 ENDIF
218 ENDIF
219
220
221 dam(i) =
min(pla(i)/
max(em6,epsf),one)
222
223
224
225 IF (((eta < third) .AND.(dam1<one)).OR.
226 . ((eta >= third).AND.(dam2<one))) THEN
227 fld_idx(i) = 1
228
229 ELSEIF ((eta >= third).AND.(dam2>=one).AND.(dam1<one)) THEN
230 fld_idx(i) = 2
231
232 ELSEIF (dam1 >= one) THEN
233 fld_idx(i) = 3
234 ENDIF
235
236 ENDIF
237 ENDDO
238
239
240 IF (nindx > 0) THEN
241 DO j=1,nindx
242 i = indx(j)
243#include "lockon.inc"
244 IF(condition1(j) >= 1) THEN
245 WRITE(iout, 2000) ngl(i),condition1(j),time
246 WRITE(istdo,2000) ngl(i),condition1(j),time
247 ENDIF
248 IF(condition3(j)==1) THEN
249 WRITE(istdo,3000) ngl(i)
250 WRITE(iout, 3000) ngl(i)
251 ENDIF
252 IF(condition2(j)==1) THEN
253 WRITE(istdo,2200) ngl(i), time
254 WRITE(iout, 2200) ngl(i), time
255 ENDIF
256#include "lockoff.inc"
257 END DO
258 END IF
259
260
261 2000 FORMAT(1x,'FOR SHELL ELEMENT (HC-DSSE)',i10,1x,'LAYER',i3,':',/,
262 . 1x,'STRESS TENSOR SET TO ZERO',1x,'AT TIME :',1pe12.4)
263
264 2200 FORMAT(1x,' *** RUPTURE OF SHELL ELEMENT (HC-DSSE)',i10,1x,
265 . ' AT TIME :',1pe12.4)
266 3000 FORMAT(1x,'FOR SHELL ELEMENT (HC-DSSE)',i10,
267 . 1x,'INSTABILITY REACHED.')
268
269
270
271
272
273
274
275 RETURN