47
49 USE elbufdef_mod
50
51
52
53#include "implicit_f.inc"
54#include "comlock.inc"
55
56
57
58#include "mvsiz_p.inc"
59
60
61
62#include "units_c.inc"
63#include "scr17_c.inc"
64#include "param_c.inc"
65#include "com08_c.inc"
66#include "impl1_c.inc"
67
68
69
70 INTEGER JFT, JLT, MTN, IOFC, IPLA,NEL,IR,IS,ISMSTR,KFTS
71 INTEGER MAT(*), PID(*), NPF(*),NGL(*), INDX(*),IPM(NPROPMI,*)
73 my_real for(nel,5), mom(nel,3), thk(*), eint(jlt,2),pm(npropm,*),
74 . off(*), gstr(nel,8), pla(*), dir(*),viscmx(*),
75 .
area(*),tf(*),dt1c(*),
76 . exx(*), eyy(*), exy(*), exz(*), eyz(*), epsp(*),
77 . kxx(*), kyy(*), kxy(*),bufmat(*),ssp(*),rho(*),
78 . zcfac(mvsiz,2),gs(*),sigy(*),thk0(*),shf(*),f_def(mvsiz,8),
79 . a1(mvsiz),a2(mvsiz),g(mvsiz),nu(mvsiz),vol0(*),
80 . epinchzz(mvsiz),epinchxz(mvsiz),epinchyz(mvsiz),
81 . forp(nel),momp(nel,2),aldt(mvsiz),ezzavg(mvsiz),areapinch(mvsiz)
82 TYPE(TTABLE) TABLE(*)
83 TYPE(ELBUF_STRUCT_), TARGET :: ELBUF_STR
84
85
86
87 INTEGER IGTYP, I, NUVAR ,NINDX,IGEO(NPROPGI,*),
88 . MX,IOFF_DUCT(MVSIZ),ISRATE
89
91 . degmb(mvsiz) ,degfx(mvsiz) ,thkn(mvsiz),
92 . depsxx(mvsiz),depsyy(mvsiz),depszz(mvsiz),
93 . depsxy(mvsiz),depsyz(mvsiz),depszx(mvsiz),
94 . sigoxx(mvsiz),sigoyy(mvsiz),sigozz(mvsiz),
95 . sigoxy(mvsiz),sigoyz(mvsiz),sigozx(mvsiz),
96 . signxx(mvsiz),signyy(mvsiz),signzz(mvsiz),
97 . signxy(mvsiz),signyz(mvsiz),signzx(mvsiz),
98 . sigvxx(mvsiz),sigvyy(mvsiz),sigvxy(mvsiz),
99 . eps_m2,eps_k2,young, visc, vol2, asrate
101 . depbxx(mvsiz),depbyy(mvsiz),depbxy(mvsiz),
102 . deppxz(mvsiz),deppyz(mvsiz),
103 . momoxx(mvsiz),momoyy(mvsiz),momoxy(mvsiz),
104 . momopxz(mvsiz),momopyz(mvsiz),
105 . momnxx(mvsiz),momnyy(mvsiz),momnxy(mvsiz),
106 . momnpxz(mvsiz),momnpyz(mvsiz),
107 . etse(mvsiz) ,epspl(mvsiz),epsp_loc(mvsiz)
109 . DIMENSION(:) ,POINTER :: uvar
110 my_real,
DIMENSION(MVSIZ) :: dt_inv
111
112 TYPE(L_BUFEL_) ,POINTER :: LBUF
113
114
115
116 lbuf => elbuf_str%BUFLY(1)%LBUF(ir,is,1)
117
118
119 igtyp = igeo(11,pid(1))
120 nuvar = elbuf_str%BUFLY(1)%NVAR_MAT
121 uvar =>elbuf_str%BUFLY(1)%MAT(ir,is,1)%VAR
122 ioff_duct(1:mvsiz) = 0
123 viscmx(1:mvsiz) = zero
124
125 DO i=jft,jlt
126 degmb(i) =
for(i,1)*exx(i)+
for(i,2)*eyy(i)+
for(i,3)*exy(i)
127 . +
for(i,4)*eyz(i)+
for(i,5)*exz(i)
128 . + half*forp(i)*epinchzz(i)
129 degfx(i) = mom(i,1)*kxx(i)+mom(i,2)*kyy(i)+mom(i,3)*kxy(i)
130 . + half*momp(i,1)*epinchxz(i)+half*momp(i,2)*epinchyz(i)
131 ENDDO
132
133 DO i=jft,jlt
134 thkn(i) = thk(i)
135 ENDDO
136! compute
the inverse of dt and
save the result
137 DO i=jft,jlt
138 dt_inv(i) = dt1c(i)/
max(dt1c(i)**2,em20)
139 ENDDO
140
141
142
143
144
145
146 mx = mat(jft)
147
148 israte = ipm(3,mx)
149 IF (israte == 1) THEN
150!#include "vectorize.inc"
151 DO i=jft,jlt
152 eps_k2 = (kxx(i)*kxx(i)+kyy(i)*kyy(i)+kxx(i)*kyy(i)
153 . + fourth*kxy(i)*kxy(i)) *thk(i)*thk(i)*one_over_9
154 eps_m2 = four_over_3*(exx(i)*exx(i)+eyy(i)*eyy(i)+exx(i)*eyy(i)
155 . + fourth*exy(i)*exy(i))
156 epsp_loc(i) = sqrt(eps_k2 + eps_m2)*dt_inv(i)
157 END DO
158 ELSE IF (israte > 1) THEN
159
160 DO i=jft,jlt
161 eps_m2 = four_over_3*(exx(i)*exx(i)+eyy(i)*eyy(i)+exx(i)*eyy(i)
162 . + fourth*exy(i)*exy(i))
163 epsp_loc(i) = sqrt(eps_m2)*dt_inv(i)
164 END DO
165 ENDIF
166
167 IF (israte > 0) THEN
168 DO i=jft,jlt
169 asrate =
min(one,pm(9,mx)*dt1c(i))
170 epsp(i) = asrate*epsp_loc(i) + (one-asrate)*epsp(i)
171 epspl(i) = epsp(i)
172 END DO
173 ENDIF
174
175
176
177 DO i=jft,jlt
178 signxx(i) = zero
179 signyy(i) = zero
180 signzz(i) = zero
181 signxy(i) = zero
182 signyz(i) = zero
183 signzx(i) = zero
184 ENDDO
185 DO i=jft,jlt
186 momnxx(i) = zero
187 momnyy(i) = zero
188 momnxy(i) = zero
189 momnpxz(i)=zero
190 momnpyz(i)=zero
191 ENDDO
192
193 DO i=jft,jlt
194 depsxx(i) = exx(i)
195 depsyy(i) = eyy(i)
196 depszz(i) = epinchzz(i)
197 depsxy(i) = exy(i)
198 depsyz(i) = eyz(i)
199 depszx(i) = exz(i)
202 sigozz(i) = forp(i)
206 ENDDO
207
208 DO i=jft,jlt
209 depbxx(i) = kxx(i)
210 depbyy(i) = kyy(i)
211 depbxy(i) = kxy(i)
212 deppxz(i) = epinchxz(i)
213 deppyz(i) = epinchyz(i)
214 momoxx(i) = mom(i,1)
215 momoyy(i) = mom(i,2)
216 momoxy(i) = mom(i,3)
217 momopxz(i)= momp(i,1)
218 momopyz(i)= momp(i,2)
219 ENDDO
220
221 IF(mtn == 1) THEN
223 1 jft ,jlt ,g ,thkn ,off ,
224 2 gs ,a1 ,a2 ,nu ,thk0 ,
225 3 nel ,ssp ,rho ,
226 4 depsxx ,depsyy ,depszz ,
227 5 depsxy ,depsyz ,depszx ,
228 6 depbxx ,depbyy ,depbxy ,
229 7 deppxz ,deppyz ,
230 8 sigoxx ,sigoyy ,sigozz ,
231 9 sigoxy ,sigoyz ,sigozx ,
232 a momoxx ,momoyy ,momoxy ,
233 b momopxz ,momopyz ,
234 c signxx ,signyy ,signzz ,
235 d signxy ,signyz ,signzx ,
236 e momnxx ,momnyy ,momnxy ,
237 f momnpxz ,momnpyz)
238 ELSEIF(mtn == 91) THEN
240 1 jft ,jlt ,nuvar ,bufmat ,rho ,
241 2 thkn ,thk0 ,nel ,ssp ,
area ,
242 3 depsxx ,depsyy ,depszz ,
243 4 depsxy ,depsyz ,depszx ,
244 5 depbxx ,depbyy ,depbxy ,
245 6 deppxz ,deppyz ,
246 7 sigoxx ,sigoyy ,sigozz ,
247 8 sigoxy ,sigoyz ,sigozx ,
248 9 momoxx ,momoyy ,momoxy ,
249 a momopxz ,momopyz ,
250 b signxx ,signyy ,signzz ,
251 c signxy ,signyz ,signzx ,
252 d momnxx ,momnyy ,momnxy ,
253 e momnpxz ,momnpyz ,tt ,uvar ,dt_inv ,
254 f viscmx ,aldt ,vol0 ,ipm ,mat ,
255 g pla ,degmb ,degfx ,
256 h ngl ,ezzavg ,areapinch)
257 ENDIF
258
259
260
261 DO i=jft,jlt
267 mom(i,1) = momnxx(i)
268 mom(i,2) = momnyy(i)
269 mom(i,3) = momnxy(i)
270 forp(i) = signzz(i)
271 momp(i,1)= momnpxz(i)
272 momp(i,2)= momnpyz(i)
273 ENDDO
274
275 DO i=jft,jlt
281 mom(i,1) =mom(i,1)*off(i)
282 mom(i,2) =mom(i,2)*off(i)
283 mom(i,3) =mom(i,3)*off(i)
284 forp(i) =forp(i)*off(i)
285 momp(i,1) =momp(i,1)*off(i)
286 momp(i,2) =momp(i,2)*off(i)
287 ENDDO
288
289 DO i=jft,jlt
290 degmb(i) = degmb(i)+
for(i,1)*exx(i)+
for(i,2)*eyy(i)+
for(i,3)*exy(i)
291 . +
for(i,4)*eyz(i)+
for(i,5)*exz(i)
292 . + half*forp(i)*epinchzz(i)
293 degfx(i) = degfx(i)+mom(i,1)*kxx(i)+mom(i,2)*kyy(i)+mom(i,3)*kxy(i)
294 . + half*momp(i,1)*epinchxz(i)+half*momp(i,2)*epinchyz(i)
295 vol2 = half*thk0(i)*
area(i)*off(i)
296 eint(i,1) = eint(i,1) + degmb(i)*vol2
297 eint(i,2) = eint(i,2) + degfx(i)*thk0(i)*vol2
298 ENDDO
299
300
301
302
303 nindx=0
304 DO i=jft,jlt
305 IF (off(i) == four_over_5 . and. ioff_duct(i) == 0) THEN
306 off(i)= zero
307 nindx=nindx+1
308 indx(nindx)=i
309 ENDIF
310 ENDDO
311 IF (nindx > 0) THEN
312 idel7nok = 1
313 IF (imconv == 1) THEN
314 DO i = 1, nindx
315#include "lockon.inc"
316 WRITE(iout, 1000) ngl(indx(i))
317 WRITE(istdo,1100) ngl(indx(i)),tt
318#include "lockoff.inc"
319 ENDDO
320 ENDIF
321 ENDIF
322 iofc = nindx
323
324 1000 FORMAT(1x,'-- RUPTURE OF SHELL ELEMENT NUMBER ',i10)
325 1100 FORMAT(1x,'-- RUPTURE OF SHELL ELEMENT :',i10,' AT TIME :',g11.4)
326
327 RETURN
end diagonal values have been computed in the(sparse) matrix id.SOL
subroutine area(d1, x, x2, y, y2, eint, stif0)
for(i8=*sizetab-1;i8 >=0;i8--)
subroutine sigeps01gpinch(jft, jlt, g, thk, off, gs, a1, a2, nu, thk0, nel, ssp, rho, depsxx, depsyy, depszz, depsxy, depsyz, depszx, depbxx, depbyy, depbxy, deppxz, deppyz, sigoxx, sigoyy, sigozz, sigoxy, sigoyz, sigozx, momoxx, momoyy, momoxy, momopxz, momopyz, signxx, signyy, signzz, signxy, signyz, signzx, momnxx, momnyy, momnxy, momnpxz, momnpyz)
subroutine sigeps91gpinch(jft, jlt, nuvar, uparam, rho0, thk, thk0, nel, ssp, area, depsxx, depsyy, depszz, depsxy, depsyz, depszx, depbxx, depbyy, depbxy, deppxz, deppyz, sigoxx, sigoyy, sigozz, sigoxy, sigoyz, sigozx, momoxx, momoyy, momoxy, momopxz, momopyz, signxx, signyy, signzz, signxy, signyz, signzx, momnxx, momnyy, momnxy, momnpxz, momnpyz, time, uvar, dt_inv, viscmx, aldt, vol0, ipm, mat, pla, degmb, degfx, ngl, ezzavg, areapinch)