OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps58c.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!|| sigeps58c ../engine/source/materials/mat/mat058/sigeps58c.F
25!||--- called by ------------------------------------------------------
26!|| mulawc ../engine/source/materials/mat_share/mulawc.f90
27!||--- calls -----------------------------------------------------
28!|| finter ../engine/source/tools/curve/finter.F
29!||--- uses -----------------------------------------------------
30!|| interface_table_mod ../engine/share/modules/table_mod.f
31!|| sensor_mod ../common_source/modules/sensor_mod.F90
32!|| table_mod ../engine/share/modules/table_mod.F
33!||====================================================================
34 SUBROUTINE sigeps58c(
35 1 NEL0 ,NUPARAM,NUVAR ,MFUNC ,KFUNC ,
36 2 NPF ,NPT0 ,IPT ,NSENSOR,
37 2 TF ,TIME ,TIMESTEP,UPARAM ,RHO0 ,
38 3 AREA ,EINT ,THKLY ,NIPARAM ,IPARAM ,
39 4 DEPSXX ,DEPSYY ,DEPSXY ,DEPSYZ ,DEPSZX ,
40 5 EPSXX ,EPSYY ,EPSXY ,EPSYZ ,EPSZX ,
41 6 SIGOXX ,SIGOYY ,SIGOXY ,SIGOYZ ,SIGOZX ,
42 7 SIGNXX ,SIGNYY ,SIGNXY ,SIGNYZ ,SIGNZX ,
43 8 SIGVXX ,SIGVYY ,SIGVXY ,SIGVYZ ,SIGVZX ,
44 9 SOUNDSP,VISCMAX,THK ,PLA ,UVAR ,
45 A OFF ,NGL ,PM ,MAT ,ETSE ,
46 B SHF ,YLD ,TAN_PHI ,ALDT ,
47 C SENSOR_TAB,ISMSTR,TABLE ,OFFGG )
48C-----------------------------------------------
49C M o d u l e s
50C-----------------------------------------------
51 USE table_mod
52 USE interface_table_mod
53 USE sensor_mod
54C-----------------------------------------------
55C I m p l i c i t T y p e s
56C-----------------------------------------------
57#include "implicit_f.inc"
58C-----------------------------------------------
59C G l o b a l P a r a m e t e r s
60C-----------------------------------------------
61#include "mvsiz_p.inc"
62C---------+---------+---+---+--------------------------------------------
63C VAR | SIZE |TYP| RW| DEFINITION
64C---------+---------+---+---+--------------------------------------------
65C NEL0 | 1 | I | R | SIZE OF THE ELEMENT GROUP NEL0
66C NUPARAM | 1 | I | R | SIZE OF THE USER PARAMETER ARRAY
67C NUVAR | 1 | I | R | NUMBER OF USER ELEMENT VARIABLES
68C---------+---------+---+---+--------------------------------------------
69C NFUNC | 1 | I | R | NUMBER FUNCTION USED FOR THIS USER LAW
70C IFUNC | NFUNC | I | R | FUNCTION INDEX
71C NPF | * | I | R | FUNCTION ARRAY
72C NPT0 | 1 | I | R | NUMBER OF LAYERS OR INTEGRATION POINTS
73C IPT | 1 | I | R | LAYER OR INTEGRATION POINT NUMBER
74C TF | * | F | R | FUNCTION ARRAY
75C---------+---------+---+---+--------------------------------------------
76C TIME | 1 | F | R | CURRENT TIME
77C TIMESTEP| 1 | F | R | CURRENT TIME STEP
78C UPARAM | NUPARAM | F | R | USER MATERIAL PARAMETER ARRAY
79C RHO0 | NEL0 | F | R | INITIAL DENSITY
80C AREA | NEL0 | F | R | AREA
81C EINT | 2*NEL0 | F | R | INTERNAL ENERGY(MEMBRANE,BENDING)
82C THKLY | NEL0 | F | R | LAYER THICKNESS
83C DEPSXX | NEL0 | F | R | STRAIN INCREMENT XX
84C DEPSYY | NEL0 | F | R | STRAIN INCREMENT YY
85C ... | | | |
86C EPSXX | NEL0 | F | R | STRAIN XX TRUE
87C EPSYY | NEL0 | F | R | STRAIN YY TRUE
88C ... | | | |
89C SIGOXX | NEL0 | F | R | OLD ELASTO PLASTIC STRESS XX
90C SIGOYY | NEL0 | F | R | OLD ELASTO PLASTIC STRESS YY
91C ... | | | |
92C---------+---------+---+---+--------------------------------------------
93C SIGNXX | NEL0 | F | W | NEW ELASTO PLASTIC STRESS XX
94C SIGNYY | NEL0 | F | W | NEW ELASTO PLASTIC STRESS YY
95C ... | | | |
96C SIGVXX | NEL0 | F | W | VISCOUS STRESS XX
97C SIGVYY | NEL0 | F | W | VISCOUS STRESS YY
98C ... | | | |
99C SOUNDSP | NEL0 | F | W | SOUND SPEED (NEEDED FOR TIME STEP)
100C VISCMAX | NEL0 | F | W | MAXIMUM DAMPING MODULUS(NEEDED FOR TIME STEP)
101C---------+---------+---+---+--------------------------------------------
102C THK | NEL0 | F |R/W| THICKNESS
103C PLA | NEL0 | F |R/W| PLASTIC STRAIN
104C UVAR |NEL0*NUVAR| F |R/W| USER ELEMENT VARIABLE ARRAY
105C OFF | NEL0 | F |R/W| DELETED ELEMENT FLAG (=1. ON, =0. OFF)
106C---------+---------+---+---+--------------------------------------------
107C C o m m o n B l o c k s
108C-----------------------------------------------
109#include "param_c.inc"
110#include "com04_c.inc"
111#include "com01_c.inc"
112#include "com08_c.inc"
113C-----------------------------------------------
114C I N P U T A r g u m e n t s
115C-----------------------------------------------
116 INTEGER ,INTENT(IN) :: NSENSOR
117 INTEGER ,INTENT(IN) :: NIPARAM
118 INTEGER ,INTENT(IN) :: NUPARAM
119 INTEGER :: NEL0, NUVAR,NPT0,IPT,ISMSTR
120 INTEGER :: NGL(NEL0),MAT(NEL0)
121 INTEGER ,DIMENSION(NIPARAM) ,INTENT(IN) :: IPARAM
122 my_real ,DIMENSION(NUPARAM) ,INTENT(IN) :: UPARAM
123 my_real TIME,TIMESTEP(NEL0),
124 . AREA(NEL0),RHO0(NEL0),EINT(NEL0,2),
125 . THKLY(NEL0),PLA(NEL0),OFFGG(NEL0),TAN_PHI(NEL0),
126 . DEPSXX(NEL0),DEPSYY(NEL0),
127 . DEPSXY(NEL0),DEPSYZ(NEL0),DEPSZX(NEL0),
128 . EPSXX(NEL0) ,EPSYY(NEL0) ,
129 . epsxy(nel0) ,epsyz(nel0) ,epszx(nel0) ,
130 . sigoxx(nel0),sigoyy(nel0),
131 . sigoxy(nel0),sigoyz(nel0),sigozx(nel0),
132 . pm(npropm,*),shf(*),aldt(nel0)
133C-----------------------------------------------
134C O U T P U T A r g u m e n t s
135C-----------------------------------------------
136 my_real
137 . signxx(nel0),signyy(nel0),
138 . signxy(nel0),signyz(nel0),signzx(nel0),
139 . sigvxx(nel0),sigvyy(nel0),
140 . sigvxy(nel0),sigvyz(nel0),sigvzx(nel0),
141 . soundsp(nel0),viscmax(nel0),etse(nel0)
142C-----------------------------------------------
143C I N P U T O U T P U T A r g u m e n t s
144C-----------------------------------------------
145 my_real uvar(nel0,nuvar), off(nel0),thk(nel0),yld(nel0)
146 TYPE(ttable) TABLE(*)
147 my_real, DIMENSION(1) :: XX
148 TYPE (SENSOR_STR_) ,DIMENSION(NSENSOR) ,INTENT(IN) :: SENSOR_TAB
149C-----------------------------------------------
150C VARIABLES FOR FUNCTION INTERPOLATION
151C-----------------------------------------------
152 INTEGER NPF(*), MFUNC, KFUNC(MFUNC)
153 my_real FINTER ,TF(*)
154 EXTERNAL FINTER
155C Y = FINTER(IFUNC(J),X,NPF,TF,DYDX)
156C Y : y = f(x)
157C X : x
158C DYDX : f'(x) = dy/dx
159C IFUNC(J): FUNCTION INDEX
160C J : FIRST(J=1), SECOND(J=2) .. FUNCTION USED FOR THIS LAW
161C NPF,TF : FUNCTION PARAMETER
162C-----------------------------------------------
163C L o c a l V a r i a b l e s
164C-----------------------------------------------
165 INTEGER I,J,NC,NT,ISENS,ITER,NITER,NINDX,NUMOFF,UNLOAD,IERR,
166 . ifuncc,ifunct,idxcnt,idx1,idx2,ierrt,ierrc,index(mvsiz),
167 . indxoff(mvsiz),tagoldc,tagoldt,tagc(0:1),tagt(0:1)
168 my_real
169 . kc,kt,kfc,kft,g0,gt,gb,tan_lock,visce,viscg,lc0,lt0,
170 . kbc,kbt,dc0,dt0,hc0,ht0,ec0,et0,dc,dt,trace,gfrot,gsh,
171 . ac,at,bc,bt,bc2,bt2,dcc,dtt,func,funt,deric,derit,dtinv,
172 . tfrot,sigg,y,ec2,et2,sigc,sigt,sig0,udc,udt,fc,ft,
173 . fpc,fpt,ccl,ttl,dec,det,hc,ht,kf,kg,kmax,hdc,hdt,
174 . sinn,tan2,damp,rfac,rfat,v1,v2,areamin,dareamin,
175 . aa,vv,zerostress,dsig,lmin,tstart,yfac(6),yy,phi,gxy,
176 . flex1,flex2,facc,fact,kfc0,kft0,etc,ett,phin,phii,sxyi,
177 . ddec,ddet,coeful,coefrl,dcc1,dtt1,dcc0,dtt0,phii2,sxyi2,
178 . epsi1,sigi1,epsi2,sigi2,kkc,kkt,eminc,emaxc,sminc,smaxc,
179 . emint,emaxt,smint,smaxt,face_c,face_t,facs_c,facs_t,
180 . face,facf,fate,fatf,fcl,ftl,fclp,ftlp,test,dphi,phio,kc0,kt0,
181 . dyc,dyt,mass
182 my_real
183 . ec(mvsiz),et(mvsiz),lc(mvsiz),lt(mvsiz),fn(mvsiz),sigg0(mvsiz),
184 . dtang(mvsiz),tfold(mvsiz),yc(mvsiz),yt(mvsiz),
185 . eminrlc(mvsiz),emaxrlc(mvsiz),epsmaxc(mvsiz),
186 . eminrlt(mvsiz),emaxrlt(mvsiz),epsmaxt(mvsiz),
187 . epsnc(mvsiz),epsnt(mvsiz),epsns(mvsiz),
188 . epsnrlc(mvsiz),epsnrlt(mvsiz),smaxrlc(mvsiz),smaxrlt(mvsiz),
189 . sminrlc(mvsiz),sigmaxc(mvsiz),sminrlt(mvsiz),sigmaxt(mvsiz),
190 . phimin(mvsiz),phimax(mvsiz),phirlmax(mvsiz),sxymax(mvsiz),
191 . sxymin(mvsiz),sxymaxrl(mvsiz),cvisc(mvsiz),cvist(mvsiz)
192C-----------------------------------------------
193C S t a t e V a r i a b l e s (U V A R)
194C-----------------------------------------------
195c UVAR(1) = SIGNXX + SIGVXX ! total engineering stress in 1st direction (elastic + viscous)
196c UVAR(2) = SIGNYY + SIGVYY ! total engineering stress in 2nd direction (elastic + viscous)
197c UVAR(3) = SIGNXY + SIGVXY ! total shear stress (elastic + friction)
198c UVAR(4) = EC ! total true element strain in 1st fiber direction
199c UVAR(5) = ET ! total true element strain in 2nd fiber direction
200c UVAR(6) = DEPSXY ! tan(alpha) of shear angle
201c UVAR(7) = YC ! total flex element elongation in 1st direction
202c UVAR(8) = YT ! total flex element elongation in 2nd direction
203c UVAR(9) = SIGVXY ! Friction stress (shear)
204c UVAR(10) = SIG0 ! Initial shear stress (if initial shear angle > 0)
205c UVAR(11) = SIGNXX ! elastic stress in 1st direction (used for zerostress relaxation)
206c UVAR(12) = SIGNYY ! elastic stress in 2nd direction (used for zerostress relaxation)
207c UVAR(13) = SIGNXY ! elastic stress in shear (used for zerostress relaxation)
208c UVAR(14) = ALDT ! initial characteristic length of the element
209c UVAR(15) = DC ! total fiber elongation in 1st direction
210c UVAR(16) = DT ! total fiber elongation in 2nd direction
211C======================================================================|
212C--- Initialisations
213C----------------------------------------------------------------------
214 niter = 3
215!
216 unload = iparam(1)
217 isens = iparam(2)
218 nc = iparam(3)
219 nt = iparam(4)
220!
221 lc0 = uparam( 1)
222 lt0 = uparam( 2)
223 dc0 = uparam( 3)
224 dt0 = uparam( 4)
225 hc0 = uparam( 5)
226 ht0 = uparam( 6)
227 kc0 = uparam( 9) ! young dir1
228 kt0 = uparam(10) ! young dir2
229 kkc = uparam(40) ! initial slope dir1
230 kkt = uparam(41) ! initial slope dir2
231 kfc0 = uparam(11)
232 kft0 = uparam(12)
233 g0 = uparam(13)
234 gt = uparam(14)
235 gb = uparam(15)
236 tan_lock = uparam(16)
237 visce = uparam(17)
238 viscg = uparam(18)
239 kbc = uparam(19)
240 kbt = uparam(20)
241 gfrot = uparam(21)
242 areamin = uparam(22)
243 dareamin = uparam(23)
244 zerostress = uparam(24)
245 flex1 = uparam(26)
246 flex2 = uparam(27)
247 gsh = uparam(32)*shf(1)
248 IF (zerostress > zero .and. isens > 0) THEN
249 tstart = sensor_tab(isens)%TSTART
250 ELSE
251 tstart = zero
252 ENDIF
253c
254 yfac(1) = uparam(27 + 1 )
255 yfac(2) = uparam(27 + 2 )
256 yfac(3) = uparam(27 + 3 )
257 yfac(4) = uparam(33)
258 yfac(5) = uparam(34)
259 yfac(6) = uparam(42)
260c
261 kfc = kfc0
262 kft = kft0
263 kc = kc0
264 kt = kt0
265 nindx = 0
266 numoff = 0
267C-----------------------------------------------------------
268C DEACTIVATION FOR VERY WARPED ELEMENT (QBAT)
269C-----------------------------------------------------------
270 IF (npt0 ==1) THEN
271 DO i=1,nel0
272 IF (uvar(i,40)==one) cycle
273 depsxx(i) = depsxx(i)*uvar(i,40)
274 depsyy(i) = depsyy(i)*uvar(i,40)
275 aa =uvar(i,10)/g0
276 depsxy(i) = aa*(one-uvar(i,40))+depsxy(i)*uvar(i,40)
277 ENDDO
278 ENDIF
279C
280 IF (unload == 0) THEN !no hysteresis
281c
282 ! KF = KFC + KFT
283 ! CCL : eps value when derivative = 0
284 IF (kbc == zero) THEN
285 ccl = ep30
286 ELSE
287 ccl = kc / kbc ! condition pour derivee nulle
288 ENDIF
289
290 IF (kbt == zero) THEN
291 ttl = ep30
292 ELSE
293 ttl = kt / kbt
294 ENDIF
295 DO i=1,nel0
296 IF(offgg(i) > zero) THEN
297 nindx = nindx+1
298 index(nindx)=i
299 ELSE
300 numoff = numoff+1
301 indxoff(numoff)=i
302 ENDIF
303 ENDDO
304#include "vectorize.inc"
305 DO j=1,nindx
306 i=index(j)
307 yc(i) = uvar(i,7)
308 yt(i) = uvar(i,8)
309 tfold(i) = uvar(i,9)
310 mass = rho0(i)*area(i)*thkly(i)*fourth ! 2m = rho*V/4
311 dtinv = timestep(i)/max(timestep(i)*timestep(i),em20)
312 cvisc(i) = sqrt(mass*kfc0)*dtinv*third
313 cvist(i) = sqrt(mass*kft0)*dtinv*third
314 ENDDO
315
316C--- integration des deformations
317#include "vectorize.inc"
318 DO j=1,nindx
319 i=index(j)
320 etc =uvar(i,4)+depsxx(i)
321 ett =uvar(i,5)+depsyy(i)
322 uvar(i,4) = etc
323 uvar(i,5) = ett
324 ec(i) = exp(etc) - one ! eng strain
325 et(i) = exp(ett) - one
326 ENDDO
327#include "vectorize.inc"
328 DO j=1,nindx
329 i=index(j)
330 lc(i) = lc0 * (one + ec(i))
331 lt(i) = lt0 * (one + et(i))
332 ENDDO
333c
334c--- Calcul YC, YT
335 ac = kc*dc0
336 ac = ac*ac
337 at = kt*dt0
338 at = at*at
339#include "vectorize.inc"
340 DO j=1,nindx
341 i=index(j)
342 fn(i) = zero
343c--- uncoupled model (compression of fiber)
344 dyc = zero
345 dyt = zero
346 DO iter = 1, niter
347 hc = hc0 + yc(i) !longueur ressort
348 ht = ht0 + yt(i)
349 dc = sqrt(lc(i)*lc(i) + hc*hc) ! long fibre
350 dt = sqrt(lt(i)*lt(i) + ht*ht)
351 udc = one / dc
352 udt = one / dt
353 hdc = hc * udc !sinus(alpha)
354 hdt = ht * udt
355 dcc = dc - dc0
356 dtt = dt - dt0
357 IF (kfunc(1) /= 0) THEN
358 fc = yfac(1)*finter(kfunc(1),dcc,npf,tf,fpc)
359 fpc = fpc *yfac(1)
360 kc = fpc
361 kfc = flex1*fpc*hc0/dc0 ! rigidite de flexion
362 IF(kfc == zero)kfc=kfc0
363 fpc = fpc*hdc
364 ELSEIF (dcc >= ccl ) THEN
365 fc = half * kc*ccl
366 fpc = zero
367 ELSE
368 fc = (kc - half*kbc * dcc) * dcc
369 fpc = (kc - kbc*dcc) * hdc! derivee % yc
370 ENDIF
371 IF (kfunc(2) /= 0) THEN
372 ft = yfac(2)*finter(kfunc(2),dtt,npf,tf,fpt)
373 fpt = fpt * yfac(2)
374 kt = fpt
375 kft = flex2*fpt*ht0/dt0
376 IF(kft == zero)kft=kft0
377 fpt = fpt*hdt
378 ELSEIF (dtt >= ttl ) THEN
379 ft = half * kt*ttl
380 fpt = zero
381 ELSE
382 ft = (kt - half*kbt * dtt) * dtt
383 fpt = (kt - kbt*dtt) * hdt
384 ENDIF
385 func = kfc*yc(i) + fc*hc/dc + cvisc(i)*dyc
386 funt = kft*yt(i) + ft*ht/dt + cvist(i)*dyt
387 deric= kfc + fpc*hdc + fc*udc*(one - hdc*hdc) + cvisc(i)
388 derit= kft + fpt*hdt + ft*udt*(one - hdt*hdt) + cvist(i)
389C
390 yc(i) = yc(i) - func / deric
391 yt(i) = yt(i) - funt / derit
392 dyc = yc(i) - uvar(i,7)
393 dyt = yt(i) - uvar(i,8)
394 fn(i) = zero
395 ENDDO !iter
396C--- Coupled model (traction fibre)
397 dyc = zero
398 dyt = zero
399 IF ((yc(i) + yt(i)) < zero) THEN !eps
400 y = half * (uvar(i,7) - uvar(i,8))
401 DO iter = 1, niter
402 hc = hc0 + y
403 ht = ht0 - y
404 dc = sqrt(lc(i)*lc(i) + hc*hc)
405 dt = sqrt(lt(i)*lt(i) + ht*ht)
406 dcc = dc - dc0
407 dtt = dt - dt0
408 udc= one / dc
409 udt= one / dt
410 hdc= hc * udc
411 hdt= ht * udt
412 IF (kfunc(1) /= 0) THEN
413 fc =yfac(1)* finter(kfunc(1),dcc,npf,tf,fpc)
414 fpc = fpc *yfac(1)
415 kc = fpc
416 kfc = flex1*fpc*hc0/dc0
417 IF(kfc == zero)kfc=kfc0
418 fpc = fpc*hdc
419 ELSEIF (dcc >= ccl ) THEN
420 fc = half * kc*ccl
421 fpc = zero
422 ELSE
423 fc = (kc - half*kbc * dcc) * dcc
424 fpc = (kc - kbc*dcc) * hdc
425 ENDIF
426 IF (kfunc(2) /= 0) THEN
427 ft = yfac(2)*finter(kfunc(2),dtt,npf,tf,fpt)
428 fpt = fpt * yfac(2)
429 kt = fpt
430 kft = flex2*fpt*ht0/dt0
431 IF(kft == zero)kft=kft0
432 fpt = fpt*hdt
433 ELSEIF (dtt >= ttl ) THEN
434 ft = half * kt*ttl
435 fpt = zero
436 ELSE
437 ft = (kt - half*kbt * dtt) * dtt
438 fpt = (kt - kbt*dtt) * hdt
439 ENDIF
440C
441 kf = kfc + kft
442 func = kf*y + fc * hdc - ft * hdt
443 . + (cvisc(i) + cvist(i))*dyc
444 deric = kf + fpc*hdc + fc*udc*(one - hdc*hdc)
445 . + fpt*hdt + ft*udt*(one - hdt*hdt)
446 . + cvisc(i) + cvist(i)
447 y = y - func / deric
448 dyc = y - half * (uvar(i,7) - uvar(i,8))
449C
450 IF (y > 0) THEN
451 y = min(y, ht0)
452 ELSE
453 y = max(y,-hc0)
454 ENDIF
455 ENDDO !iter
456c
457 yc(i) = y
458 yt(i) =-y
459 fn(i) = fc*hc/dc + ft*ht/dt
460 ENDIF
461 hc = hc0 + yc(i)
462 ht = ht0 + yt(i)
463 dc = sqrt(lc(i)*lc(i) + hc*hc)
464 dt = sqrt(lt(i)*lt(i) + ht*ht)
465 ENDDO !i,nel
466C-------------------------------------
467#include "vectorize.inc"
468 DO j=1,nindx
469 i=index(j)
470 hc = hc0 + yc(i)
471 ht = ht0 + yt(i)
472 dc = sqrt(lc(i)*lc(i) + hc*hc)
473 dt = sqrt(lt(i)*lt(i) + ht*ht)
474 dcc = dc - dc0
475 dtt = dt - dt0
476 IF (kfunc(1) /= 0) THEN
477 fc =yfac(1)* finter(kfunc(1),dcc,npf,tf,fpc)
478 fpc = fpc *yfac(1)
479 kc = fpc
480 ELSEIF (dcc >= ccl ) THEN
481 fc = half * kc*ccl
482 ELSE
483 fc = (kc - half*kbc * dcc) * dcc
484 ENDIF
485 IF(kfunc(2) /= 0)THEN
486 ft = yfac(2)*finter(kfunc(2),dtt,npf,tf,fpt)
487 fpt = fpt * yfac(2)
488 kt = fpt
489 ELSEIF (dtt >= ttl ) THEN
490 ft = half * kt*ttl
491 ELSE
492 ft = (kt - half*kbt * dtt) * dtt
493 ENDIF
494C------------------------------------------------------------------
495C Trace = Eps1 + Eps2 (Trace of true principal strain tensor)
496C Trace = ec_true + ec2_true
497C ec2_eng = exp(Tr) / (ec_eng + 1) - 1
498C rfac = 2*Nc / (ec2_eng+1)
499C ec2_eng+1 > 0
500C------------------------------------------------------------------
501 trace = exp(epsxx(i) + epsyy(i))! =exp(tr)
502 ec2 = max(trace / (ec(i) + one), em6)
503 et2 = max(trace / (et(i) + one), em6)
504 rfac = nc / ec2
505 rfat = nt / et2
506C--- contraintes membrane
507 sigc = fc * lc(i) / dc
508 sigt = ft * lt(i) / dt
509 signxx(i) = sigc * rfac
510 signyy(i) = sigt * rfat
511c
512 uvar(i,1) = sigc
513 uvar(i,2) = sigt
514 uvar(i,15) = dc
515 uvar(i,16) = dt
516C--- contraintes cisaillement
517C------------------------------------------------------------------
518C------------------------------------------------------------------
519c++++++++++
520c G0 = GT / (ONE + TAN_LOCK*TAN_LOCK)
521c SIGG = G0*TAN_PHI - SIGG0
522c KG = G0*TAN_PHI*TAN_PHI
523c++++++++++
524 sig0 = uvar(i,10)
525 tan_phi(i)= depsxy(i)
526 tan2 = tan_phi(i)*tan_phi(i)
527 IF (kfunc(3) /= 0)THEN
528 phi = atan(tan_phi(i))*hundred80/pi
529 signxy(i) = yfac(3)* finter(kfunc(3),phi,npf,tf,gxy)
530 kg = gxy*yfac(3)*tan2
531 ELSEIF (tan_phi(i) > tan_lock) THEN
532 signxy(i) = gt*tan_phi(i) + gb - sig0
533 kg = gt*tan2
534 ELSEIF (tan_phi(i) < -tan_lock) THEN
535 signxy(i) = gt*tan_phi(i) - gb - sig0
536 kg = gt*tan2
537 ELSE
538 signxy(i) = g0*tan_phi(i) - sig0
539 kg = g0*tan2
540 ENDIF
541C--- contraintes visc fibres
542 dtinv = timestep(i)/max(timestep(i)*timestep(i),em20)
543 damp = sqrt(rho0(i)*area(i)*thkly(i)*half)
544 v1 = visce*damp*sqrt(nc*kc)
545 v2 = visce*damp*sqrt(nt*kt)
546 sigvxx(i) = dtinv*(depsxx(i))*v1
547 sigvyy(i) = dtinv*(depsyy(i))*v2
548C--- contraintes frottement cisaillement
549 IF (fn(i) > zero) THEN
550 tfrot = two_third*viscg*fn(i)*(hc0+ht0)/(lc(i)+lt(i))
551 dtang(i) = depsxy(i) - tan_phi(i)
552 sigg = tfold(i) + gfrot*dtang(i)
553 IF (abs(sigg) > tfrot) THEN
554 sigvxy(i) = sign(tfrot,sigg)
555 ELSE
556 sigvxy(i) = sigg
557 ENDIF
558 ENDIF
559C---
560 sinn = tan_phi(i) / sqrt(one + tan2)
561 kmax = max(kc0*rfac,kt0*rfat)*(one+sinn) + kg
562 lmin = min(dc/dc0,dt/dt0)*uvar(i,14)
563c if (RHO0(I)==ZERO.OR.LMIN==ZERO)
564 soundsp(i) = sqrt(kmax/(rho0(i)))*aldt(i)/lmin
565 viscmax(i) = max(v1,v2)
566 etse(i) = one
567 ENDDO
568C---
569#include "vectorize.inc"
570 DO j=1,nindx
571 i=index(j)
572 uvar(i,3) = signxy(i) + sigvxy(i)
573 tan_phi(i)= depsxy(i)
574 uvar(i,6) = depsxy(i)
575 uvar(i,7) = yc(i)
576 uvar(i,8) = yt(i)
577 uvar(i,9) = sigvxy(i)
578c
579 signyz(i) = sigoyz(i) + gsh *depsyz(i)
580 signzx(i) = sigozx(i) + gsh *depszx(i)
581 ENDDO
582!=========================================================
583!=========================================================
584 ELSE ! LOADING/UNLOADING WITH HYSTERESIS (UNLOAD = 1)
585!=========================================================
586!=========================================================
587 niter = 5
588 epsi1 = uparam(36)
589 sigi1 = uparam(37)
590 epsi2 = uparam(38)
591 sigi2 = uparam(39)
592 phii = uparam(43)
593 sxyi = uparam(44)
594 phii2 = uparam(45)
595 sxyi2 = uparam(46)
596 IF (isigi == 0) THEN
597 IF (time == zero)THEN
598 DO i=1,nel0
599 uvar(i,15)= dc0
600 uvar(i,16)= dt0
601 ENDDO
602 ENDIF
603 ENDIF
604 DO i=1,nel0
605 IF(offgg(i) > zero) THEN
606 nindx = nindx+1
607 index(nindx)=i
608 ELSE
609 numoff = numoff+1
610 indxoff(numoff)=i
611 ENDIF
612 ENDDO
613#include "vectorize.inc"
614 DO j=1,nindx
615 i=index(j)
616 yc(i) = uvar(i,7)
617 yt(i) = uvar(i,8)
618 tfold(i) = uvar(i,9)
619 ENDDO
620C--- integration des deformations
621#include "vectorize.inc"
622 DO j=1,nindx
623 i=index(j)
624 etc =uvar(i,4)+depsxx(i)
625 ett =uvar(i,5)+depsyy(i)
626 uvar(i,4) = etc
627 uvar(i,5) = ett
628 ec(i) = exp(etc) - one ! eng strain
629 et(i) = exp(ett) - one
630 ENDDO
631#include "vectorize.inc"
632 DO j=1,nindx
633 i=index(j)
634 lc(i) = one + ec(i)
635 lt(i) = one + et(i)
636 epsmaxc(i) = uvar(i,17)
637 eminrlc(i) = uvar(i,18)
638 emaxrlc(i) = uvar(i,19)
639 sigmaxc(i) = uvar(i,20)
640 sminrlc(i) = uvar(i,25)
641 smaxrlc(i) = uvar(i,27)
642 epsmaxt(i) = uvar(i,21)
643 eminrlt(i) = uvar(i,22)
644 emaxrlt(i) = uvar(i,23)
645 sigmaxt(i) = uvar(i,24)
646 sminrlt(i) = uvar(i,26)
647 smaxrlt(i) = uvar(i,28)
648 ENDDO
649#include "vectorize.inc"
650 DO j=1,nindx
651 i=index(j)
652 fn(i) = zero
653c--- uncoupled model (compression of fiber)
654 hc = hc0 + yc(i) !longueur ressort
655 ht = ht0 + yt(i)
656 dc = sqrt(lc(i)*lc(i) + hc*hc) ! long fibre
657 dt = sqrt(lt(i)*lt(i) + ht*ht)
658 dcc = uvar(i,15) - dc0
659 dtt = uvar(i,16) - dt0
660 ddec = dc - uvar(i,15)
661 ddet = dt - uvar(i,16)
662 IF (dcc >zero ) THEN
663 IF ((ddec >=zero .AND. dcc> = uvar(i,17)) .OR.
664 . (ddec >=zero .AND. uvar(i,17)==uvar(i,18)) ) THEN
665 !----------------------LOAD
666 DO iter = 1, niter
667 hc = hc0 + yc(i) !longueur ressort
668 dc = sqrt(lc(i)*lc(i) + hc*hc) ! long fibre
669 dcc = dc - dc0
670 dcc = max(uvar(i,18),dcc)
671 dc = dcc + dc0
672 udc = one / dc
673 hdc = hc * udc !sinus(alpha)
674 IF (kfunc(1) == 0) THEN
675 fc = kc0 * dcc
676 fpc = kc0 * hdc
677 ELSE
678 fc = yfac(1) * finter(kfunc(1),dcc,npf,tf,fpc)
679 fpc = fpc * yfac(1)
680 kc = fpc
681 fpc = fpc * hdc
682 ENDIF
683 epsmaxc(i) = dcc
684 sigmaxc(i) = fc
685 smaxrlc(i) = sigmaxc(i)
686 emaxrlc(i) = epsmaxc(i)! UVAR(I,19))
687 func = kfc*yc(i) + fc*hc/dc
688 deric= kfc + fpc*hdc + fc*udc*(one - hdc*hdc)
689 yc(i) = yc(i) - func / deric
690 eminrlc(i) = dcc
691 sminrlc(i) = fc
692 ENDDO !iter
693 !TAGOLDC = 1
694 ELSEIF(ddec >=zero ) THEN
695 !----------------------RELOAD
696 DO iter = 1, niter
697 hc = hc0 + yc(i) !longueur ressort
698 dc = sqrt(lc(i)*lc(i) + hc*hc) ! long fibre
699 dcc = dc - dc0
700 dcc = max(uvar(i,18),dcc)
701 dc = dcc + dc0
702 udc = one / dc
703 hdc = hc * udc !sinus(alpha)
704 epsnrlc(i) = uvar(i,17)*(dcc - uvar(i,18)) /(uvar(i,17)-uvar(i,18))
705 IF (kfunc(1) == 0) THEN
706 fc = kc0 * dcc
707 fpc = kc0 * hdc
708 ELSE
709 fc = yfac(1) * finter(kfunc(1),epsnrlc(i),npf,tf,fpc)
710 coefrl = (uvar(i,20) - uvar(i,25))/uvar(i,20)
711 fc = uvar(i,25) + coefrl * fc
712 fpc = fpc * yfac(1)
713 fpc = fpc * uvar(i,17)* (uvar(i,20) - uvar(i,25))
714 . /(uvar(i,17)-uvar(i,18))/uvar(i,20)
715 kc = fpc
716 fpc = fpc * hdc
717 ENDIF
718 emaxrlc(i)= dcc ! UVAR(I,19))
719 smaxrlc(i) = fc
720 func = kfc*yc(i) + fc*hc/dc
721 deric= kfc + fpc*hdc + fc*udc*(one - hdc*hdc)
722 yc(i) = yc(i) - func / deric
723 ENDDO !iter
724 IF(dcc>uvar(i,17))then!emax
725 epsmaxc(i) = dcc
726 sigmaxc(i) = fc
727 ENDIF
728 !TAGOLDC = 2
729 ELSE
730c
731 IF (uvar(i,19) > zero) THEN !emaxrl
732 DO iter = 1, niter
733 hc = hc0 + yc(i) !longueur ressort
734 dc = sqrt(lc(i)*lc(i) + hc*hc) ! long fibre
735 dcc = dc - dc0
736 dcc = max(zero,dcc)
737 dc = dcc + dc0
738 udc = one / dc
739 hdc = hc * udc !sinus(alpha)
740 epsnc(i) = dcc*epsi1/uvar(i,19)
741 IF (kfunc(1) == 0) THEN
742 fc = kc0 * dcc
743 fpc = kc0 * hdc
744 ELSE
745 fc = yfac(4) * finter(kfunc(4),epsnc(i),npf,tf,fpc)
746 fpc = fpc*yfac(4)
747 coeful = fc/sigi1
748 fc = coeful *uvar(i,27) !SMAXRLC(I)
749 fpc = fpc*uvar(i,27)*epsi1/uvar(i,19)/sigi1
750 kc = fpc
751 fpc = fpc * hdc
752 ENDIF
753 eminrlc(i)= dcc
754 sminrlc(i) = fc
755 func = kfc*yc(i) + fc*hc/dc
756 deric= kfc + fpc*hdc + fc*udc*(one - hdc*hdc)
757 yc(i) = yc(i) - func / deric
758 ENDDO !iter
759 else!(UVAR(I,19) >ZERO
760 DO iter = 1, niter
761 hc = hc0 + yc(i) !longueur ressort
762 dc = sqrt(lc(i)*lc(i) + hc*hc) ! long fibre
763 udc = one / dc
764 hdc = hc * udc !sinus(alpha)
765 dcc = dc - dc0
766 fc = kkc * dcc
767 fpc = kkc * hdc! derivee % yc
768 kc = kkc
769 func = kfc*yc(i) + fc*hc/dc
770 deric= kfc + fpc*hdc + fc*udc*(one - hdc*hdc)
771 yc(i) = yc(i) - func / deric
772 fn(i) = zero
773 ENDDO !iter
774 epsmaxc(i) = zero
775 sigmaxc(i) = zero
776 smaxrlc(i) = zero
777 emaxrlc(i)= zero
778 eminrlc(i) = zero
779 sminrlc(i) = zero
780 endif!(UVAR(I,19) >ZERO)
781 !TAGOLDC = 3
782 ENDIF !LOAD RELOD OR UNLOAD
783 ELSE ! DCC < ZERO domaine de compression
784 DO iter = 1, niter
785 hc = hc0 + yc(i) !longueur ressort
786 dc = sqrt(lc(i)*lc(i) + hc*hc) ! long fibre
787 udc = one / dc
788 hdc = hc * udc !sinus(alpha)
789 dcc = dc - dc0
790 ! compression fibre
791 fc = kkc * dcc
792 kc = kkc
793 fpc = kkc * hdc! derivee % yc
794 func = kfc*yc(i) + fc*hc/dc
795 deric= kfc + fpc*hdc + fc*udc*(one - hdc*hdc)
796 yc(i) = yc(i) - func / deric
797 fn(i) = zero
798 ENDDO !iter
799 epsmaxc(i) = zero
800 eminrlc(i) = zero
801 emaxrlc(i) = zero
802 sigmaxc(i) = zero
803 sminrlc(i) = zero
804 smaxrlc(i) = zero
805 ENDIF ! COMPRESSION OR TENSION FOR DIRECTION CHAINE
806c----------------------------------
807c----------------direction trame---
808c----------------------------------
809 IF(dtt >zero ) THEN !(DDET >=ZERO .AND. UVAR(I,22)==ZERO).OR.
810 IF((ddet >=zero .AND.dtt>uvar(i,21))
811 . .OR.(ddet >=zero .AND.uvar(i,21)==uvar(i,22))) THEN
812 !----------------------LOAD
813 DO iter = 1, niter
814 ht = ht0 + yt(i) !longueur ressort
815 dt = sqrt(lt(i)*lt(i) + ht*ht) ! long fibre
816 dtt = dt - dt0
817 dtt = max(uvar(i,22),dtt)
818 dt = dtt + dt0
819 udt = one / dt
820 hdt = ht * udt !sinus(alpha)
821 IF (kfunc(2) == 0) THEN
822 ft = kt0 * dtt
823 fpt = kt0 * hdt
824 ELSE
825 ft = yfac(2)*finter(kfunc(2),dtt,npf,tf,fpt)
826 fpt = fpt * yfac(2)
827 kt = fpt
828 fpt = fpt *hdt
829 ENDIF
830 epsmaxt(i) = dtt
831 sigmaxt(i) = ft
832 smaxrlt(i) = sigmaxt(i)
833 emaxrlt(i)= epsmaxt(i)!
834 eminrlt(i)= zero
835 sminrlt(i) = zero
836 funt = kft*yt(i) + ft*ht/dt
837 derit= kft + fpt*hdt + ft*udt*(one - hdt*hdt)
838 yt(i) = yt(i) - funt / derit
839 eminrlt(i)= dtt
840 sminrlt(i) = ft
841 ENDDO !iter
842 !TAGOLDT = 1
843 ELSEIF(ddet >=zero) THEN
844 !----------------------RELOAD
845 DO iter = 1, niter
846 ht = ht0 + yt(i) !longueur ressort
847 dt = sqrt(lt(i)*lt(i) + ht*ht) ! long fibre
848 dtt = dt - dt0
849 dtt = max(uvar(i,22),dtt)
850 dt = dtt + dt0
851 udt = one / dt
852 hdt = ht * udt !sinus(alpha)
853 epsnrlt(i) = uvar(i,21)*(dtt - uvar(i,22)) /(uvar(i,21)-uvar(i,22))
854 IF (kfunc(2) == 0) THEN
855 ft = kt0 * dtt
856 fpt = kt0 * hdt
857 ELSE
858 ft = yfac(2)*finter(kfunc(2),epsnrlt(i),npf,tf,fpt)
859 coefrl = (uvar(i,24) - uvar(i,26))/ uvar(i,24)
860 ft = uvar(i,26) + coefrl * ft
861
862 fpt = fpt * yfac(2)
863 fpt = fpt * uvar(i,21) * (uvar(i,24) - uvar(i,26))
864 . /(uvar(i,21)-uvar(i,22))/uvar(i,24)
865 kt = fpt
866 fpt = fpt * hdt
867 ENDIF
868 funt = kft*yt(i) + ft*ht/dt
869 derit= kft + fpt*hdt + ft*udt*(one - hdt*hdt)
870 yt(i) = yt(i) - funt / derit
871 ENDDO !iter
872 emaxrlt(i) = dtt
873 smaxrlt(i) = ft
874 IF(dtt>uvar(i,21))THEN
875 epsmaxt(i) = dtt
876 sigmaxt(i) = ft
877 !EMINRLT(I)= ZERO
878 !SMINRLT(I) = ZERO
879 ENDIF
880 !TAGOLDT = 2
881 ELSE !unload
882 !-----------------------UNLOAD
883 IF (uvar(i,23)>zero)THEN
884 DO iter = 1, niter
885 ht = ht0 + yt(i) !longueur ressort
886 dt = sqrt(lt(i)*lt(i) + ht*ht) ! long fibre
887 dtt = dt - dt0
888 dtt = max(zero,dtt)
889 dt = dtt + dt0
890 udt = one / dt
891 hdt = ht * udt !sinus(alpha)
892 epsnt(i) = epsi2*dtt/uvar(i,23)
893 IF (kfunc(2) == 0) THEN
894 ft = kt0 * dtt
895 fpt = kt0 * hdt
896 ELSE
897 ft = yfac(5)*finter(kfunc(5),epsnt(i),npf,tf,fpt)
898 fpt = fpt * yfac(5)
899 coeful = ft /sigi2
900 ft = coeful * uvar(i,28)
901 fpt = fpt * epsi2*uvar(i,28)/uvar(i,23)/sigi2
902 kt = fpt
903 fpt = fpt * hdt
904 ENDIF
905 eminrlt(i)= dtt
906 sminrlt(i) = ft
907 funt = kft*yt(i) + ft*ht/dt
908 derit= kft + fpt*hdt + ft*udt*(one - hdt*hdt)
909 yt(i) = yt(i) - funt / derit
910 ENDDO !iter
911 ELSE !(UVAR(I,23)>ZERO)
912 DO iter = 1, niter
913 ht = ht0 + yt(i) !longueur ressort
914 dt = sqrt(lt(i)*lt(i) + ht*ht) ! long fibre
915 udt = one / dt
916 hdt = ht * udt !sinus(alpha)
917 dtt = dt - dt0
918 ft = kkt * dtt
919 fpt = kkt * hdt! derivee % yT
920 kt = kkt
921 funt = kft*yt(i) + ft*ht/dt
922 derit= kft + fpt*hdt + ft*udt*(one - hdt*hdt)
923 yt(i) = yt(i) - funt / derit
924 ENDDO !iter
925 epsmaxt(i) = zero
926 sigmaxt(i) = zero
927 smaxrlt(i) = zero
928 emaxrlt(i)= zero
929 eminrlt(i)= zero
930 sminrlt(i) = zero
931 ENDIF !(UVAR(I,23)>ZERO)
932 !TAGOLDT = 3
933 ENDIF !LOAD RELOD OR UNLOAD
934 ELSE ! dtt < zero domaine de compression
935 DO iter = 1, niter
936 ht = ht0 + yt(i) !longueur ressort
937 dt = sqrt(lt(i)*lt(i) + ht*ht) ! long fibre
938 udt = one / dt
939 hdt = ht * udt !sinus(alpha)
940 dtt = dt - dt0
941 ft = kkt * dtt
942 fpt = kkt * hdt! derivee % yT
943 kt = kkt
944 funt = kft*yt(i) + ft*ht/dt
945 derit= kft + fpt*hdt + ft*udt*(one - hdt*hdt)
946 yt(i) = yt(i) - funt / derit
947 ENDDO !iter
948 !TAGOLDT = 5
949 epsmaxt(i) = zero
950 eminrlt(i) = zero
951 emaxrlt(i) = zero
952 sigmaxt(i) = zero
953 sminrlt(i) = zero
954 smaxrlt(i) = zero
955 ENDIF ! COMPRESSION OR TENSION FOR DIRECTION CHAINE
956 fn(i) = zero
957
958C ========================================================
959C--- Coupled model (traction fibre)
960C ========================================================
961 IF ((yc(i) + yt(i)) < zero) THEN !eps
962 epsmaxc(i) = uvar(i,17)
963 eminrlc(i) = uvar(i,18)
964 emaxrlc(i) = uvar(i,19)
965 sigmaxc(i) = uvar(i,20)
966 sminrlc(i) = uvar(i,25)
967 smaxrlc(i) = uvar(i,27)
968 epsmaxt(i) = uvar(i,21)
969 eminrlt(i) = uvar(i,22)
970 emaxrlt(i) = uvar(i,23)
971 sigmaxt(i) = uvar(i,24)
972 sminrlt(i) = uvar(i,26)
973 smaxrlt(i) = uvar(i,28)
974 tagoldc = nint(uvar(i,30))
975 tagoldt = nint(uvar(i,31))
976 hc = hc0 + uvar(i,7)
977 ht = ht0 + uvar(i,8)
978 dc = sqrt(lc(i)*lc(i) + hc*hc)
979 dt = sqrt(lt(i)*lt(i) + ht*ht)
980 ddec = dc - uvar(i,15)
981 ddet = dt - uvar(i,16)
982
983 y = half * (uvar(i,7) - uvar(i,8))
984 hc = hc0 + y
985 ht = ht0 - y
986 dc = sqrt(lc(i)*lc(i) + hc*hc)
987 dt = sqrt(lt(i)*lt(i) + ht*ht)
988 dcc0 = dc - dc0
989 dtt0 = dt - dt0
990 dcc1 = uvar(i,15) - dc0
991 dtt1 = uvar(i,16) - dt0
992 kf = kfc + kft
993
994 idxcnt = 1
995 idx1 = mod(idxcnt,2)
996 idx2 = mod(idxcnt+1,2)
997c---------------initialization
998 IF(dcc1 >=zero ) THEN
999c----------------------------------LOAD
1000 IF((ddec >=zero .AND.dcc0>uvar(i,17))
1001 . .OR.(ddec >=zero .AND.uvar(i,17)==uvar(i,18)) ) THEN
1002 tagc(idx1) = 1
1003 tagc(idx2) = 1
1004 ifuncc = kfunc(1)
1005 facc = yfac(1)
1006 eminc = zero
1007 sminc = zero
1008 face = one
1009 facf = one
1010c----------------------------------RELOAD
1011 ELSEIF(ddec >=zero)THEN !RELOAD
1012 tagc(idx1) = 2
1013 tagc(idx2) = 2
1014 ifuncc= kfunc(1)
1015 facc = yfac(1)
1016 eminc = uvar(i,18)
1017 sminc = uvar(i,25)
1018 face = uvar(i,17)/(uvar(i,17)-uvar(i,18))
1019 facf = (uvar(i,20) - uvar(i,25))/uvar(i,20)
1020 IF (face <= zero) THEN
1021c !print*,'ERROR1: CYCLE,ELEM,FACE=',NCYCLE,NGL(I),FACE
1022c !print*,'a EMAXC,EMINRLC=',UVAR(I,17),UVAR(I,18)
1023 face = abs(face)
1024 ENDIF
1025 IF (facf <= zero) THEN
1026c !print*,'ERROR1: CYCLE,ELEM,FACF=',NCYCLE,NGL(I),FACF
1027c !print*,'a SMAXC,SMINRLC=',UVAR(I,20),UVAR(I,25)
1028 facf = abs(facf)
1029 ENDIF
1030 ELSE !UNLOAd
1031c----------------------------------UNLOAD
1032 IF(uvar(i,19) >zero)THEN
1033 tagc(idx1) = 3
1034 tagc(idx2) = 3
1035 ifuncc= kfunc(4)
1036 facc = yfac(4)
1037 eminc = zero
1038 sminc = zero
1039 face = epsi1 /uvar(i,19)
1040 facf = uvar(i,27)/sigi1
1041 else!UNLOAd UVAR(I,19) <=ZERO
1042 tagc(idx1) = 5
1043 tagc(idx2) = 5
1044 endif!UNLOAd UVAR(I,19) <=ZERO
1045 endif!UNLOAd RELOAD LOAD
1046 ELSE ! DOMAINE POSITIF
1047 tagc(idx1) = 5
1048 tagc(idx2) = 5
1049 ENDIF ! (DCC1 >=ZERO )
1050c ---------------direction trame ---------------------
1051 IF(dtt1 >=zero ) THEN
1052 IF((ddet >=zero .AND.dtt0>uvar(i,21))
1053 . .OR.(ddet >=zero .AND.uvar(i,21)==uvar(i,22)) ) THEN
1054c----------------------------------LOAD
1055 tagt(idx1) = 1
1056 tagt(idx2) = 1
1057 ifunct = kfunc(2)
1058 fact = yfac(2)
1059 emint = zero
1060 smint = zero
1061 fate = one
1062 fatf = one
1063 ELSEIF(ddet >=zero)THEN
1064c--------------------------------RELOAD
1065 tagt(idx1) = 2
1066 tagt(idx2) = 2
1067 ifunct= kfunc(2)
1068 fact = yfac(2)
1069 emint = uvar(i,22)
1070 smint = uvar(i,26)
1071 fate = uvar(i,21) /(uvar(i,21)-uvar(i,22))
1072 fatf = (uvar(i,24) - uvar(i,26))/uvar(i,24)
1073 IF (fate <= zero) THEN
1074c !print*,'ERROR1: CYCLE,ELEM,FATE=',NCYCLE,NGL(I),FATE
1075c !print*,'a EMAXT,EMINRLT=',UVAR(I,21),UVAR(I,22)
1076 fate = abs(fate)
1077 ENDIF
1078 IF (fatf <= zero) THEN
1079c !print*,'ERROR1: CYCLE,ELEM,FATF=',NCYCLE,NGL(I),FATF
1080c !print*,'a SMAXT,SMINRLT=',UVAR(I,24),UVAR(I,26)
1081 fatf = abs(fatf)
1082 ENDIF
1083 ELSE !
1084c--------------------------------UNLOAD
1085 IF(uvar(i,23) >zero)THEN
1086 tagt(idx1) = 3
1087 tagt(idx2) = 3
1088 ifunct = kfunc(5)
1089 fact = yfac(5)
1090 emint = zero
1091 smint = zero
1092 fate = epsi2 /uvar(i,23)
1093 fatf = uvar(i,28)/sigi2
1094 ELSE !unload (UVAR(I,23) >ZERO)
1095 tagt(idx1) = 5
1096 tagt(idx2) = 5
1097 ENDIF
1098 ENDIF
1099 ELSE ! DOMAINE POSITIF
1100 ! compression fibre
1101 tagt(idx1) = 5
1102 tagt(idx2) = 5
1103 ENDIF
1104c-------------------------------------------
1105 ierr = 1
1106 DO WHILE (ierr == 1)
1107c-------------------------------------------
1108 DO iter = 1, niter
1109c--- chaine
1110 hc = hc0 + y
1111 ht = ht0 - y
1112 dc = sqrt(lc(i)*lc(i) + hc*hc)
1113 dt = sqrt(lt(i)*lt(i) + ht*ht)
1114 dcc = dc - dc0
1115 dtt = dt - dt0
1116c
1117 SELECT CASE (tagc(idx1))
1118 CASE(1) ! tension load or reload
1119 dcc = max(uvar(i,18),dcc)
1120 dc = dcc + dc0
1121 udc = one / dc
1122 hdc = hc * udc
1123 epsnrlc(i) = (dcc -eminc) *face
1124 IF (ifuncc == 0) THEN
1125 fc = kc0 * dcc
1126 fpc = kc0 * hdc
1127 ELSE
1128 fc = facc * finter(ifuncc,epsnrlc(i),npf,tf,fpc)
1129 fc = sminc + fc *facf
1130 fpc = fpc * facc
1131 fpc = fpc *facf * face
1132 kc = fpc
1133 fpc = fpc * hdc
1134 ENDIF
1135 IF ( iter==niter) then! DCC > UVAR(I,19) .AND.
1136 epsmaxc(i) = dcc
1137 sigmaxc(i) = fc
1138 emaxrlc(i) = dcc
1139 smaxrlc(i) = fc
1140 eminrlc(i) = dcc
1141 sminrlc(i) = fc
1142 ENDIF
1143 !IF (DCC > UVAR(I,17) .AND. ITER==NITER) THEN
1144 ! EPSMAXC(I) = DCC
1145 ! SIGMAXC(I) = FC
1146 !EMINRLC(I) = ZERO
1147 !SMINRLC(I) = ZERO
1148 !ENDIF
1149 CASE(2) ! tension load or reload
1150 dcc = max(uvar(i,18),dcc)
1151 dc = dcc + dc0
1152 udc = one / dc
1153 hdc = hc * udc
1154 epsnrlc(i) = (dcc -eminc) *face
1155 IF (ifuncc == 0) THEN
1156 fc = kc0 * dcc
1157 fpc = kc0 * hdc
1158 ELSE
1159 fc = facc * finter(ifuncc,epsnrlc(i),npf,tf,fpc)
1160 fc = sminc + fc *facf
1161 fpc = fpc * facc
1162 fpc = fpc *facf * face
1163 kc = fpc
1164 fpc = fpc * hdc
1165 ENDIF
1166 IF ( iter==niter) THEN ! DCC > UVAR(I,19) .AND.
1167 emaxrlc(i) = dcc
1168 smaxrlc(i) = fc
1169 ENDIF
1170 IF (dcc > uvar(i,17) .AND. iter==niter) THEN
1171 epsmaxc(i) = dcc
1172 sigmaxc(i) = fc
1173 !EMINRLC(I) = ZERO
1174 !SMINRLC(I) = ZERO
1175 ENDIF
1176 CASE(3) ! tension unload
1177 dcc = max(zero,dcc)
1178 dc = dcc + dc0
1179 udc = one / dc
1180 hdc = hc * udc
1181 epsnrlc(i) = (dcc -eminc) *face
1182 IF (ifuncc == 0) THEN
1183 fc = kc0 * dcc
1184 fpc = kc0 * hdc
1185 ELSE
1186 fc = facc * finter(ifuncc,epsnrlc(i),npf,tf,fpc)
1187 fc = sminc + fc *facf
1188 fpc = fpc * facc
1189 fpc = fpc *facf * face
1190 kc = fpc
1191 fpc = fpc * hdc
1192 ENDIF
1193 IF (dcc < uvar(i,18)) THEN
1194 eminrlc(i) = dcc
1195 sminrlc(i) = fc
1196 ENDIF
1197 CASE(5) ! compression
1198 udc= one / dc
1199 hdc= hc * udc
1200 fc = kkc * dcc
1201 fpc = kkc * hdc! derivee % yc
1202 kc = kkc
1203 epsmaxc(i) = zero
1204 eminrlc(i) = zero
1205 emaxrlc(i) = zero
1206 sigmaxc(i) = zero
1207 sminrlc(i) = zero
1208 smaxrlc(i) = zero
1209 END SELECT
1210c
1211c -----------------------direction trame
1212 SELECT CASE (tagt(idx1))
1213 CASE(1) ! tension load or reload
1214 dtt = max(uvar(i,22),dtt)
1215 dt = dtt + dt0
1216 udt = one / dt
1217 hdt = ht * udt !sinus(alpha)
1218 epsnrlt(i) = (dtt - emint)*fate
1219 IF (ifunct == 0) THEN
1220 ft = kt0 * dtt
1221 fpt = kt0 * hdt
1222 ELSE
1223 ft = fact*finter(ifunct,epsnrlt(i),npf,tf,fpt)
1224 ft = smint + ft * fatf
1225 fpt = fpt * fact
1226 fpt = fpt *fatf * fate
1227 kt = fpt
1228 fpt = fpt * hdt
1229 ENDIF
1230 IF (iter==niter) THEN !DTT > UVAR(I,23) .AND.
1231 emaxrlt(i) = dtt
1232 smaxrlt(i) = ft
1233 epsmaxt(i) = dtt
1234 sigmaxt(i) = ft
1235 eminrlt(i) = dtt
1236 sminrlt(i) = ft
1237 ENDIF
1238 !IF(DTT > UVAR(I,21).AND.ITER==NITER)THEN
1239 ! EPSMAXT(I) = DTT
1240 ! SIGMAXT(I) = FT
1241 ! EMINRLT(I) = ZERO
1242 ! SMINRLT(I) = ZERO
1243 !ENDIF
1244
1245 CASE(2) ! tension load or reload
1246 dtt = max(uvar(i,22),dtt)
1247 dt = dtt + dt0
1248 udt = one / dt
1249 hdt = ht * udt !sinus(alpha)
1250 epsnrlt(i) = (dtt - emint)*fate
1251 IF (ifunct == 0) THEN
1252 ft = kt0 * dtt
1253 fpt = kt0 * hdt
1254 ELSE
1255 ft = fact*finter(ifunct,epsnrlt(i),npf,tf,fpt)
1256 ft = smint + ft * fatf
1257 fpt = fpt * fact
1258 fpt = fpt *fatf * fate
1259 kt = fpt
1260 fpt = fpt * hdt
1261 ENDIF
1262 IF (iter==niter) THEN !DTT > UVAR(I,23) .AND.
1263 emaxrlt(i) = dtt
1264 smaxrlt(i) = ft
1265 ENDIF
1266 IF(dtt > uvar(i,21).AND.iter==niter)THEN
1267 epsmaxt(i) = dtt
1268 sigmaxt(i) = ft
1269 !EMINRLT(I) = ZERO
1270 !SMINRLT(I) = ZERO
1271 ENDIF
1272
1273 CASE(3) ! tension unload
1274 dtt = max (zero,dtt)
1275 dt = dtt + dt0
1276 udt = one / dt
1277 hdt = ht * udt !sinus(alpha)
1278 epsnrlt(i) = (dtt - emint)*fate
1279 IF (ifunct == 0) THEN
1280 ft = kt0 * dtt
1281 fpt = kt0 * hdt
1282 ELSE
1283 ft = fact*finter(ifunct,epsnrlt(i),npf,tf,fpt)
1284 ft = smint + ft * fatf
1285 fpt = fpt * fact
1286 fpt = fpt *fatf * fate
1287 kt = fpt
1288 fpt = fpt * hdt
1289 ENDIF
1290 IF (dtt < uvar(i,22)) THEN
1291 eminrlt(i) = dtt
1292 sminrlt(i) = ft
1293 ENDIF
1294 CASE(5) ! compression
1295 udt= one / dt
1296 hdt= ht * udt
1297 ft = kkt * dtt
1298 fpt = kkt * hdt! derivee % yT
1299 kt = kkt
1300 epsmaxt(i) = zero
1301 eminrlt(i) = zero
1302 emaxrlt(i) = zero
1303 sigmaxt(i) = zero
1304 sminrlt(i) = zero
1305 smaxrlt(i) = zero
1306 END SELECT
1307c----
1308 func = kf*y + fc * hdc - ft * hdt
1309 deric = kf + fpc*hdc + fc*udc*(one - hdc*hdc)
1310 . + fpt*hdt + ft*udt*(one - hdt*hdt)
1311 y = y - func / deric
1312
1313 IF (y > 0) THEN
1314 y = min(y, ht0)
1315 ELSE
1316 y = max(y,-hc0)
1317 ENDIF
1318
1319 ENDDO ! ITER = 1, NITER
1320 yc(i) = y
1321 yt(i) =-y
1322 fn(i) = fc*hc/dc + ft*ht/dt
1323c----------------------------------
1324 hc = hc0 + yc(i)
1325 ht = ht0 + yt(i)
1326 dc = sqrt(lc(i)*lc(i) + hc*hc)
1327 dt = sqrt(lt(i)*lt(i) + ht*ht)
1328C------------------------------------------------------------------
1329c POST testing
1330C------------------------------------------------------------------
1331 ierr = 0
1332 ierrt = 0
1333 ierrc = 0
1334 IF (tagc(idx1) == 3 .and. ((eminrlc(i) > uvar(i,18) .and. tagoldc == 3)
1335 . .OR. eminrlc(i) >uvar(i,17) .OR. eminrlc(i) >uvar(i,19))) THEN
1336 ierr = 1
1337 ierrc= 1
1338 eminrlc(i)= uvar(i,18)
1339 sminrlc(i) = uvar(i,25)
1340 dc = uvar(i,15)
1341 y = half * (uvar(i,7) - uvar(i,8))
1342 ifuncc = kfunc(1)
1343 facc = yfac(1)
1344 IF (uvar(i,18)==zero .OR. uvar(i,17)==uvar(i,18))THEN
1345 tagc(idx2) = 1
1346 eminc = zero
1347 sminc = zero
1348 face = one
1349 facf = one
1350 ELSE
1351 tagc(idx2) = 2
1352 eminc = uvar(i,18)
1353 sminc = uvar(i,25)
1354 face = uvar(i,17) /(uvar(i,17)-uvar(i,18))
1355 facf = (uvar(i,20) - uvar(i,25))/uvar(i,20)
1356 IF (face <= zero) THEN
1357c !print*,'ERROR: CYCLE,ELEM,FACE=',NCYCLE,NGL(I),FACE
1358c !print*,' EMAXC,EMINRLC=',UVAR(I,17),UVAR(I,18)
1359 face = abs(face)
1360 ENDIF
1361 IF (facf <= zero) THEN
1362c !print*,'ERROR: CYCLE,ELEM,FACF=',NCYCLE,NGL(I),FACF
1363c !print*,' SMAXC,SMINRLC=',UVAR(I,20),UVAR(I,25)
1364 facf = abs(facf)
1365 ENDIF
1366 ENDIF
1367 ELSEIF ((tagc(idx1) == 1.or.tagc(idx1) == 2) .and. (dcc < uvar(i,19).and. tagoldc /= 3 )) THEN
1368 ierr = 1
1369 ierrc= 1
1370 epsmaxc(i) = uvar(i,17)
1371 eminrlc(i) = uvar(i,18)
1372 emaxrlc(i) = uvar(i,19)
1373 sigmaxc(i) = uvar(i,20)
1374 sminrlc(i) = uvar(i,25)
1375 smaxrlc(i) = uvar(i,27)
1376 dc = uvar(i,15)
1377 y = half * (uvar(i,7) - uvar(i,8))
1378 tagc(idx2) = 3
1379 ifuncc= kfunc(4)
1380 facc = yfac(4)
1381 eminc = zero
1382 sminc = zero
1383 face = epsi1 /uvar(i,19)
1384 facf = uvar(i,27)/sigi1
1385 ENDIF !IF (TAGC(IDX1)
1386C---------------------------TRAME------------------------------------
1387 IF (tagt(idx1) == 3 .and. ((eminrlt(i) > uvar(i,22) .and.tagoldt == 3)
1388 . .OR. eminrlt(i) >uvar(i,21).OR. eminrlt(i) >uvar(i,23))) THEN
1389 ierr = 1
1390 ierrt= 1
1391 eminrlt(i)= uvar(i,22)
1392 sminrlt(i) = uvar(i,26)
1393 dt = uvar(i,16)
1394 y = half * (uvar(i,7) - uvar(i,8))
1395 ifunct = kfunc(2)
1396 fact = yfac(2)
1397 IF(uvar(i,22)==zero.OR.uvar(i,21)==uvar(i,22))THEN
1398 tagt(idx2) = 1
1399 emint = zero
1400 smint = zero
1401 fate = one
1402 fatf = one
1403 ELSE
1404 tagt(idx2) = 2
1405 emint = uvar(i,22)
1406 smint = uvar(i,26)
1407 fate = uvar(i,21)/(uvar(i,21)-uvar(i,22))
1408 fatf = (uvar(i,24) - uvar(i,26))/uvar(i,24)
1409 IF (fate <= zero) THEN
1410c !print*,'ERROR: CYCLE,ELEM,FATE=',NCYCLE,NGL(I),FATE
1411c !print*,' EMAXT,EMINRLT=',UVAR(I,21),UVAR(I,22)
1412 fate = abs(fate)
1413 ENDIF
1414 IF (fatf <= zero) THEN
1415c !print*,'ERROR: CYCLE,ELEM,FATF=',NCYCLE,NGL(I),FATF
1416c !print*,' SMAXT,SMINRLT=',UVAR(I,24),UVAR(I,26)
1417 fatf = abs(fatf)
1418 ENDIF
1419 ENDIF
1420
1421 ELSEIF ((tagt(idx1) == 1.or.tagt(idx1) == 2).and.(dtt < uvar(i,23).and. tagoldt /= 3 ))THEN
1422 ierr = 1
1423 ierrt= 1
1424 epsmaxt(i) = uvar(i,21)
1425 eminrlt(i) = uvar(i,22)
1426 emaxrlt(i) = uvar(i,23)
1427 sigmaxt(i) = uvar(i,24)
1428 sminrlt(i) = uvar(i,26)
1429 smaxrlt(i) = uvar(i,28)
1430 dt = uvar(i,16)
1431 y = half * (uvar(i,7) - uvar(i,8))
1432 tagt(idx2) = 3
1433 ifunct= kfunc(5)
1434 fact = yfac(5)
1435 emint = zero
1436 smint = zero
1437 fate = epsi2 /uvar(i,23)
1438 fatf = uvar(i,28)/sigi2
1439 ENDIF ! IF (TAGT(IDX1)
1440 idxcnt = idxcnt + 1
1441
1442 idx1 = mod(idxcnt,2)
1443
1444 idx2 = mod(idxcnt+1,2)
1445
1446 IF (idxcnt > 3) THEN
1447 !!print*, 'IDXCNT > 3'
1448 IF (ierrc == 1)THEN
1449 eminrlc(i) = uvar(i,18)
1450 sminrlc(i) = uvar(i,25)
1451 emaxrlc(i) = uvar(i,19)
1452 smaxrlc(i) = uvar(i,27)
1453 epsmaxc(i) = uvar(i,17)
1454 sigmaxc(i) = uvar(i,20)
1455 IF (dcc > uvar(i,19)) THEN
1456 emaxrlc(i) = dcc
1457 smaxrlc(i) = fc
1458 ENDIF
1459 IF (dcc > uvar(i,17)) THEN
1460 epsmaxc(i) = dcc
1461 sigmaxc(i) = fc
1462 eminrlc(i) = zero
1463 sminrlc(i) = zero
1464 ELSEIF (dcc < uvar(i,18)) THEN
1465 eminrlc(i) = dcc
1466 sminrlc(i) = fc
1467 ENDIF
1468 ENDIF !IF (IERRC == 1)
1469 IF (ierrt==1) THEN
1470 eminrlt(i) = uvar(i,22)
1471 sminrlt(i) = uvar(i,26)
1472 epsmaxt(i) = uvar(i,21)
1473 sigmaxt(i) = uvar(i,24)
1474 epsmaxt(i) = uvar(i,21)
1475 emaxrlt(i) = dtt
1476 smaxrlt(i) = ft
1477 IF (dtt > uvar(i,23)) THEN
1478 emaxrlt(i) = dtt
1479 smaxrlt(i) = ft
1480 ENDIF
1481 IF (dtt > uvar(i,21)) THEN
1482 epsmaxt(i) = dtt
1483 sigmaxt(i) = ft
1484 eminrlt(i) = zero
1485 sminrlt(i) = zero
1486 ELSEIF (dtt < uvar(i,22)) THEN
1487 eminrlt(i) = dtt
1488 sminrlt(i) = ft
1489 ENDIF
1490 ENDIF !IF (IERRT==1)
1491 EXIT
1492 ENDIF
1493
1494 ENDDO ! WHILE (IERR == 1)
1495 uvar(i,30) = tagc(idx2)
1496 uvar(i,31) = tagt(idx2)
1497c
1498 ENDIF ! unload option
1499C=======================================================================
1500c SHEAR
1501C=======================================================================
1502 trace = exp(epsxx(i) + epsyy(i))! =exp(tr)
1503 ec2 = max(trace / (ec(i) + one), em6)
1504 et2 = max(trace / (et(i) + one), em6)
1505c
1506 rfac = one / ec2
1507 rfat = one / et2
1508C--- contraintes membrane
1509c
1510 sigc = fc * lc(i) / dc
1511 sigt = ft * lt(i) / dt
1512 signxx(i) = sigc * rfac
1513 signyy(i) = sigt * rfat
1514c
1515 uvar(i,1) = sigc
1516 uvar(i,2) = sigt
1517 uvar(i,15) = dc
1518 uvar(i,16) = dt
1519C------------------------------------------------------------------
1520C--- contraintes cisaillement
1521C------------------------------------------------------------------
1522 sig0 = uvar(i,10)
1523 tan_phi(i)= depsxy(i)
1524 uvar(i,6)= depsxy(i)
1525 tan2 = tan_phi(i)*tan_phi(i)
1526 phio = uvar(i,29)
1527 phi = atan(tan_phi(i))*hundred80/pi
1528 dphi = phi - phio
1529 uvar(i,29) = phi
1530
1531 phimin(i) = uvar(i,32)
1532 phimax(i) = uvar(i,33)
1533 phirlmax(i) = uvar(i,34)
1534 sxymax(i) = uvar(i,35)
1535 sxymin(i) = uvar(i,36)
1536 sxymaxrl(i) = uvar(i,37)
1537
1538 IF (kfunc(3) == 0) THEN ! analytic shear without UNLOAD
1539 IF (tan_phi(i) > tan_lock) THEN
1540 signxy(i) = gt*tan_phi(i) + gb - sig0
1541 kg = gt*tan2
1542 ELSEIF (tan_phi(i) < -tan_lock) THEN
1543 signxy(i) = gt*tan_phi(i) - gb - sig0
1544 kg = gt*tan2
1545 ELSE
1546 signxy(i) = g0*tan_phi(i) - sig0
1547 kg = g0*tan2
1548 ENDIF
1549c-----
1550 ELSE ! KFUNC(3) > 0 => load/UNLOAD
1551c-----
1552 IF (phi >=zero)THEN
1553 IF ((dphi >=zero .AND. uvar(i,32) ==zero).OR.(dphi >=zero .AND.
1554 . phi >uvar(i,33)).OR.(dphi >=zero .AND. uvar(i,33)==uvar(i,32)) .OR.
1555 . (dphi >=zero .AND. uvar(i,35)==zero)) THEN
1556 !----------------------LOAD
1557 signxy(i) = yfac(3)*finter(kfunc(3),phi,npf,tf,gxy)
1558 kg = gxy*yfac(3)*tan2
1559 phimax(i) = phi
1560 phirlmax(i) = phi
1561 sxymax(i) = signxy(i)
1562 sxymaxrl(i) = signxy(i)
1563 ELSEIF(dphi>=zero)THEN
1564 !----------------------RELOAD
1565 phin = uvar(i,33)*(phi - uvar(i,32))/(uvar(i,33) - uvar(i,32))
1566 signxy(i) = yfac(3)*finter(kfunc(3),phin,npf,tf,gxy)
1567 signxy(i) = uvar(i,36) + signxy(i) *(uvar(i,35) - uvar(i,36))/uvar(i,35)
1568 kg = gxy*yfac(3)*tan2
1569 phirlmax(i) = phi
1570 sxymaxrl(i) = signxy(i)
1571 IF (phi >uvar(i,33))THEN
1572 phimax(i) = phi
1573 sxymax(i) = signxy(i)
1574 ENDIF
1575 ELSE
1576 !----------------------UNLOAD
1577 phin = phi * phii/uvar(i,34)
1578 signxy(i) = yfac(6)*finter(kfunc(6),phin,npf,tf,gxy)
1579 signxy(i) = signxy(i)*uvar(i,37) / sxyi
1580 kg = gxy*yfac(6)*tan2
1581 phimin(i) = phi
1582 sxymin(i) = signxy(i)
1583 endif!DPHI LOAD - RELOAD - UNLOAD
1584 ELSE ! PHI <ZERO
1585 IF((dphi <=zero .AND. uvar(i,32) ==zero).OR.(dphi <=zero .AND. phi <uvar(i,33))
1586 . .OR.(dphi <=zero .AND. uvar(i,33)==uvar(i,32)) ) THEN
1587 !----------------------LOAD
1588 signxy(i) = yfac(3)*finter(kfunc(3),phi,npf,tf,gxy)
1589 kg = gxy*yfac(3)*tan2
1590 phimax(i) = phi
1591 phirlmax(i) = phi
1592 sxymax(i) = signxy(i)
1593 sxymaxrl(i) = signxy(i)
1594 ELSEIF(dphi<=zero)THEN
1595 !----------------------RELOAD
1596 phin = uvar(i,33)*(phi - uvar(i,32))/(uvar(i,33) - uvar(i,32))
1597 signxy(i) = yfac(3)*finter(kfunc(3),phin,npf,tf,gxy)
1598 signxy(i) = uvar(i,36) + signxy(i) *(uvar(i,35) - uvar(i,36))/uvar(i,35)
1599 kg = gxy*yfac(3)*tan2
1600 phirlmax(i) = phi
1601 sxymaxrl(i) = signxy(i)
1602 IF (phi <uvar(i,33))THEN
1603 phimax(i) = phi
1604 sxymax(i) = signxy(i)
1605 ENDIF
1606 ELSE
1607 !----------------------UNLOAD
1608 phin = phi * phii2/uvar(i,34)
1609 signxy(i) = yfac(6)*finter(kfunc(6),phin,npf,tf,gxy)
1610 signxy(i) = signxy(i)*uvar(i,37) / sxyi2
1611 kg = gxy*yfac(6)*tan2
1612 phimin(i) = phi
1613 sxymin(i) = signxy(i)
1614 endif!DPHI LOAD - RELOAD - UNLOAD
1615 ENDIF
1616c
1617 ENDIF ! KFUNC(3) > 0
1618c---------------------------------------------------------------------
1619C--- contraintes visc fibres
1620c
1621 dtinv = timestep(i)/max(timestep(i)*timestep(i),em20)
1622 damp = sqrt(rho0(i)*area(i)*thkly(i)*half)
1623 IF (kc < zero) THEN
1624 kc = one
1625 ENDIF
1626 IF (kt < zero) THEN
1627 kt = one
1628 ENDIF
1629 v1 = visce*damp*sqrt(kc)
1630 v2 = visce*damp*sqrt(kt)
1631 sigvxx(i) = dtinv*(depsxx(i))*v1
1632 sigvyy(i) = dtinv*(depsyy(i))*v2
1633C--- contraintes frottement cisaillement
1634 IF (fn(i) > zero) THEN
1635 tfrot = two_third*viscg*fn(i)*(hc0+ht0)/(lc(i)+lt(i))
1636 dtang(i) = depsxy(i) - tan_phi(i)
1637 sigg = tfold(i) + gfrot*dtang(i)
1638 IF (abs(sigg) > tfrot) THEN
1639 sigvxy(i) = sign(tfrot,sigg)
1640 ELSE
1641 sigvxy(i) = sigg
1642 ENDIF
1643 ENDIF
1644C---
1645 sinn = tan_phi(i) / sqrt(one + tan2)
1646 kmax = max(kc0*rfac,kt0*rfat)*(one+sinn) + kg
1647 lmin = min(dc/dc0,dt/dt0)*uvar(i,14)
1648 soundsp(i) = sqrt(kmax/(rho0(i)))*aldt(i)/lmin
1649 viscmax(i) = max(v1,v2)
1650 etse(i) = one
1651 ENDDO
1652C---
1653#include "vectorize.inc"
1654 DO j=1,nindx
1655 i=index(j)
1656 uvar(i,3) = signxy(i) + sigvxy(i)
1657 tan_phi(i) = depsxy(i)
1658 uvar(i,6) = depsxy(i)
1659 uvar(i,7) = yc(i)
1660 uvar(i,8) = yt(i)
1661 uvar(i,9) = sigvxy(i)
1662c
1663 uvar(i,17) = epsmaxc(i)
1664 uvar(i,18) = eminrlc(i)
1665 uvar(i,19) = emaxrlc(i)
1666 uvar(i,20) = sigmaxc(i)
1667 uvar(i,21) = epsmaxt(i)
1668 uvar(i,22) = eminrlt(i)
1669 uvar(i,23) = emaxrlt(i)
1670 uvar(i,24) = sigmaxt(i)
1671 uvar(i,25) = sminrlc(i)
1672 uvar(i,26) = sminrlt(i)
1673 uvar(i,27) = smaxrlc(i)
1674 uvar(i,28) = smaxrlt(i)
1675 uvar(i,32) = phimin(i)
1676 uvar(i,33) = phimax(i)
1677 uvar(i,34) = phirlmax(i)
1678 uvar(i,35) = sxymax(i)
1679 uvar(i,36) = sxymin(i)
1680 uvar(i,37) = sxymaxrl(i)
1681c
1682 signyz(i) = sigoyz(i) + gsh *depsyz(i)
1683 signzx(i) = sigozx(i) + gsh *depszx(i)
1684 ENDDO
1685
1686 ENDIF ! UNLOAD
1687C-----------------------------------------------------------
1688C REF-STATE ZEROSTRESS OPTION
1689C-----------------------------------------------------------
1690 IF (zerostress /= zero)THEN
1691 IF (tt <= tstart) THEN
1692#include "vectorize.inc"
1693 DO j=1,nindx
1694 i=index(j)
1695 uvar(i,11) = signxx(i)
1696 uvar(i,12) = signyy(i)
1697 uvar(i,13) = signxy(i)
1698 signxx(i) = zero
1699 signyy(i) = zero
1700 signxy(i) = zero
1701 ENDDO
1702 ELSE
1703#include "vectorize.inc"
1704 DO j=1,nindx
1705 i=index(j)
1706 dsig = signxx(i) - sigoxx(i) - uvar(i,11)
1707 IF((uvar(i,11) > zero).AND.(dsig < zero))THEN
1708 uvar(i,11) = max(zero,uvar(i,11)+zerostress*dsig)
1709 ELSEIF((uvar(i,11) < zero).AND.(dsig > zero))THEN
1710 uvar(i,11) = min(zero,uvar(i,11)+zerostress*dsig)
1711 ENDIF
1712 dsig = signyy(i) - sigoyy(i) - uvar(i,12)
1713 IF((uvar(i,12) > zero).AND.(dsig < zero))THEN
1714 uvar(i,12) = max(zero,uvar(i,12)+zerostress*dsig)
1715 ELSEIF((uvar(i,12) < zero).AND.(dsig > zero))THEN
1716 uvar(i,12) = min(zero,uvar(i,12)+zerostress*dsig)
1717 ENDIF
1718 dsig = signxy(i) - sigoxy(i) - uvar(i,13)
1719 IF((uvar(i,13) > zero).AND.(dsig < zero))THEN
1720 uvar(i,13) = max(zero,uvar(i,13)+zerostress*dsig)
1721 ELSEIF((uvar(i,13) < zero).AND.(dsig > zero))THEN
1722 uvar(i,13) = min(zero,uvar(i,13)+zerostress*dsig)
1723 ENDIF
1724 signxx(i) = signxx(i) - uvar(i,11)
1725 signyy(i) = signyy(i) - uvar(i,12)
1726 signxy(i) = signxy(i) - uvar(i,13)
1727 ENDDO
1728 ENDIF
1729 ENDIF
1730c
1731 DO j=1,numoff
1732 i = indxoff(j)
1733 kmax = max(kc0,kt0)*two
1734 soundsp(i) = sqrt(kmax/(rho0(i)))
1735 ENDDO
1736C-----------------------------------------------------------
1737C DESACTIVATION SI SURFACE < SURFACE MIN
1738C-----------------------------------------------------------
1739 IF (areamin > zero) THEN
1740#include "vectorize.inc"
1741 DO j=1,nindx
1742 i=index(j)
1743 aa = (one+ec(i)+et(i) - areamin) * dareamin
1744 aa = min(max(aa,zero),one)
1745 signxx(i) = signxx(i)*aa
1746 signyy(i) = signyy(i)*aa
1747 signxy(i) = signxy(i)*aa
1748 signyz(i) = signyz(i)*aa
1749 signzx(i) = signzx(i)*aa
1750 ENDDO
1751 ENDIF
1752 !IF(ncycle== 2)!stop
1753C-----------------------------------------------------------
1754C DEACTIVATION FOR VERY WARPED ELEMENT (QBAT)
1755C-----------------------------------------------------------
1756 IF (npt0 ==1) THEN
1757#include "vectorize.inc"
1758 DO j=1,nindx
1759 i=index(j)
1760 IF (uvar(i,40)==one) cycle
1761 aa =uvar(i,40)
1762 signxx(i) = signxx(i)*aa
1763 signyy(i) = signyy(i)*aa
1764 signxy(i) = signxy(i)*aa
1765 IF (uvar(i,40)<em02) THEN
1766 kmax = max(kc0,kt0)*two
1767 soundsp(i) = sqrt(kmax/(rho0(i)))
1768 END IF
1769 ENDDO
1770 ENDIF
1771C-----------
1772 RETURN
1773 END SUBROUTINE sigeps58c
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
interface_table_mod
Definition table_mod.F:238
table_mod
Definition table_mod.F:112
sigeps58c
subroutine sigeps58c(nel0, nuparam, nuvar, mfunc, kfunc, npf, npt0, ipt, nsensor, tf, time, timestep, uparam, rho0, area, eint, thkly, niparam, iparam, depsxx, depsyy, depsxy, depsyz, depszx, epsxx, epsyy, epsxy, epsyz, epszx, sigoxx, sigoyy, sigoxy, sigoyz, sigozx, signxx, signyy, signxy, signyz, signzx, sigvxx, sigvyy, sigvxy, sigvyz, sigvzx, soundsp, viscmax, thk, pla, uvar, off, ngl, pm, mat, etse, shf, yld, tan_phi, aldt, sensor_tab, ismstr, table, offgg)
Definition sigeps58c.F:48
table_mod::ttable
Definition table_mod.F:125