41
42
43
44#include "implicit_f.inc"
45#include "comlock.inc"
46
47
48
49#include "mvsiz_p.inc"
50
51
52
53#include "com01_c.inc"
54#include "com04_c.inc"
55#include "com08_c.inc"
56#include "scr02_c.inc"
57#include "scr07_c.inc"
58#include "scr17_c.inc"
59#include "scr18_c.inc"
60#include "param_c.inc"
61#include "cong1_c.inc"
62#include "units_c.inc"
63#include "impl1_c.inc"
64#include "sms_c.inc"
65
66
67
68 INTEGER, INTENT(IN) :: ITY
69 INTEGER, INTENT(IN) :: JTUR
70 INTEGER, INTENT(IN) :: JTHE
71 INTEGER, INTENT(IN) :: JSMS
72 INTEGER :: NEL
73
75 . pm(npropm,*), off(mvsiz) , rho(nel) , rk(mvsiz) , t(mvsiz),
76 . re(mvsiz) , sti(*) , eint(nel),
77 . d1(mvsiz,8) , d2(mvsiz,8), d3(mvsiz,8),
78 . vol(mvsiz) , dvol(mvsiz),
79 . vd2(mvsiz) ,deltax(mvsiz),vis(mvsiz),qold(nel), ssp(mvsiz),
80 . geo(npropg,*) , dt2t , offg(*), mssa(*), dmels(*)
81 INTEGER MAT(MVSIZ),NC(8,MVSIZ),NGL(MVSIZ), PID(*), NELTST,ITYPTST
82
83
84
85 INTEGER I,J, MT,IPT
86
88 . dd(mvsiz), al(mvsiz), dtx(mvsiz), dty(mvsiz),
89 . ad(mvsiz), qx(mvsiz), cx(mvsiz),qvis(mvsiz),
90 . visi, facq, qa, qb,
91 . cns1, cns2, sph, ak1, bk1, ak2, bk2, tli, akk, xmu, tmu, rpr,
92 . atu
93
94
95 DO i=1, nel
96 al(i)=zero
97 ENDDO
98 DO 2 i=1,nel
99 2 dd(i)=zero
100 DO 5 ipt=1,8
101 DO 5 i=1,nel
102 5 dd(i)=dd(i)- one_over_8*(d1(i,ipt)+d2(i,ipt)+d3(i,ipt))
103
104 IF(impl==zero)THEN
105 DO 10 i=1,nel
106 10 cx(i)=ssp(i)+sqrt(vd2(i))
107 visi=1.
108 facq=1.
109 IF(impl_s>0.AND.idyna==0)THEN
110 visi=zero
111 facq=zero
112 ENDIF
113 IF(imconv<0) THEN
114 DO i=1,nel
115 vol(i)= abs(vol(i))
116 ENDDO
117 ENDIF
118 ELSE
119 DO 15 i=1,nel
120 15 cx(i)=sqrt(vd2(i))
121 visi=zero
122 facq=zero
123 ENDIF
124
125
126 DO 20 i=1,nel
127 ad(i)=zero
128 IF(off(i)<one.OR.offg(i)<=zero) GO TO 20
129 al(i)=vol(i)** third
130 IF(n2d>0) al(i)=sqrt(vol(i))
131 ad(i)=
max(zero,dd(i))*al(i)
132 20 CONTINUE
133
134 IF(invstr>=35)THEN
135 DO i=1,nel
136 qa =facq*geo(14,pid(i))
137 qb =facq*geo(15,pid(i))
138 cns1=geo(16,pid(i))
139 cns2=geo(17,pid(i))*ssp(i)*al(i)*rho(i)
140 qx(i)=(qb+cns1)*ssp(i)+qa*qa*deltax(i) *
max(zero,dd(i))
141 . + visi*(two*vis(i)+cns2)
142 . /
max(em20,rho(i)*deltax(i))
143 qvis(i)=rho(i)*ad(i)*(qa*qa*ad(i)+qb*ssp(i))
144 ENDDO
145 ELSE
146 mt = mat(1)
147 qa =facq*pm(2,mt)
148 qb =facq*pm(3,mt)
149 cns1=pm(93,mt)
150 DO i=1,nel
151 cns2=pm(94,mt)*ssp(i)*al(i)*rho(i)
152
153
154 qx(i)=(qb+cns1)*ssp(i)+qa*qa*deltax(i) *
max(zero,dd(i))
155 . + visi*(two*vis(i)+cns2)
156 . /
max(em20,rho(i)*deltax(i))
157 qvis(i)=rho(i)*ad(i)*(qa*qa*ad(i)+qb*ssp(i))
158 ENDDO
159 ENDIF
160
161 DO 30 i=1,nel
162 30 dtx(i)=deltax(i)/
163 .
max(em20,qx(i)+sqrt(qx(i)*qx(i)+cx(i)*cx(i)))
164
165 IF(jthe > 0 )THEN
166 mt = mat(1)
167 sph = pm(69,mt)
168 ak1 = pm(75,mt)
169 bk1 = pm(76,mt)
170 ak2 = pm(77,mt)
171 bk2 = pm(78,mt)
172 tli = pm(80,mt)
173 DO 40 i=1,nel
174 IF(t(i)<tli)THEN
175 akk=ak1+bk1*t(i)
176 ELSE
177 akk=ak2+bk2*t(i)
178 ENDIF
179 IF(jtur/=0)THEN
180 xmu = rho(i)*pm(24,mt)
181 tmu = pm(81,mt)
182 rpr = pm(95,mt)
183 atu=rpr*tmu*rk(i)**2/(
max(em15,re(i)*vol(i))*xmu)
184 akk=akk*(1.+atu)
185 ENDIF
186 dtx(i) =
min(dtx(i),half*deltax(i)**2*sph/akk)
187 40 CONTINUE
188 ENDIF
189
190
191
192 IF(.NOT.(idtmins==2.AND.jsms/=0))THEN
193 DO i=1,nel
194 eint(i)=eint(i)+half*off(i)*dvol(i)*(-qold(i)-qvis(i))
195 qold(i)=qvis(i)
196
197
198 sti(i) = off(i) * rho(i) * vol(i) / dtx(i)**2
199 dtx(i)= dtfac1(ity)*dtx(i)
200 ENDDO
201
202 IF(nodadt==0)THEN
203 DO i = 1,nel
204 IF(off(i)==zero.OR.offg(i)<zero) cycle
205 IF(dtx(i)<dt2t)THEN
206 dt2t = dtx(i)
207 neltst = ngl(i)
208 ityptst = ity
209 END IF
210 ENDDO
211 END IF
212
213 ELSE
214 DO i=1,nel
215 dty(i)= dtx(i)
216 dtx(i)= dtfac1(ity)*dtx(i)
217 END DO
218 END IF
219
220 IF(imconv==1)THEN
221 IF(idtmin(ity)==1)THEN
222 DO 70 i=1,nel
223 IF(dtx(i)>dtmin1(ity).OR.off(i)==zero.OR.
224 . offg(i)<zero)GO TO 70
225 tstop = tt
226#include "lockon.inc"
227 WRITE(iout,*)
228 . ' **ERROR : TIME STEP LESS OR EQUAL DTMIN FOR SOLID ELEMENTS'
229 WRITE(istdo,*)
230 . ' **ERROR : TIME STEP LESS OR EQUAL DTMIN FOR SOLID ELEMENTS'
231#include "lockoff.inc"
232 70 CONTINUE
233 ELSEIF(idtmin(ity)==2)THEN
234 DO 75 i=1,nel
235 IF(dtx(i)>dtmin1(ity).OR.off(i)==zero.OR.
236 . offg(i)<zero)GO TO 75
237 off(i) = zero
238#include "lockon.inc"
239 WRITE(iout,*)
240 . ' -- DELETE SOLID ELEMENTS',ngl(i)
241 WRITE(istdo,*)
242 . ' -- DELETE SOLID ELEMENTS',ngl(i)
243#include "lockoff.inc"
244 idel7nok = 1
245 75 CONTINUE
246 ELSEIF(idtmin(ity)==5)THEN
247 DO 570 i=1,nel
248 IF(dtx(i)>dtmin1(ity).OR.off(i)==zero.OR.
249 . offg(i)<zero)GO TO 570
250 mstop = 2
251#include "lockon.inc"
252 WRITE(iout,*)
253 . ' **ERROR : TIME STEP LESS OR EQUAL DTMIN FOR SOLID ELEMENTS'
254 WRITE(istdo,*)
255 . ' **ERROR : TIME STEP LESS OR EQUAL DTMIN FOR SOLID ELEMENTS'
256#include "lockoff.inc"
257 570 CONTINUE
258 ENDIF
259 END IF
260
261 IF(idtmins == 2 .AND. jsms /= 0)THEN
262 DO i=1,nel
263 eint(i)=eint(i)+half*off(i)*dvol(i)*(-qold
264 qold(i)=qvis(i)
265
266
267 sti(i) = off(i) * rho(i) * vol(i) / dtx(i)**2
268
269
270
271
272
273 dmels(i)=
max(dmels(i),
274 . two*mssa(i)*((dtmins/(dtfacs*dty(i)))**2 - one))
275 dtx(i)=dtfacs*sqrt(one+dmels(i)/(two*mssa(i)))*dty(i)
276 ENDDO
277
278 DO i = 1,nel
279 IF(off(i)==zero.OR.offg(i)<zero) cycle
280 IF(dtx(i)<dt2t)THEN
281 dt2t = dtx(i)
282 neltst = ngl(i)
283 ityptst = ity
284 END IF
285 ENDDO
286
287 END IF
288
289 RETURN