37
38
39
40 use matparam_def_mod
41
42
43
44#include "implicit_f.inc"
45#include "comlock.inc"
46
47
48
49#include "units_c.inc"
50#include "param_c.inc"
51#include "scr17_c.inc"
52#include "com08_c.inc"
53#include "impl1_c.inc"
54
55
56
57 INTEGER ,INTENT(IN) :: NEL,NUVAR
58 INTEGER ,INTENT(IN) :: PID(NEL),(NEL)
59 my_real :: off(*), pla(*),al(nel),
60 .
for(nel,3), mom(nel,3), eint(nel,2), geo(npropg,*),
61 . exx(nel),exy(nel),exz(nel),kxx(nel),kyy(nel),
62 . kzz(nel),fa1(nel),fa2(nel),fa3(nel),
63 . ma1(nel),ma2(nel),ma3(nel),a1(nel)
64 my_real ,
DIMENSION(NEL,NUVAR),
INTENT(INOUT) :: uvar
65 my_real ,
DIMENSION(NEL) ,
INTENT(INOUT) :: sigy
66 type (matparam_struct_) ,intent(in) :: mat_param
67
68
69
70 INTEGER :: I,J,NINDX,IFORM,ICC,VP
71 INTEGER :: INDX(NEL)
73 my_real :: degfx(nel),esp(nel),g(nel),
74 . yldtmp(nel),ym(nel),b1(nel), b2(nel), b3(nel), degmb(nel),
75 . dmpm(nel), dmpf(nel), shf(nel),
76 . f1(nel), m1(nel), m2(nel), m3(nel),
77 . degsh(nel), yeq(nel), yld(nel), dwpla(nel),
78 . epmx(nel), dwelm(nel), dwelf(nel),z3(nel),z4(nel),
79 . ca(nel), cb(nel), cn(nel),
ymax(nel),
80 . rr(nel), epsp(nel),
81 . sh(nel), yma2(nel), sh10(nel), sh20(nel), sh0(nel),
82 . sh1(nel), sh2(nel),plap(nel),dpla(nel),
83 . tmp1(nel), tmp2(nel), tmp3(nel)
84
85 IF (impl_s == 0 .OR. idyna > 0) THEN
86 DO i=1,nel
87 dmpm(i)=geo(16,pid(i))*al(i)
88 dmpf(i)=geo(17,pid(i))*al(i)
89 ENDDO
90 ELSE
91 DO i=1,nel
92 dmpm(i)=zero
93 dmpf(i)=zero
94 ENDDO
95 ENDIF
96
97 iform = mat_param%iparam(1)
98 icc = mat_param%iparam(2)
99 vp = mat_param%iparam(3)
100 cc = mat_param%uparam(6)
101 epdr = mat_param%uparam(7)
102 if (vp == 1) then
103 epdr =
max(epdr, em20)
104 else
105 epdr =
max(epdr*dt1, em20)
106 endif
107 asrate = mat_param%uparam(9)
108 asrate =
min(one, asrate*dt1)
109 DO i=1,nel
110 g(i) = mat_param%shear
111 ym(i) = mat_param%young
112
113 ca(i) = mat_param%uparam(1)
114 cb(i) = mat_param%uparam(2)
115 cn(i) = mat_param%uparam(3)
116 epmx(i)= mat_param%uparam(4)
117 ymax(i)= mat_param%uparam(5)
118
119 a1(i) =geo(1,pid(i))
120 b1(i) =geo(2,pid(i))
121 b2(i) =geo(18,pid(i))
122 b3(i) =geo(4,pid(i))
123 shf(i) =geo(37,pid(i))
124 z3(i) = mat_param%uparam(10)
125 z4(i) = mat_param%uparam(11)
126 ENDDO
127
128
129
130 DO i=1,nel
131 esp(i) = (eint(i,1)+eint(i,2))/al(i)/a1(i)
132 ENDDO
133
134 DO i=1,nel
135 degmb(i) =
for(i,1)*exx(i)
136 degsh(i) =
for(i,2)*exy(i)+
for(i,3)*exz(i)
137 degfx(i) = mom(i,1)*kxx(i)+mom(i,2)*kyy(i)+mom(i,3)*kzz(i)
138 ENDDO
139
140
141
142 DO i=1,nel
143 sh(i)=five_over_6*g(i)*a1(i)
144 yma2(i)=twelve*ym(i)/al(i)**2
145 sh10(i)=yma2(i)*b1(i)
146 sh20(i)=yma2(i)*b2(i)
147 sh0(i)=(one - shf(i))*sh(i)
148 sh1(i)=sh0(i)*sh10(i)/(sh(i)+sh10(i)) + shf(i)*sh10(i)
149 sh2(i)=sh0(i)*sh20(i)/(sh(i)+sh20(i)) + shf(i)*sh20(i)
150
151 f1(i) =
for(i,1)+ exx(i)*a1(i)*ym(i)
152 for(i,2)=
for(i,2)+ exy(i)*sh2(i)
153 for(i,3)=
for(i,3)+ exz(i)*sh1(i)
154 m1(i) =mom(i,1)+ kxx(i)*g(i) *b3(i)
155 m2(i) =mom(i,2)+ kyy(i)*ym(i)*b1(i)
156 m3(i) =mom(i,3)+ kzz(i)*ym(i)*b2(i)
157 ENDDO
158
159
160
161 DO i=1,nel
162 yeq(i)= f1(i)*f1(i) + three * a1(i) *
163 + ( m1(i)*m1(i) /
max(b3(i),em20)
164 + + m2(i)*m2(i) /
max(b1(i),em20)
165 + + m3(i)*m3(i) /
max(b2(i),em20) )
166 yeq(i)= sqrt(yeq(i))/a1(i)
167 ENDDO
168
169
170
171 IF (cc /= zero) THEN
172 DO i=1,nel
173 IF(vp == 1)THEN
174 plap(i) = uvar(i,1)
175 plap(i) =
max(plap(i),epdr)
176 epsp(i) = log(plap(i)/epdr)
177 ELSE
178 epsp(i)=abs(degmb(i)+degfx(i))/(yeq(i)+ em20)/a1(i)
179 tmp2(i)=abs(degmb(i)+degfx(i))
180 tmp3(i)=epsp(i)
181 epsp(i)=
max(epsp(i),epdr)
182 epsp(i)= log(epsp(i)/epdr)
183 ENDIF
184 ENDDO
185 IF (iform == 1) THEN
186 DO i=1,nel
187 epsp(i)=(one + cc * epsp(i))
188 tmp1(i)=epsp(i)
189 epsp(i)= cc*exp((-z3(i)+z4(i) * epsp(i))*esp(i))
190 IF(icc==1)
ymax(i)=
ymax(i) + epsp(i)
191 ca(i) = ca(i) + epsp(i)
192 epsp(i)=one
193 ENDDO
194 ELSE
195 DO i=1,nel
196 epsp(i)=(one + cc * epsp(i))
197 IF (icc == 1)
ymax(i) =
ymax(i) * epsp(i)
198 ENDDO
199 ENDIF
200 ELSE
201 DO i=1,nel
202 epsp(i)=one
203 ENDDO
204 ENDIF
205
206
207
208 DO i=1,nel
209 yld(i)= ca(i) + cb(i) * exp(cn(i) * log(pla(i)+em30))
210 yldtmp(i)=yld(i)
211 yld(i)=
min(yld(i)*epsp(i),
ymax(i))
212 rr(i) =
min(one,yld(i)/(yeq(i)+ em20))
213 sigy(i) = yld(i)
214 ENDDO
215
216 DO i=1,nel
217 f1(i) = f1(i)*rr(i)
218 dwelm(i) =(f1(i)+
for(i,1))*(f1(i)-
for(i,1))/ym(i)/a1(i)
219 degmb(i) = degmb(i) + f1(i)*exx(i)
220 ENDDO
221
222 DO i=1,nel
223 m1(i) = m1(i)*rr(i)
224 m2(i) = m2(i)*rr(i)
225 m3(i) = m3(i)*rr(i)
226 dwelf(i) =(m1(i)+mom(i,1))*(m1(i)-mom(i,1))/ g(i)/b3(i)+
227 . (m2(i)+mom(i,2))*(m2(i)-mom(i,2))/ym(i)/b1(i)+
228 . (m3(i)+mom(i,3))*(m3(i)-mom(i,3))/ym(i)/b2(i)
229 degfx(i) = degfx(i)+ m1(i)*kxx(i)+m2(i)*kyy(i)+m3(i)*kzz(i)
230 ENDDO
231
232 DO i=1,nel
233 dwpla(i)= degmb(i)+degfx(i)-dwelm(i)-dwelf(i)
234 ENDDO
235
236
237
238 DO i=1,nel
239 tmp1(i)=epsp(i)*dwpla(i)/yld(i)
240
241 dpla(i) =
max(zero,half*tmp1(i)/a1(i))
242 pla(i) = pla(i)+off(i) * dpla(i)
243 ENDDO
244 DO i=1,nel
245 IF (vp == 1) THEN
246 plap1 = dpla(i) /
max(em20,dt1)
247 uvar(i,1) = asrate * plap1 + (one - asrate) * plap(i)
248 ENDIF
249 ENDDO
250
251
252
253
254 DO i=1,nel
255 IF (off(i) < em01) off(i) = zero
256 IF (off(i) < one ) off(i) = off(i)*four_over_5
257 ENDDO
258
259 nindx = 0
260 DO i=1,nel
261 IF (off(i) < one) cycle
262 IF (pla(i) < epmx(i)) cycle
263 off(i)=four_over_5
264 idel7nok = 1
265 nindx=nindx+1
266 indx(nindx)=i
267 ENDDO
268
269 IF (nindx > 0 .AND. imconv == 1) THEN
270 DO j=1,nindx
271#include "lockon.inc"
272 WRITE(iout, 1000) ngl(indx(j))
273 WRITE(istdo,1100) ngl(indx(j)),tt
274#include "lockoff.inc"
275 ENDDO
276 ENDIF
277
278 DO i=1,nel
279 fa1(i) = f1(i)
282 ma1(i) = m1(i)
283 ma2(i) = m2(i)
284 ma3(i) = m3(i)
285 ENDDO
286
287 1000 FORMAT(1x,'-- RUPTURE OF BEAM ELEMENT NUMBER ',i10)
288 1100 FORMAT(1x,'-- RUPTURE OF BEAM ELEMENT :',i10,' AT TIME :',g11.4)
289
290 RETURN
subroutine ymax(idn, fac, npc, pld, stiffmin, stiffmax, stiffini, stiffavg)
for(i8=*sizetab-1;i8 >=0;i8--)
1.15.0