OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
mat104_nodam_newton.F
Go to the documentation of this file.
1Copyright> OpenRadioss
2Copyright> Copyright (C) 1986-2025 Altair Engineering Inc.
3Copyright>
4Copyright> This program is free software: you can redistribute it and/or modify
5Copyright> it under the terms of the GNU Affero General Public License as published by
6Copyright> the Free Software Foundation, either version 3 of the License, or
7Copyright> (at your option) any later version.
8Copyright>
9Copyright> This program is distributed in the hope that it will be useful,
10Copyright> but WITHOUT ANY WARRANTY; without even the implied warranty of
11Copyright> MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12Copyright> GNU Affero General Public License for more details.
13Copyright>
14Copyright> You should have received a copy of the GNU Affero General Public License
15Copyright> along with this program. If not, see <https://www.gnu.org/licenses/>.
16Copyright>
17Copyright>
18Copyright> Commercial Alternative: Altair Radioss Software
19Copyright>
20Copyright> As an alternative to this open-source version, Altair also offers Altair Radioss
21Copyright> software under a commercial license. Contact Altair to discuss further if the
22Copyright> commercial version may interest you: https://www.altair.com/radioss/.
23!||====================================================================
24!|| mat104_nodam_newton ../engine/source/materials/mat/mat104/mat104_nodam_newton.F
25!||--- called by ------------------------------------------------------
26!|| sigeps104 ../engine/source/materials/mat/mat104/sigeps104.F
27!||====================================================================
29 1 NEL ,NGL ,NUPARAM ,NUVAR ,VOLUME ,FHEAT ,
30 2 TIME ,TIMESTEP,UPARAM ,UVAR ,JTHE ,OFF ,
31 3 RHO0 ,RHO ,PLA ,DPLA ,EPSD ,SOUNDSP ,
32 4 DEPSXX ,DEPSYY ,DEPSZZ ,DEPSXY ,DEPSYZ ,DEPSZX ,
33 5 SIGOXX ,SIGOYY ,SIGOZZ ,SIGOXY ,SIGOYZ ,SIGOZX ,
34 6 SIGNXX ,SIGNYY ,SIGNZZ ,SIGNXY ,SIGNYZ ,SIGNZX ,
35 7 SIGY ,ET ,TEMP ,SEQ ,INLOC ,DPLANL ,
36 8 JLAG )
37 !=======================================================================
38 ! Implicit types
39 !=======================================================================
40#include "implicit_f.inc"
41 !=======================================================================
42 ! Common
43 !=======================================================================
44 !=======================================================================
45 ! Dummy arguments
46 !=======================================================================
47 INTEGER NEL,NUPARAM,NUVAR,JTHE,INLOC
48 INTEGER ,INTENT(IN) :: JLAG
49 INTEGER ,DIMENSION(NEL), INTENT(IN) :: NGL
50 my_real
51 . TIME,TIMESTEP
52 my_real,DIMENSION(NUPARAM), INTENT(IN) ::
53 . UPARAM
54 my_real ,DIMENSION(NEL), INTENT(IN) :: VOLUME
55 my_real ,DIMENSION(NEL), INTENT(INOUT) :: fheat
56 my_real ,DIMENSION(NEL), INTENT(IN) :: rho0,rho,dplanl,
57 . depsxx,depsyy,depszz,depsxy,depsyz,depszx,
58 . sigoxx,sigoyy,sigozz,sigoxy,sigoyz,sigozx
59c
60 my_real ,DIMENSION(NEL), INTENT(OUT) ::
61 . soundsp,sigy,et,
62 . signxx,signyy,signzz,signxy,signyz,signzx
63c
64 my_real ,DIMENSION(NEL), INTENT(INOUT) ::
65 . pla,dpla,epsd,off,temp,seq
66 my_real ,DIMENSION(NEL,NUVAR), INTENT(INOUT) ::
67 . uvar
68 !=======================================================================
69 ! Local Variables
70 !=======================================================================
71 INTEGER I,II,ITER,NITER,NINDX,INDEX(NEL)
72c
73 my_real
74 . young,bulk,lam,g,g2,nu,cdr,kdr,hard,yld0,qvoce,bvoce,jcc,
75 . epsp0,mtemp,tini,tref,eta,cp,dpis,dpad,asrate,afiltr
76 my_real
77 . ldav,
78 . dtemp,dpdt,dlam,ddep
79 my_real
80 . dsdrdj2,dsdrdj3,
81 . dj3dsxx,dj3dsyy,dj3dszz,dj3dsxy,dj3dsyz,dj3dszx,
82 .
83 .
84 . normxx,normyy,normzz,normxy,normyz,normzx,
85 . sig_dfdsig,dfdsig2,
86 . dphi_dsig,dphi_dyld,dphi_dfdr,
87 . dpdt_nl,
88 . dyld_dpla,dyld_dtemp,dtemp_dlam,normsig
89c
90 my_real, DIMENSION(NEL) ::
91 . trsig,
92 . sxx,syy,szz,sxy,syz,szx,sigm,j2,j3,sigdr,yld,weitemp,
93 . hardp,fhard,frate,ftherm,fdr,
94 . phi,ft,dpla_dlam,dphi_dlam,
95 . dpxx,dpyy,dpxy,dpyz,dpzx,dpzz,sigdr2,yld2i
96 !=======================================================================
97 ! DRUCKER - VOCE - JOHNSON-COOK MATERIAL LAW
98 !=======================================================================
99c
100 !=======================================================================
101 ! - INITIALISATION OF COMPUTATION ON TIME STEP
102 !=======================================================================
103 ! Recovering model parameters
104 ! Elastic parameters
105 young = uparam(1) ! Young modulus
106 bulk = uparam(2) ! Bulk modulus
107 g = uparam(3) ! Shear modulus
108 g2 = uparam(4) ! 2*Shear modulus
109 lam = uparam(5) ! Lambda Hooke parameter
110 nu = uparam(6) ! Poisson ration
111 ! Plastic criterion and hardening parameters [Drucker, 1948]
112 cdr = uparam(12) ! Drucker coefficient
113 kdr = uparam(13) ! Drucker 1/K coefficient
114 tini = uparam(14) ! Initial temperature
115 hard = uparam(15) ! Linear hardening
116 yld0 = uparam(16) ! Initial yield stress
117 qvoce = uparam(17) ! 1st Voce parameter
118 bvoce = uparam(18) ! 2nd Voce parameter
119 ! Strain-rate dependence parameters
120 jcc = uparam(20) ! Johnson-Cook strain rate coefficient
121 epsp0 = uparam(21) ! Johnson-Cook reference strain rate
122 ! Thermal softening and self-heating parameters
123 mtemp = uparam(22) ! Thermal softening slope
124 tref = uparam(23) ! Reference temperature
125 eta = uparam(24) ! Taylor-Quinney coefficient
126 cp = uparam(25) ! Thermal mass capacity
127 dpis = uparam(26) ! Isothermal plastic strain rate
128 dpad = uparam(27) ! Adiabatic plastic strain rate
129 ! Plastic strain-rate filtering parameters
130 asrate = uparam(28) ! Plastic strain rate filtering frequency
131 afiltr = asrate*timestep/(asrate*timestep + one)
132c
133 ! Recovering internal variables
134 DO i=1,nel
135 IF (off(i) < em01) off(i) = zero
136 IF (off(i) < one) off(i) = off(i)*four_over_5
137 ! Standard inputs
138 dpla(i) = zero ! Initialization of the plastic strain increment
139 et(i) = one ! Initialization of hourglass coefficient
140 hardp(i) = zero ! Initialization of hardening modulus
141 ENDDO
142c
143 ! Initialization of temperature and self-heating weight factor
144 IF (time == zero .and. jthe == 0) THEN
145 temp(1:nel) = tini
146 ENDIF
147 IF (cp > zero .OR. jthe /= 0) THEN
148 IF (inloc == 0) THEN
149 DO i=1,nel
150 IF (epsd(i) < dpis) THEN
151 weitemp(i) = zero
152 ELSEIF (epsd(i) > dpad) THEN
153 weitemp(i) = one
154 ELSE
155 weitemp(i) = ((epsd(i)-dpis)**2 )
156 . * (three*dpad - two*epsd(i) - dpis)
157 . / ((dpad-dpis)**3)
158 ENDIF
159 ENDDO
160 ENDIF
161 ENDIF
162c
163 ! Computation of the initial yield stress
164 DO i = 1,nel
165 ! a) - Hardening law
166 fhard(i) = yld0 + hard*pla(i) + qvoce*(one-exp(-bvoce*pla(i)))
167 ! b) - Correction factor for strain-rate dependence (Johnson-Cook)
168 frate(i) = one
169 IF (epsd(i) > epsp0) frate(i) = frate(i) + jcc*log(epsd(i)/epsp0)
170 ! c) - Correction factor for thermal softening
171 ftherm(i) = one - mtemp * (temp(i) - tref)
172 ! d) - Computation of the yield function and value check
173 yld(i) = fhard(i)*frate(i)*ftherm(i)
174 ! e) - Checking values
175 yld(i) = max(em10, yld(i))
176 ENDDO
177c
178 !========================================================================
179 ! - COMPUTATION OF TRIAL VALUES
180 !========================================================================
181 DO i=1,nel
182 ! Computation of the trial stress tensor
183 ldav = (depsxx(i) + depsyy(i) + depszz(i)) * lam
184 signxx(i) = sigoxx(i) + depsxx(i)*g2 + ldav
185 signyy(i) = sigoyy(i) + depsyy(i)*g2 + ldav
186 signzz(i) = sigozz(i) + depszz(i)*g2 + ldav
187 signxy(i) = sigoxy(i) + depsxy(i)*g
188 signyz(i) = sigoyz(i) + depsyz(i)*g
189 signzx(i) = sigozx(i) + depszx(i)*g
190 ! Computation of the trace of the trial stress tensor
191 trsig(i) = signxx(i) + signyy(i) + signzz(i)
192 sigm(i) =-trsig(i) * third
193 ! Computation of the deviatoric trial stress tensor
194 sxx(i) = signxx(i) + sigm(i)
195 syy(i) = signyy(i) + sigm(i)
196 szz(i) = signzz(i) + sigm(i)
197 sxy(i) = signxy(i)
198 syz(i) = signyz(i)
199 szx(i) = signzx(i)
200 ! Second deviatoric invariant
201 j2(i) = half*(sxx(i)**2 + syy(i)**2 + szz(i)**2 )
202 . + sxy(i)**2 + syz(i)**2 + szx(i)**2
203 ! Third deviatoric invariant
204 j3(i) = sxx(i) * syy(i) * szz(i)
205 . + sxy(i) * syz(i) * szx(i) * two
206 . - sxx(i) * syz(i)**2
207 . - syy(i) * szx(i)**2
208 . - szz(i) * sxy(i)**2
209 ! Drucker equivalent stress
210 fdr(i) = j2(i)**3 - cdr*(j3(i)**2)
211 ! Checking equivalent stress values
212 IF (fdr(i) > zero) THEN
213 sigdr(i) = kdr * exp((one/six)*log(fdr(i))) ! FDR(I)**(1/6)
214 ELSE
215 sigdr(i) = zero
216 ENDIF
217 ENDDO
218c
219 !========================================================================
220 ! - COMPUTATION OF YIELD FONCTION
221 !========================================================================
222 phi(1:nel) = (sigdr(1:nel) / yld(1:nel))**2 - one
223 ! Checking plastic behavior for all elements
224 nindx = 0
225 DO i=1,nel
226 IF ((phi(i) >= zero).AND.(off(i) == one)) THEN
227 nindx=nindx+1
228 index(nindx)=i
229 ENDIF
230 ENDDO
231c
232 !====================================================================
233 ! - PLASTIC CORRECTION WITH CUTTING PLANE METHOD (SEMI-IMPLICIT)
234 !====================================================================
235c
236 ! Number of iterations
237 niter = 3
238c
239 ! Loop over yielding elements
240#include "vectorize.inc"
241 DO ii=1,nindx
242c
243 ! Number of the element with plastic behaviour
244 i = index(ii)
245c
246 ! Initialization of the iterative Newton procedure
247 sigdr2(i) = sigdr(i)**2
248 yld2i(i) = one / yld(i)**2
249 dpxx(i) = zero
250 dpyy(i) = zero
251 dpzz(i) = zero
252 dpxy(i) = zero
253 dpyz(i) = zero
254 dpzx(i) = zero
255 ENDDO
256c
257 ! Loop over the iterations
258 DO iter = 1, niter
259#include "vectorize.inc"
260 ! Loop over yielding elements
261 DO ii=1,nindx
262 i = index(ii)
263c
264 ! note: in this part, the purpose is to compute for each iteration
265 ! a plastic multiplier allowing to update internal variables to satisfy
266 ! the consistency condition using the cutting plane semi-implicit
267 ! iterative procedure.
268 ! its expression at each iteration is : dlambda = - phi/dphi_dlambda
269 ! -> PHI : current value of yield function (known)
270 ! -> DPHI_DLAMBDA : derivative of PHI with respect to DLAMBDA by taking
271 ! into account of internal variables kinetic :
272 ! plasticity, temperature and damage (to compute)
273c
274 ! 1 - Computation of DPHI_DSIG the normal to the yield surface
275 !-------------------------------------------------------------
276c
277 ! Derivative with respect to the equivalent stress
278 yld2i(i) = one/(yld(i)**2)
279 dphi_dsig = two*sigdr(i)*yld2i(i)
280c
281 ! Computation of the Eulerian norm of the stress tensor
282 normsig = sqrt(signxx(i)*signxx(i)
283 . + signyy(i)*signyy(i)
284 . + signzz(i)*signzz(i)
285 . + two*signxy(i)*signxy(i)
286 . + two*signyz(i)*signyz(i)
287 . + two*signzx(i)*signzx(i))
288 normsig = max(normsig,one)
289c
290 ! Derivative with respect to Fdr
291 fdr(i) = (j2(i)/(normsig**2))**3 - cdr*((j3(i)/(normsig**3))**2)
292 dphi_dfdr = dphi_dsig*kdr*(one/six)*exp(-(five/six)*log(fdr(i)))
293 dsdrdj2 = dphi_dfdr*three*(j2(i)/(normsig**2))**2
294 dsdrdj3 = -dphi_dfdr*two*cdr*(j3(i)/(normsig**3))
295 ! dJ3/dSig
296 dj3dsxx = two_third*(syy(i)*szz(i)-syz(i)**2)/(normsig**2)
297 . - third*(sxx(i)*szz(i)-szx(i)**2)/(normsig**2)
298 . - third*(sxx(i)*syy(i)-sxy(i)**2)/(normsig**2)
299 dj3dsyy = - third*(syy(i)*szz(i)-syz(i)**2)/(normsig**2)
300 . + two_third*(sxx(i)*szz(i)-szx(i)**2)/(normsig**2)
301 . - third*(sxx(i)*syy(i)-sxy(i)**2)/(normsig**2)
302 dj3dszz = - third*(syy(i)*szz(i)-syz(i)**2)/(normsig**2)
303 . - third*(sxx(i)*szz(i)-szx(i)**2)/(normsig**2)
304 . + two_third*(sxx(i)*syy(i)-sxy(i)**2)/(normsig**2)
305 dj3dsxy = two*(sxx(i)*sxy(i) + sxy(i)*syy(i) + szx(i)*syz(i))/(normsig**2)
306 dj3dsyz = two*(sxy(i)*szx(i) + syy(i)*syz(i) + syz(i)*szz(i))/(normsig**2)
307 dj3dszx = two*(sxx(i)*szx(i) + sxy(i)*syz(i) + szx(i)*szz(i))/(normsig**2)
308 ! dPhi/dSig
309 normxx = dsdrdj2*sxx(i)/normsig + dsdrdj3*dj3dsxx
310 normyy = dsdrdj2*syy(i)/normsig + dsdrdj3*dj3dsyy
311 normzz = dsdrdj2*szz(i)/normsig + dsdrdj3*dj3dszz
312 normxy = two*dsdrdj2*sxy(i)/normsig + dsdrdj3*dj3dsxy
313 normyz = two*dsdrdj2*syz(i)/normsig + dsdrdj3*dj3dsyz
314 normzx = two*dsdrdj2*szx(i)/normsig + dsdrdj3*dj3dszx
315c
316 ! 2 - Computation of DPHI_DLAMBDA
317 !---------------------------------------------------------
318c
319 ! a) Derivative with respect stress increments tensor DSIG
320 ! --------------------------------------------------------
321 dfdsig2 = normxx * normxx * g2
322 . + normyy * normyy * g2
323 . + normzz * normzz * g2
324 . + normxy * normxy * g
325 . + normyz * normyz * g
326 . + normzx * normzx * g
327c
328 ! b) Derivatives with respect to plastic strain P
329 ! ------------------------------------------------
330c
331 ! i) Derivative of the yield stress with respect to plastic strain dYLD / dPLA
332 ! ----------------------------------------------------------------------------
333 hardp(i) = hard + qvoce*bvoce*exp(-bvoce*pla(i))
334 dyld_dpla = hardp(i)*frate(i)*ftherm(i)
335c
336 ! ii) Derivative of dPLA with respect to DLAM
337 ! -------------------------------------------
338 sig_dfdsig = signxx(i) * normxx
339 . + signyy(i) * normyy
340 . + signzz(i) * normzz
341 . + signxy(i) * normxy
342 . + signyz(i) * normyz
343 . + signzx(i) * normzx
344 dpla_dlam(i) = sig_dfdsig / yld(i)
345c
346 ! c) Derivatives with respect to the temperature TEMP
347 ! ---------------------------------------------------
348 IF (jthe == 0 .AND. cp > zero .AND. inloc == 0) THEN
349 ! i) Derivative of the yield stress with respect to temperature dYLD / dTEMP
350 ! ---------------------------------------------------------------------------
351 dyld_dtemp = -fhard(i)*frate(i)*mtemp
352 ! ii) Derivative of the temperature TEMP with respect to DLAM
353 ! -----------------------------------------------------------
354 dtemp_dlam = weitemp(i)*eta/(rho0(i)*cp)*sig_dfdsig
355 ELSE
356 dyld_dtemp = zero
357 dtemp_dlam = zero
358 ENDIF
359c
360 ! d) Derivative with respect to the yield stress
361 ! ----------------------------------------------
362 dphi_dyld = -two*sigdr2(i)/(yld(i)**3)
363c
364 ! 3 - Computation of plastic multiplier and variables update
365 !----------------------------------------------------------
366c
367 ! Derivative of PHI with respect to DLAM
368 dphi_dlam(i) = - dfdsig2 + (dphi_dyld*dyld_dpla*dpla_dlam(i))
369 IF (jthe == 0 .AND. cp > zero .AND. inloc == 0) THEN
370 dphi_dlam(i) = dphi_dlam(i) + dphi_dyld*dyld_dtemp*dtemp_dlam
371 ENDIF
372 dphi_dlam(i) = sign(max(abs(dphi_dlam(i)),em20) ,dphi_dlam(i))
373c
374 ! Computation of the plastic multiplier
375 dlam = -phi(i)/dphi_dlam(i)
376c
377 ! Plastic strains tensor update
378 dpxx(i) = dlam * normxx
379 dpyy(i) = dlam * normyy
380 dpzz(i) = dlam * normzz
381 dpxy(i) = dlam * normxy
382 dpyz(i) = dlam * normyz
383 dpzx(i) = dlam * normzx
384c
385 ! Elasto-plastic stresses update
386 signxx(i) = signxx(i) - dpxx(i)*g2
387 signyy(i) = signyy(i) - dpyy(i)*g2
388 signzz(i) = signzz(i) - dpzz(i)*g2
389 signxy(i) = signxy(i) - dpxy(i)*g
390 signyz(i) = signyz(i) - dpyz(i)*g
391 signzx(i) = signzx(i) - dpzx(i)*g
392 trsig(i) = signxx(i) + signyy(i) + signzz(i)
393 sigm(i) =-trsig(i) * third
394 sxx(i) = signxx(i) + sigm(i)
395 syy(i) = signyy(i) + sigm(i)
396 szz(i) = signzz(i) + sigm(i)
397 sxy(i) = signxy(i)
398 syz(i) = signyz(i)
399 szx(i) = signzx(i)
400c
401 ! Cumulated plastic strain and strain rate update
402 ddep = (dlam/yld(i))*sig_dfdsig
403 dpla(i) = max(zero, dpla(i) + ddep)
404 pla(i) = pla(i) + ddep
405c
406 ! Drucker equivalent stress update
407 j2(i) = half*(sxx(i)**2 + syy(i)**2 + szz(i)**2 )
408 . + sxy(i)**2 + syz(i)**2 + szx(i)**2
409 j3(i) = sxx(i) * syy(i) * szz(i) + two * sxy(i) * syz(i) * szx(i)
410 . - sxx(i) * syz(i)**2 - syy(i) * szx(i)**2 - szz(i) * sxy(i)**2
411 fdr(i) = j2(i)**3 - cdr*(j3(i)**2)
412 sigdr(i) = kdr * exp((one/six)*log(fdr(i)))
413c
414 ! Temperature update
415 IF (jthe == 0 .AND. cp > zero .AND. inloc == 0) THEN
416 dtemp = weitemp(i)*yld(i)*ddep*eta/(rho0(i)*cp)
417 temp(i) = temp(i) + dtemp
418 ENDIF
419 ftherm(i) = one - mtemp*(temp(i) - tref)
420c
421 ! Hardening law update
422 fhard(i) = yld0 + hard*pla(i) + qvoce*(one-exp(-bvoce*pla(i)))
423c
424 ! Yield stress update
425 yld(i) = fhard(i)*frate(i)*ftherm(i)
426 yld(i) = max(yld(i), em10)
427c
428 ! Yield function value update
429 sigdr2(i) = sigdr(i)**2
430 yld2i(i) = one / yld(i)**2
431 phi(i) = sigdr2(i) * yld2i(i) - one
432c
433 ENDDO
434 ! End of the loop over the iterations
435 ENDDO
436 ! End of the loop over the iterations
437 !===================================================================
438 ! - END OF PLASTIC CORRECTION WITH CUTTING PLANE ITERATIVE METHOD
439 !===================================================================
440c
441 ! Storing new values
442 DO i=1,nel
443 ! USR Outputs
444 seq(i) = sigdr(i)
445 ! Plastic strain-rate (filtered)
446 dpdt = dpla(i) / max(em20,timestep)
447 epsd(i) = afiltr * dpdt + (one - afiltr) * epsd(i)
448 ! Coefficient for hourglass
449 IF (dpla(i) > zero) THEN
450 et(i) = hardp(i)*frate(i) / (hardp(i)*frate(i) + young)
451 ELSE
452 et(i) = one
453 ENDIF
454 ! Computation of the sound speed
455 soundsp(i) = sqrt((bulk + four_over_3*g)/rho0(i))
456 ! Storing the yield stress
457 sigy(i) = yld(i)
458 ! Non-local temperature
459 IF (inloc > 0 .AND. off(i) == one) THEN
460 IF (cp > zero .AND. jthe == 0) THEN
461 ! Computation of the weight factor
462 dpdt_nl = max(dplanl(i),zero)/max(timestep,em20)
463 IF (dpdt_nl < dpis) THEN
464 weitemp(i) = zero
465 ELSEIF (dpdt_nl > dpad) THEN
466 weitemp(i) = one
467 ELSE
468 weitemp(i) = ((dpdt_nl-dpis)**2 )
469 . * (three*dpad - two*dpdt_nl - dpis)
470 . / ((dpad-dpis)**3)
471 ENDIF
472 dtemp = weitemp(i)*dplanl(i)*yld(i)*eta/(rho0(i)*cp)
473 temp(i) = temp(i) + dtemp
474 ENDIF
475 ENDIF
476 ENDDO
477c
478C-----------------------------------------------
479C plastic work dissipated as heat - in case of /heat/mat
480C-----------------------------------------------
481 IF (jthe < 0 .AND. jlag /= 0) THEN
482 DO i=1,nel
483 fheat(i) = fheat(i) + eta*(one-ft(i))*weitemp(i)*yld(i)*dpla(i)*volume(i)
484 ENDDO
485 END IF
486!-----------
487 RETURN
488 END
end diagonal values have been computed in the(sparse) matrix id.SOL
#define max(a, b)
Definition macros.h:21
subroutine mat104_nodam_newton(nel, ngl, nuparam, nuvar, volume, fheat, time, timestep, uparam, uvar, jthe, off, rho0, rho, pla, dpla, epsd, soundsp, depsxx, depsyy, depszz, depsxy, depsyz, depszx, sigoxx, sigoyy, sigozz, sigoxy, sigoyz, sigozx, signxx, signyy, signzz, signxy, signyz, signzx, sigy, et, temp, seq, inloc, dplanl, jlag)
for(i8=*sizetab-1;i8 >=0;i8--)