41
42
43
44 USE elbufdef_mod
45
46
47
48#include "implicit_f.inc"
49#include "comlock.inc"
50
51
52
53#include "mvsiz_p.inc"
54
55
56
57#include "units_c.inc"
58#include "param_c.inc"
59#include "scr17_c.inc"
60#include "com08_c.inc"
61#include "impl1_c.inc"
62
63
64
65 INTEGER JFT,,IOFC,NEL,IOFF_DUCT(*),MX
66 INTEGER NGL(MVSIZ),INDX(MVSIZ)
68 . pm(npropm,*),
69 . off(*),dt1c(*),gs(*),thk(*),epsp(mvsiz)
71 . exz(mvsiz),eyz(mvsiz),kxx(mvsiz),kyy(mvsiz),kxy(mvsiz),
72 . nu(mvsiz),thk0(mvsiz),degmb(mvsiz),degfx(mvsiz)
74 .
for(nel,5),mom(nel,3),
75 . depsxx(nel),depsyy(nel),depsxy(nel),depsyz(nel),depszx(nel),
76 . depbxx(nel),depbyy(nel),depbxy(nel),
77 . sigoxx(nel),sigoyy(nel),sigoxy(nel),sigoyz(nel),sigozx(nel),
78 . signxx(nel),signyy(nel),signxy(nel),signyz(nel),signzx(nel),
79 . momoxx(nel),momoyy(nel),momoxy(nel),
80 . momnxx(nel),momnyy(nel),momnxy(nel)
81 TYPE(ELBUF_STRUCT_), TARGET :: ELBUF_STR
82
83
84
85 INTEGER ICC(MVSIZ),I,J,ICC_1,NINDX
87 . f1(mvsiz),f2(mvsiz),f3(mvsiz),f4(mvsiz),f5(mvsiz),
88 . ym(mvsiz),fp1(mvsiz),fp2(mvsiz),fp3(mvsiz),
89 . m1(mvsiz),m2(mvsiz),m3(mvsiz),mp1(mvsiz),mp2(mvsiz),mp3(mvsiz),
90 . dwelm(mvsiz),dwelf(mvsiz),ca(mvsiz),cb(mvsiz),cn(mvsiz),
91 . epmx(mvsiz),
ymax(mvsiz),yeq(mvsiz),dwpla(mvsiz),hh(mvsiz),
92 . rr(mvsiz),c1(mvsiz),c2(mvsiz),c3(mvsiz),cc(mvsiz),
93 . b1(mvsiz),b2(mvsiz),b3(mvsiz),epdr(mvsiz),epsl(mvsiz),hl(mvsiz),
94 . yldl(mvsiz),a1(mvsiz),a2(mvsiz),g(mvsiz),
95 . degmb_loc(mvsiz),degsh_loc(mvsiz),degfx_loc(mvsiz),
96 . alpe(mvsiz),yld(mvsiz)
98 . depsl,ymi,a1i,thk12,ezz,dpla,aaa,bbb,ccc,
99 . p1,p2,p3,q1,q2,q3,ym_1,a1_1,c1_1,c2_1,c3_1,ca_1,cb_1,
100 . cn_1,epmx_1,ymax_1,cc_1,epsl_1,hl_1,yldl_1
101
102 TYPE(G_BUFEL_) ,POINTER :: GBUF
103
104 gbuf => elbuf_str%GBUF
105
106 DATA p1/ 0.9659258/, p2/-0.2588190/, p3/ 1.7320508/
107 DATA q1/ 1.1153548/, q2/+0.2988584/, q3/ 0.5773503/
108
109 ym_1 = pm(20,mx)
110 a1_1 = pm(24,mx)
111 c1_1 = pm(28,mx)
112 c2_1 = pm(29,mx)
113 c3_1 = pm(30,mx)
114 ca_1 = pm(38,mx)
115 cb_1 = pm(39,mx)
116 cn_1 = pm(40,mx)
117 epmx_1 = pm(41,mx)
118 ymax_1 = pm(42,mx)
119 cc_1 = pm(43,mx)
120 epsl_1 = pm(45,mx)
121 hl_1 = pm(46,mx)
122 yldl_1 = pm(47,mx)
123 icc_1 = nint(pm(49,mx))
124
125 DO i=jft,jlt
126 ym(i) = ym_1
127 a1(i) = a1_1
128 c1(i) = c1_1
129 c2(i) = c2_1
130 c3(i) = c3_1
131 ca(i) = ca_1
132 cb(i) = cb_1
133 cn(i) = cn_1
134 epmx(i) = epmx_1
136 cc(i) = cc_1
137 epdr(i) =
max(em20,pm(44,mx)*dt1c(i))
138 epsl(i) = epsl_1
139 hl(i) = hl_1
140 yldl(i) = yldl_1
141 icc(i) = icc_1
142 ENDDO
143
144
145
146
147 DO i=jft,jlt
148 yld(i) = ca(i)+cb(i)*exp(cn(i) * log(gbuf%PLA(i)+em30))
150 depsl =
max(zero,gbuf%PLA(i)-epsl(i))
151 yld(i) =
min(yld(i),yldl(i)+hl(i)*depsl)
152 yld(i) =
max(yld(i),zero)
153 alpe(i)=
min(one,yld(i)/(yld(i)+ym(i)*depsl))
154 ymi = alpe(i)*ym(i)
155 g(i) = half*ymi/(one+nu(i))
156 a1i = ymi/(one-nu(i)**2)
157 alpe(i)=
max(em30,a1i/a1(i))
158 a1(i) = a1i
159 a2(i) = nu(i)*a1(i)
160 thk12 = thk0(i)/twelve
161 b1(i) = a1(i)*thk12
162 b2(i) = a2(i)*thk12
163 b3(i) = g(i) *thk12
164 ENDDO
165
166#include "vectorize.inc"
167 DO i=jft,jlt
168 degsh_loc(i) =
for(i,4)*eyz(i)+
for(i,5)*exz(i)
169 degmb_loc(i) = degmb(i) - degsh_loc(i)
170 degfx_loc(i) = degfx(i)
171 ENDDO
172
173
174
175#include "vectorize.inc"
176 DO i=jft,jlt
177 f1(i) = sigoxx(i) + a1(i)*depsxx(i)+a2(i)*depsyy(i)
178 f2(i) = sigoyy(i) + a1(i)*depsyy(i)+a2(i)*depsxx(i)
179 f3(i) = sigoxy(i) + g(i) *depsxy(i)
180 f4(i) = sigoyz(i) + alpe(i)*gs(i)*depsyz(i)
181 f5(i) = sigozx(i) + alpe(i)*gs(i)*depszx(i)
182
183 fp1(i) = p1*f1(i) + p2*f2(i)
184 fp2(i) = p2*f1(i) + p1*f2(i)
185 fp3(i) = p3*f3(i)
186 yeq(i) = fp1(i)**2 + fp2(i)**2
187 yeq(i) = yeq(i) + fp3(i)**2
188
189 m1(i) = momoxx(i) + b1(i)*depbxx(i)+b2(i)*depbyy(i)
190 m2(i) = momoyy(i) + b1(i)*depbyy(i)+b2(i)*depbxx(i)
191 m3(i) = momoxy(i) + b3(i)*depbxy(i)
192
193 mp1(i) = (p1*m1(i) + p2*m2(i))*four
194 mp2(i) = (p2*m1(i) + p1*m2(i))*four
195 mp3(i) = p3*m3(i) * four
196 yeq(i) = yeq(i) + mp1(i)**2+mp2(i)**2
197 yeq(i) = sqrt(yeq(i) + mp3(i)**2)
198 ENDDO
199
200
201
202 DO i=jft,jlt
203 epsp(i) = abs(degmb_loc(i)+degfx_loc(i)*thk0(i))/(yeq(i)+em20)
204 epsp(i) =
max(epsp(i),epdr(i))
205 yld(i) = yld(i)*(one+cc(i)*log(epsp(i)/epdr(i)))
206 IF (icc(i) == 2) yld(i)=
min(yld(i),
ymax(i))
207 rr(i) =
min(one,yld(i)/(yeq(i)+em20))
208 ENDDO
209
210
211
212#include "vectorize.inc"
213 DO i=jft,jlt
214 f1(i) =(q1*fp1(i) + q2*fp2(i))*rr(i)
215 f2(i) =(q2*fp1(i) + q1*fp2(i))*rr(i)
216 f3(i) = q3*fp3(i)*rr(i)
217
218 dwelm(i) =(f1(i)+sigoxx(i))*(c1(i)*(f1(i)-sigoxx(i))+
219 . c2(i)*(f2(i)-sigoyy(i)))+
220 . (f2(i)+sigoyy(i))*(c2(i)*(f1(i)-sigoxx(i))+
221 . c1(i)*(f2(i)-sigoyy(i)))+
222 . (f3(i)+sigoxy(i))*(c3(i)*(f3(i)-sigoxy(i)))
223 degmb_loc(i) = degmb_loc(i) + f1(i)*depsxx(i)+f2(i)*depsyy(i)+f3(i)*depsxy(i)
224 ENDDO
225
226 DO i=jft,jlt
227 m1(i) =(q1*mp1(i) + q2*mp2(i))*rr(i)*fourth
228 m2(i) =(q2*mp1(i) + q1*mp2(i))*rr(i)*fourth
229 m3(i) = q3*mp3(i)*rr(i)*fourth
230
231 dwelf(i) =(m1(i)+momoxx(i))*twelve*(c1(i)*(m1(i)-momoxx(i))+
232 . c2(i)*(m2(i)-momoyy(i)))+
233 . (m2(i)+momoyy(i))*twelve*(c2(i)*(m1(i)-momoxx(i))+
234 . c1(i)*(m2(i)-momoyy(i)))+
235 . (m3(i)+momoxy(i))*twelve*(c3(i)*(m3(i)-momoxy(i)))
236 degfx_loc(i) = degfx_loc(i) + m1(i)*depbxx(i)+m2(i)*depbyy(i)+m3(i)*depbxy(i)
237 ENDDO
238
239#include "vectorize.inc"
240 DO i=jft,jlt
241 dwpla(i) = degmb_loc(i) + degfx_loc(i)*thk0(i)-dwelm(i)-dwelf(i)
242 ENDDO
243
244
245
246 DO i=jft,jlt
247 dpla = off(i) *
max(zero,half*dwpla(i)/
max(em20,yld(i)))
248 gbuf%PLA(i) = gbuf%PLA(i) + dpla
249 aaa = abs(dwelm(i)+dwelf(i))
250 bbb = abs(dwpla(i))
251 ccc =
max(em20,aaa+bbb)
252 ezz = - (depsxx(i) + depsyy(i)) * (nu(i)*aaa/(one-nu(i)) + bbb) / ccc
253 thk(i) = thk(i) * (one + ezz*off(i))
254 ENDDO
255!--------------------------------
256
257
258 DO i=jft,jlt
259 IF (off(i) < em01) off(i) = zero
260 IF (off(i) < one) off(i) = off(i)*four_over_5
261 ENDDO
262
263 nindx=0
264
265 DO i=jft,jlt
266 IF (off(i) < one) cycle
267 IF (gbuf%PLA(i) < epmx(i)) cycle
268 off(i) = four_over_5
269 nindx = nindx + 1
270 indx(nindx) = i
271 ioff_duct(i) = 1
272 ENDDO
273
274 IF (nindx > 0) THEN
275 IF (inconv == 1) THEN
276 DO j=1,nindx
277#include "lockon.inc"
278 WRITE(iout, 1000) ngl(indx(j))
279 WRITE(istdo,1100) ngl(indx(j)),tt
280#include "lockoff.inc"
281 ENDDO
282 ENDIF
283 ENDIF
284
285 iofc = nindx
286
287 DO i=jft,jlt
288 signxx(i) = f1(i)
289 signyy(i) = f2(i)
290 signxy(i) = f3(i)
291 signyz(i) = f4(i)
292 signzx(i) = f5(i)
293 momnxx(i) = m1(i)
294 momnyy(i) = m2(i)
295 momnxy(i) = m3(i)
296 ENDDO
297
298 1000 FORMAT(1x,'-- RUPTURE OF SHELL ELEMENT NUMBER ',i10)
299 1100 FORMAT(1x,'-- RUPTURE OF SHELL ELEMENT :',i10,' AT TIME :',g11.4)
300
301 RETURN
subroutine ymax(idn, fac, npc, pld, stiffmin, stiffmax, stiffini, stiffavg)
for(i8=*sizetab-1;i8 >=0;i8--)
1.15.0