34 1 NEL ,NFUNC ,IFUNC ,NPF ,TF ,
36 3 EPSXX ,EPSYY ,EPSZZ ,EPSXY ,EPSYZ ,
37 4 EPSZX ,OFF ,DFMAX ,TDELE ,
38 5 DELTAX,NUVAR ,UVAR )
45#include "implicit_f.inc"
78 INTEGER NEL, NUPARAM,NUVAR
80 INTEGER,
INTENT(IN) :: NFUNC
81 INTEGER,
DIMENSION(NFUNC),
INTENT(IN) :: IFUNC
82 my_real TIME,UPARAM(*),
83 . epsxx(nel) ,epsyy(nel) ,epszz(nel) ,uvar(nel,nuvar),
84 . epsxy(nel) ,epsyz(nel) ,epszx(nel) ,deltax(nel)
88 my_real off(nel), dfmax(nel),tdele(nel)
93 my_real finter ,tf(*),dydx
105 INTEGER I,II,J,IFAIL,NINDX1,NINDX2,NUM,DEN,
106 . FAIL_ORT,ORDI,,NINDXOFF,IDEL
108 . VOL_STRAIN_LIMIT,MAX_COMP_STRAIN,RATIO_2,,DAMAGE_SP(NEL)
110 . i1,i2,i3,e11,e22,e33,exx,eyy,ezz,exy,ezx,eyz,
111 . q,r,r_inter,phi,ratio(nel),damage(nel),e00,fac,
112 . vol_strain,numerator,denominator,ordinate,
114 INTEGER,
DIMENSION(NEL) :: INDX1,INDX2,INDX3,INDXOFF
155 vol_strain_limit = uparam(1)
158 ordi = int(uparam(4))
159 comp_dir = int(uparam(5))
160 max_comp_strain = uparam(6)
162 idel = int(uparam(8))
166 IF ((uvar(1,1) == zero).AND.(nfunc == 2))
THEN
168 lambda = deltax(i)/el_ref
169 uvar(i,1) = finter(ifunc(2),lambda,npf,tf,dydx)
176 IF (off(i) < one) off(i) = off(i)*four_over_5
177 IF (off(i) < em01) off(i) = zero
194 IF (off(i) == one)
THEN
195 nindxoff = nindxoff + 1
196 indxoff(nindxoff) = i
204#include "vectorize.inc"
223 i2 = exx*eyy + eyy*ezz + ezz*exx - exy*exy - ezx*ezx - eyz*eyz
224 i3 = exx*eyy*ezz - exx*eyz*eyz - eyy*ezx*ezx - ezz*exy*exy + two*exy*ezx*eyz
225 q = (three*i2 - i1*i1)/nine
226 r = (two*i1*i1*i1-nine*i1*i2+twenty7*i3)/cinquante4
227 r_inter =
min(r/sqrt(
max(em20,(-q**3))),one)
228 phi = acos(
max(r_inter,-one))
230 e11 = two*sqrt(
max(-q,zero))*cos( phi/three ) + third*i1
231 e22 = two*sqrt(
max(-q,zero))*cos((phi+two*pi)/three ) + third*i1
232 e33 = two*sqrt(
max(-q,zero))*cos((phi+four*pi)/three) + third*i1
251 vol_strain = e11 + e22 + e33 + e11*e22 + e11*e33 + e22*e33 + e11*e22*e33
256 IF (abs(vol_strain) >= vol_strain_limit)
THEN
261 ELSEIF (num == 2)
THEN
263 ELSEIF (num == 3)
THEN
265 ELSEIF (num == 4)
THEN
267 ELSEIF (num == 5)
THEN
275 denominator = ((exx + ezz)/two) + sqrt
276 ELSEIF (den == 2)
THEN
277 denominator = ((exx + eyy)/two) + sqrt(((exx - eyy)/two)**2 + exy**2)
278 ELSEIF (den == 3)
THEN
279 denominator = ((eyy + ezz)/two) + sqrt(((eyy - ezz)/two)**2 + eyz**2)
280 ELSEIF (den == 4)
THEN
282 ELSEIF (den == 5)
THEN
289 ratio(i) = abs(numerator/
max(denominator,em20))
290 IF (numerator > zero .AND. denominator > zero) ratio(i) = zero
291 IF (numerator < zero .AND. denominator < zero) fail_ort = 1000
292 plamax(i) = finter(ifunc(1),ratio(i),npf,tf,dydx)
293 plamax(i) = plamax(i) * fac
296 IF (fail_ort == 1)
THEN
298 damage(i)=
max(exx/(
max(plamax(i),em20)),
299 . eyy/(
max(plamax(i),em20)),
300 . ezz/(
max(plamax(i),em20)))
302 IF (exx >= plamax(i) .OR. eyy >= plamax(i) .OR. ezz >= plamax(i))
THEN
309 ELSEIF (fail_ort == 2)
THEN
310 damage(i) = exx/(
max(plamax(i),em20))
311 IF (exx >= plamax(i) .AND. exx > zero)
THEN
318 ELSEIF (fail_ort == 3)
THEN
319 damage(i) = eyy/(
max(plamax(i),em20))
320 IF (eyy >= plamax(i) .AND. eyy > zero)
THEN
327 ELSEIF (fail_ort == 4)
THEN
328 damage(i) = ezz/(
max(plamax(i),em20))
329 IF (ezz >= plamax(i) .AND. ezz > zero)
THEN
336 ELSEIF (fail_ort == 5)
THEN
337 damage(i) = e11/(
max(plamax(i),em20))
338 IF (e11 >= plamax(i) .AND. e11 > zero)
THEN
345 ELSEIF (fail_ort == 6)
THEN
346 e00 = ((exx + ezz)/two) + sqrt(((exx - ezz)/two)**2 + ezx**2)
347 damage(i) = e00/(
max(plamax(i),em20))
348 IF (e00 >= plamax(i) .AND. e00 > zero)
THEN
355 ELSEIF (fail_ort == 7)
THEN
356 e00 = ((exx + eyy)/two) + sqrt(((exx - eyy)/two)**2 + exy**2)
357 damage(i) = e00/(
max(plamax(i),em20))
358 IF (e00 >= plamax(i) .AND. e00 > zero)
THEN
365 ELSEIF (fail_ort == 8)
THEN
366 e00 = ((eyy + ezz)/two) + sqrt(((eyy - ezz)/two)**2 + eyz**2)
367 damage(i) = e00/(
max(plamax(i),em20))
368 IF (e00 >= plamax(i) .AND. e00 > zero)
THEN
381 IF (comp_dir > 0 .AND. dfmax(i) < one)
THEN
382 IF (comp_dir == 1)
THEN
383 IF (eyy < (-abs(max_comp_strain)))
THEN
392 ELSEIF (ezz < (-abs(max_comp_strain)*ratio_2))
THEN
402 ELSEIF (comp_dir == 2)
THEN
403 IF (ezz < (-abs(max_comp_strain)))
THEN
412 ELSEIF (exx < (-abs(max_comp_strain)*ratio_2))
THEN
422 ELSEIF (comp_dir == 3)
THEN
423 IF (exx < (-abs(max_comp_strain)))
THEN
432 ELSEIF (eyy < (-abs(max_comp_strain)*ratio_2))
THEN
446 damage_sp(i) = damage(i)
447 IF (dfmax(i) <= one) dfmax(i) =
max(dfmax(i),damage_sp(i))
448 dfmax(i) =
min(dfmax(i),one)
460 WRITE(iout, 1000) ngl(i),plamax(i),ratio(i),time
461 WRITE(istdo,1100) ngl(i),plamax(i),ratio(i),time
462#include "lockoff.inc"
469 IF (max_comp_strain > zero)
THEN
471 WRITE(iout, 2000) ngl(i),time
472 WRITE(istdo,2100) ngl(i),time
473#include "lockoff.inc"
476 WRITE(iout, 2002) ngl(i),time
477 WRITE(istdo,2102) ngl(i),time
478#include "lockoff.inc"
480 IF (indx3(i) == 1)
THEN
482 WRITE(iout, 2010) epsxx(i)
483 WRITE(istdo,2110) epsxx(i)
484#include "lockoff.inc"
486 IF (indx3(i) == 2)
THEN
488 WRITE(iout, 2020) epsyy(i)
489 WRITE(istdo,2120) epsyy(i)
490#include "lockoff.inc"
492 IF (indx3(i) == 3)
THEN
494 WRITE(iout, 2030) epszz(i)
495 WRITE(istdo,2130) epszz(i)
496#include "lockoff.inc"
501 1000
FORMAT(1x,
'DELETE SOLID ELEMENT NUMBER (SAHRAEI) EL :',i10,
502 .
' STRAIN LIMIT OF ',1pe10.3,
'REACHED. AT STRAIN RATIO : ',1pe10.3,
503 .
' AT TIME :',1pe20.13)
504 1100
FORMAT(1x,
'DELETE SOLID ELEMENT NUMBER (SAHRAEI) EL :',i10,
505 .
' STRAIN LIMIT OF ',1pe10.3,
' REACHED. AT STRAIN RATIO : ',1pe10.3,
506 .
' AT TIME :',1pe20.13)
507 2000
FORMAT(1x,
'SOLID ELEMENT NUMBER (SAHRAEI) EL :',i10,
508 .
' MAX. IN-PLANE-COMPRESSION STRAIN REACHED AT TIME :',1pe20.13)
509 2100
FORMAT(1x,
'SOLID ELEMENT NUMBER (SAHRAEI) EL :',i10,
510 .
' MAX. IN-PLANE-COMPRESSION STRAIN REACHED AT TIME :'
511 2002
FORMAT(1x,
'SOLID ELEMENT NUMBER (SAHRAEI) EL :',i10,
512 .
' MAX. IN-PLANE-COMPRESSION STRAIN REACHED AT TIME :',1pe20.13)
513 2102
FORMAT(1x,
'SOLID ELEMENT NUMBER (SAHRAEI) EL :',i10,
514 .
' MAX. IN-PLANE-COMPRESSION STRAIN REACHED AT TIME :',1pe20.13)
516 2010
FORMAT(1x,
'MAX PRESSURE STRAIN IN LOCAL X-DIRECTION REACHED : EPSXX= ',1pe20.13)
517 2110
FORMAT(1x,
'MAX PRESSURE STRAIN IN LOCAL X-DIRECTION REACHED : EPSXX= ',1pe20.13)
519 2020
FORMAT(1x,
'MAX PRESSURE STRAIN IN LOCAL Y-DIRECTION REACHED : EPSYY= ',1pe20.13)
520 2120
FORMAT(1x,
'MAX PRESSURE STRAIN IN LOCAL Y-DIRECTION REACHED : EPSYY= ',1pe20.13)
522 2030
FORMAT(1x,
'MAX PRESSURE STRAIN IN LOCAL Z-DIRECTION REACHED : EPSZZ= ',1pe20.13)
523 2130
FORMAT(1x,
'MAX PRESSURE STRAIN IN LOCAL Z-DIRECTION REACHED : EPSZZ= ',1pe20.13)
subroutine mmain(pm, elbuf_str, ix, nix, x, geo, iparg, nel, skew, bufmat, ipart, ipartel, nummat, matparam, imat, ipm, ngl, pid, npf, tf, mfxx, mfxy, mfxz, mfyx, mfyy, mfyz, mfzx, mfzy, mfzz, rx, ry, rz, sx, sy, sz, gama, voln, dvol, s1, s2, s3, s4, s5, s6, dxx, dyy, dzz, d4, d5, d6, wxx, wyy, wzz)