OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
cm58_refsta.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!|| cm58_refsta ../starter/source/materials/mat/mat058/cm58_refsta.F
25!||--- called by ------------------------------------------------------
26!|| cmlawi ../starter/source/elements/shell/coque/cepsini.F
27!||--- calls -----------------------------------------------------
28!|| finter ../starter/source/tools/curve/finter.F
29!||====================================================================
30 SUBROUTINE cm58_refsta(NEL ,NUPARAM,NIPARAM,UPARAM ,IPARAM ,
31 2 FOR ,EINT ,GSTR ,
32 3 THK ,DIR1 ,DIR2 ,NUVAR ,UVAR ,
33 4 NBFUNC ,IFUNC ,NPF ,TF ,AREA ,
34 5 EXX ,EYY ,EXY )
35C-----------------------------------------------
36C I m p l i c i t T y p e s
37C-----------------------------------------------
38#include "implicit_f.inc"
39C-----------------------------------------------
40C G l o b a l P a r a m e t e r s
41C-----------------------------------------------
42#include "mvsiz_p.inc"
43C-----------------------------------------------
44C D u m m y A r g u m e n t s
45C-----------------------------------------------
46 INTEGER :: NEL,NUPARAM,NIPARAM,NUVAR,NBFUNC
47 INTEGER ,INTENT(IN) :: IFUNC(NBFUNC)
48 INTEGER ,INTENT(IN) :: IPARAM(NIPARAM)
49 my_real ,INTENT(IN) :: UPARAM(NUPARAM)
50 my_real :: for(nel,5),thk(*),uvar(nel,nuvar),area(nel),
51 . eint(nel,2),gstr(nel,8),dir1(nel,*),dir2(nel,*),
52 . exx(nel),eyy(nel),exy(nel)
53C-----------------------------------------------
54 INTEGER NPF(*)
55 my_real FINTER ,TF(*)
56 EXTERNAL FINTER
57C-----------------------------------------------
58C L o c a l V a r i a b l e s
59C-----------------------------------------------
60 INTEGER I,J,MX,NC,NT,ITER,NITER,FLAGUL
61 my_real
62 . KC,KT,KFC,KFT,G0,GT,GB,TAN_LOCK,VISCE,VISCG,LC0,LT0,
63 . KBC,KBT,DC0,DT0,HC0,HT0,ETC,ETT,DC,DT,TRACE,R3R3,S3S3,
64 . ac,at,bc,bt,bc2,bt2,cc,tt,func,funt,deric,derit,
65 . y,ec2,et2,sigc,sigt,udc,udt,fc,ft,fact,sig0,
66 . fpc,fpt,ccl,ttl,dec,det,hc,ht,kf,kmax,hdc,hdt,
67 . tan_phi,rfac,rfat,r1,r2,r3,s1,s2,s3,t1,t2,t3,rs1,rs2,rs3,
68 . r12,s12,r22,s22,zerostress,flex1,flex2,yfac(6),phi,gxy,ddec,ddet,
69 . dcc1,dtt1,dcc,dtt,kkc,kkt
70 my_real
71 . ec(mvsiz),et(mvsiz),lc(mvsiz),lt(mvsiz),fn(mvsiz),
72 . yc(mvsiz),yt(mvsiz),degmb(mvsiz),
73 . depsxx(mvsiz),depsyy(mvsiz),depsxy(mvsiz),
74 . signxx(mvsiz),signyy(mvsiz),signxy(mvsiz),
75 . tens(3,mvsiz),emincrl(mvsiz),emintrl(mvsiz),eminsrl(mvsiz),
76 . epsnc(mvsiz),epsnt(mvsiz),epsns(mvsiz),emaxsl(mvsiz),emaxsrl(mvsiz),
77 . emaxcl(mvsiz),emaxcrl(mvsiz),emaxtl(mvsiz),emaxtrl(mvsiz),
78 . epsmaxc(mvsiz),epsmaxrlc(mvsiz),epsnrlc(mvsiz),
79 . slmaxrlc(mvsiz),slmaxrlt(mvsiz),slmaxs(mvsiz),slmins(mvsiz),
80 . epsmaxs(mvsiz),epsmaxrls(mvsiz),epsnrls(mvsiz),slmaxrls(mvsiz) ,
81 . slmaxc(mvsiz),epsmaxt(mvsiz),epsmaxrlt(mvsiz),epsnrlt(mvsiz),
82 . slmaxt(mvsiz),phio(mvsiz),slminc(mvsiz),slmint(mvsiz)
83C=======================================================================
84C--- Initialisations
85C----------------------------------------------------------------------
86 niter = 5
87 flagul = iparam(1)
88 nc = iparam(3)
89 nt = iparam(4)
90
91 lc0 = uparam( 1)
92 lt0 = uparam( 2)
93 dc0 = uparam( 3)
94 dt0 = uparam( 4)
95 hc0 = uparam( 5)
96 ht0 = uparam( 6)
97 kc = uparam( 9)
98 kt = uparam(10)
99 kfc = uparam(11)
100 kft = uparam(12)
101 g0 = uparam(13)
102 gt = uparam(14)
103 gb = uparam(15)
104 tan_lock = uparam(16)
105 visce = uparam(17)
106 viscg = uparam(18)
107 kbc = uparam(19)
108 kbt = uparam(20)
109 zerostress = uparam(24)
110 kkc = uparam(40)
111 kkt = uparam(41)
112 flex1 = uparam(26)
113 flex2 = uparam(27)
114 yfac(1) = uparam(28)
115 yfac(2) = uparam(29)
116 yfac(3) = uparam(30)
117c
118 IF (kbc == zero) THEN
119 ccl = ep20
120 ELSE
121 ccl = kc / kbc
122 ENDIF
123 IF (kbt == zero) THEN
124 ttl = ep20
125 ELSE
126 ttl = kt / kbt
127 ENDIF
128C
129 DO i=1,nel
130 degmb(i) = for(i,1)*exx(i)+for(i,2)*eyy(i)+for(i,3)*exy(i)
131 for(i,1)=zero
132 for(i,2)=zero
133 for(i,3)=zero
134 ENDDO
135C----------------------------------------------------------------------
136c strain rate in fiber coord sys
137 DO i=1,nel
138 r1 = dir1(i,1)
139 s1 = dir1(i,2)
140 r2 = dir2(i,1)
141 s2 = dir2(i,2)
142c
143 t1 = exx(i)
144 t2 = eyy(i)
145 t3 = half*exy(i)
146 depsxx(i) = r1*r1*t1 + s1*s1*t2 + two*r1*s1*t3 ! delta_eps_x
147 depsyy(i) = r2*r2*t1 + s2*s2*t2 + two*r2*s2*t3 ! delta_eps_y
148 depsxy(i) = (r1*r2 + s1*s2) / (r1*s2 - r2*s1) ! Tan(alpha_totale)
149C--- integration des deformations
150 etc = uvar(i,4) + depsxx(i)
151 ett = uvar(i,5) + depsyy(i)
152 uvar(i,4) = etc
153 uvar(i,5) = ett
154 ec(i) = exp(etc) - one ! eng strain
155 et(i) = exp(ett) - one
156 ENDDO
157c---
158 IF(flagul /= 1)then!no hysteresis
159c---
160 DO i=1,nel
161 lc(i) = lc0 * (one + ec(i))
162 lt(i) = lt0 * (one + et(i))
163 ENDDO
164C--- Calcul YC, YT
165 ac = kc*dc0
166 ac = ac*ac
167 at = kt*dt0
168 at = at*at
169 DO i=1,nel
170 yc(i) = uvar(i,7)
171 yt(i) = uvar(i,8)
172 fn(i) = zero
173C--- non coupled model
174 DO iter = 1, niter
175 hc = hc0 + yc(i)
176 ht = ht0 + yt(i)
177 dc = sqrt(lc(i)*lc(i) + hc*hc)
178 dt = sqrt(lt(i)*lt(i) + ht*ht)
179 udc= one / dc
180 udt= one / dt
181 hdc= hc * udc
182 hdt= ht * udt
183 cc = dc - dc0
184 tt = dt - dt0
185 IF(ifunc(1) /= 0)THEN
186 fc = yfac(1)*finter(ifunc(1),cc,npf,tf,fpc)
187 fpc = fpc *yfac(1)
188 kc = fpc
189 kfc = flex1*fpc*hc0/dc0
190 fpc = fpc*hdc
191 ELSEIF (cc >= ccl ) THEN
192 fc = half * kc*ccl
193 fpc = zero
194 ELSE
195 fc = (kc - half*kbc * cc) * cc
196 fpc = (kc - kbc*cc) * hdc
197 ENDIF
198 IF(ifunc(2) /= 0)THEN
199 ft = yfac(2)*finter(ifunc(2),tt,npf,tf,fpt)
200 fpt = fpt *yfac(2)
201 kt = fpt
202 kft = flex2*fpt*ht0/dt0
203 fpt = fpt*hdt
204 ELSEIF (tt >= ttl ) THEN
205 ft = half * kt*ttl
206 fpt = zero
207 ELSE
208 ft = (kt - half*kbt * tt) * tt
209 fpt = (kt - kbt*tt) * hdt
210 ENDIF
211 func = kfc*yc(i) + fc*hc/dc
212 funt = kft*yt(i) + ft*ht/dt
213 deric= kfc + fpc*hdc + fc*udc*(one - hdc*hdc)
214 derit= kft + fpt*hdt + ft*udt*(one - hdt*hdt)
215C-------------------------------------
216 yc(i) = yc(i) - func / deric
217 yt(i) = yt(i) - funt / derit
218 fn(i) = zero
219 ENDDO
220C--- Coupled model (traction fibre)
221 IF ((yc(i) + yt(i)) < zero) THEN
222 y = half * (uvar(i,7) - uvar(i,8))
223 DO iter = 1, niter
224 hc = hc0 + y
225 ht = ht0 - y
226 dc = sqrt(lc(i)*lc(i) + hc*hc)
227 dt = sqrt(lt(i)*lt(i) + ht*ht)
228 cc = dc - dc0
229 tt = dt - dt0
230 udc= one / dc
231 udt= one / dt
232 hdc= hc * udc
233 hdt= ht * udt
234 IF(ifunc(1) /= 0)THEN
235 fc =yfac(1)* finter(ifunc(1),cc,npf,tf,fpc)
236 fpc = fpc *yfac(1)
237 kc = fpc
238 kfc = flex1*fpc*hc0/dc0
239 fpc = fpc*hdc
240 ELSEIF (cc >= ccl ) THEN
241 fc = half * kc*ccl
242 fpc = zero
243 ELSE
244 fc = (kc - half*kbc * cc) * cc
245 fpc = (kc - kbc*cc) * hdc
246 ENDIF
247 IF(ifunc(2) /= 0)THEN
248 ft = yfac(2)*finter(ifunc(2),tt,npf,tf,fpt)
249 fpt = fpt *yfac(2)
250 kt = fpt
251 kft = flex2*fpt*ht0/dt0
252 fpt = fpt*hdt
253 ELSEIF (tt >= ttl ) THEN
254 ft = half * kt*ttl
255 fpt = zero
256 ELSE
257 ft = (kt - half*kbt * tt) * tt
258 fpt = (kt - kbt*tt) * hdt
259 ENDIF
260C
261 kf = kfc + kft
262 func = kf*y + fc * hdc - ft * hdt
263 deric = kf + fpc*hdc + fc*udc*(one - hdc*hdc)
264 . + fpt*hdt + ft*udt*(one - hdt*hdt)
265 y = y - func / deric
266C
267 IF (y > 0) THEN
268 y = min(y, ht0)
269 ELSE
270 y = max(y,-hc0)
271 ENDIF
272 ENDDO
273 yc(i) = y
274 yt(i) =-y
275 fn(i) = fc*hc/dc + ft*ht/dt
276 ENDIF
277 hc = hc0 + yc(i)
278 ht = ht0 + yt(i)
279 dc = sqrt(lc(i)*lc(i) + hc*hc)
280 dt = sqrt(lt(i)*lt(i) + ht*ht)
281 ENDDO
282C---
283 DO i=1,nel
284 hc = hc0 + yc(i)
285 ht = ht0 + yt(i)
286 dc = sqrt(lc(i)*lc(i) + hc*hc)
287 dt = sqrt(lt(i)*lt(i) + ht*ht)
288 cc = dc - dc0
289 tt = dt - dt0
290 IF(ifunc(1) /= 0)THEN
291 fc =yfac(1)* finter(ifunc(1),cc,npf,tf,fpc)
292 fpc = fpc *yfac(1)
293 kc = fpc
294 ELSEIF (cc >= ccl ) THEN
295 fc = half * kc*ccl
296 ELSE
297 fc = (kc - half*kbc * cc) * cc
298 ENDIF
299 IF(ifunc(2) /= 0)THEN
300 ft = yfac(2)*finter(ifunc(2),tt,npf,tf,fpt)
301 fpt = fpt *yfac(2)
302 kt = fpt
303 ELSEIF (tt >= ttl ) THEN
304 ft = half * kt*ttl
305 ELSE
306 ft = (kt - half*kbt * tt) * tt
307 ENDIF
308C
309 trace = exp(gstr(i,1) + gstr(i,2))
310 IF (trace /= zero) THEN
311 ec2 = trace / max((ec(i) + one), em6)
312 et2 = trace / max((et(i) + one), em6)
313 rfac = nc / ec2
314 rfat = nt / et2
315 ELSE
316 rfac = zero
317 rfat = zero
318 ENDIF
319C--- contraintes membrane
320 sigc = rfac * fc * lc(i) / dc
321 sigt = rfat * ft * lt(i) / dt
322 signxx(i) = sigc
323 signyy(i) = sigt
324C--- contraintes cisaillement
325 sig0 = uvar(i,10)
326 tan_phi = depsxy(i)
327 phi = atan(tan_phi)*hundred80/pi
328 IF (ifunc(3) /= 0)THEN
329 signxy(i) =yfac(3)* finter(ifunc(3),phi,npf,tf,gxy)
330 gxy = gxy *yfac(3)
331 ELSEIF (tan_phi > tan_lock) THEN
332 signxy(i) = gt*tan_phi + gb - sig0
333 ELSEIF (tan_phi < -tan_lock) THEN
334 signxy(i) = gt*tan_phi - gb - sig0
335 ELSE
336 signxy(i) = g0*tan_phi - sig0
337 ENDIF
338 ENDDO
339C---
340 DO i=1,nel
341 uvar(i,6) = depsxy(i)
342 uvar(i,7) = yc(i)
343 uvar(i,8) = yt(i)
344 ENDDO
345C----------------
346C---------------- FIN LAW 58 sans hysteresis
347C----------------
348 ELSE
349 DO i=1,nel
350 lc(i) = one + ec(i)
351 lt(i) = one + et(i)
352 epsmaxc(i) = uvar(i,17)
353 emincrl(i) = uvar(i,18)
354 epsmaxrlc(i) = uvar(i,19)
355 slmaxc(i) = uvar(i,20)
356 epsmaxt(i) = uvar(i,21)
357 emintrl(i) = uvar(i,22)
358 epsmaxrlt(i) = uvar(i,23)
359 slmaxt(i) = uvar(i,24)
360 slminc(i) = uvar(i,25)
361 slmint(i) = uvar(i,26)
362 slmaxrlc(i) = uvar(i,27)
363 slmaxrlt(i) = uvar(i,28)
364 epsmaxs(i) = uvar(i,30)
365 eminsrl(i) = uvar(i,31)
366 slmaxs(i) = uvar(i,32)
367 slmins(i) = uvar(i,33)
368 epsmaxrls(i) = uvar(i,34)
369 slmaxrls(i) = uvar(i,35)
370 ENDDO
371 DO i=1,nel
372 fn(i) = zero
373c--- uncoupled model (compression of fiber)
374 yc(i) = uvar(i,7)
375 yt(i) = uvar(i,8)
376 hc = hc0 + yc(i) !longueur ressort
377 ht = ht0 + yt(i)
378 dc = sqrt(lc(i)*lc(i) + hc*hc) ! long fibre
379 dt = sqrt(lt(i)*lt(i) + ht*ht)
380 dcc = uvar(i,15) - dc0
381 dtt = uvar(i,16) - dt0
382
383 ddec = dc - uvar(i,15)
384 ddet = dt - uvar(i,16)
385 IF(dcc >=zero ) THEN
386 !----------------------LOAD
387 DO iter = 1, niter
388 hc = hc0 + yc(i) !longueur ressort
389 dc = sqrt(lc(i)*lc(i) + hc*hc) ! long fibre
390 udc = one / dc
391 hdc = hc * udc !sinus(alpha)
392 dcc = dc - dc0
393 dcc = max(uvar(i,18),dcc)
394 dc = dcc + dc0
395 udc = one / dc
396 hdc = hc * udc !sinus(alpha)
397 fc = yfac(1) * finter(ifunc(1),dcc,npf,tf,fpc)
398 fpc = fpc * yfac(1)
399 kc = fpc
400 fpc = fpc * hdc
401 epsmaxc(i) = dcc
402 slmaxc(i) = fc
403 slmaxrlc(i) = slmaxc(i)
404 epsmaxrlc(i)= epsmaxc(i)! UVAR(I,19))
405 func = kfc*yc(i) + fc*hc/dc
406 deric= kfc + fpc*hdc + fc*udc*(one - hdc*hdc)
407 yc(i) = yc(i) - func / deric
408 ENDDO !iter
409 ELSE ! DCC < ZERO domaine de compression
410 DO iter = 1, niter
411 hc = hc0 + yc(i) !longueur ressort
412 dc = sqrt(lc(i)*lc(i) + hc*hc) ! long fibre
413 udc = one / dc
414 hdc = hc * udc !sinus(alpha)
415 dcc = dc - dc0
416 ! compression fibre
417 fc = kkc * dcc
418 kc = kkc
419 fpc = kkc * hdc! derivee % yc
420 func = kfc*yc(i) + fc*hc/dc
421 deric= kfc + fpc*hdc + fc*udc*(one - hdc*hdc)
422 yc(i) = yc(i) - func / deric
423 fn(i) = zero
424 ENDDO !iter
425 epsmaxc(i) = zero
426 emincrl(i) = zero
427 epsmaxrlc(i) = zero
428 slmaxc(i) = zero
429 slminc(i) = zero
430 slmaxrlc(i) = zero
431 ENDIF ! COMPRESSION OR TENSION FOR DIRECTION CHAINE
432c----------------------------------
433c----------------direction trame---
434c----------------------------------
435 IF(dtt >=zero ) THEN
436 !----------------------LOAD
437 DO iter = 1, niter
438 ht = ht0 + yt(i) !longueur ressort
439 dt = sqrt(lt(i)*lt(i) + ht*ht) ! long fibre
440 udt = one / dt
441 hdt = ht * udt !sinus(alpha)
442 dtt = dt - dt0
443 dtt = max(uvar(i,22),dtt)
444 dt = dtt + dt0
445 udt = one / dt
446 hdt = ht * udt !sinus(alpha)
447 ft = yfac(2)*finter(ifunc(2),dtt,npf,tf,fpt)
448 fpt = fpt * yfac(2)
449 kt = fpt
450 fpt = fpt *hdt
451 epsmaxt(i) = dtt
452 slmaxt(i) = ft
453 slmaxrlt(i) = slmaxt(i)
454 epsmaxrlt(i)= epsmaxt(i)!
455 emintrl(i)= zero
456 slmint(i) = zero
457 funt = kft*yt(i) + ft*ht/dt
458 derit= kft + fpt*hdt + ft*udt*(one - hdt*hdt)
459 yt(i) = yt(i) - funt / derit
460 ENDDO !iter
461 ELSE ! DTT < ZERO domaine de compression
462 DO iter = 1, niter
463 ht = ht0 + yt(i) !longueur ressort
464 dt = sqrt(lt(i)*lt(i) + ht*ht) ! long fibre
465 udt = one / dt
466 hdt = ht * udt !sinus(alpha)
467 dtt = dt - dt0
468 ft = kkt * dtt
469 fpt = kkt * hdt! derivee % yT
470 kt = kkt
471 funt = kft*yt(i) + ft*ht/dt
472 derit= kft + fpt*hdt + ft*udt*(one - hdt*hdt)
473 yt(i) = yt(i) - funt / derit
474 ENDDO !iter
475 epsmaxt(i) = zero
476 emintrl(i) = zero
477 epsmaxrlt(i) = zero
478 slmaxt(i) = zero
479 slmint(i) = zero
480 slmaxrlt(i) = zero
481 ENDIF ! COMPRESSION OR TENSION FOR DIRECTION CHAINE
482 fn(i) = zero
483
484C ========================================================
485C--- Coupled model (traction fibre)
486C ========================================================
487 IF ((yc(i) + yt(i)) < zero) THEN !eps
488 hc = hc0 + uvar(i,7)
489 ht = ht0 + uvar(i,8)
490 dc = sqrt(lc(i)*lc(i) + hc*hc)
491 dt = sqrt(lt(i)*lt(i) + ht*ht)
492 ddec = dc - uvar(i,15)
493 ddet = dt - uvar(i,16)
494
495 y = half * (uvar(i,7) - uvar(i,8))
496 hc = hc0 + y
497 ht = ht0 - y
498 dc = sqrt(lc(i)*lc(i) + hc*hc)
499 dt = sqrt(lt(i)*lt(i) + ht*ht)
500 dcc1 = uvar(i,15) - dc0
501 dtt1 = uvar(i,16) - dt0
502 kf = kfc + kft
503 DO iter = 1, niter
504 hc = hc0 + y
505 ht = ht0 - y
506 dc = sqrt(lc(i)*lc(i) + hc*hc)
507 dt = sqrt(lt(i)*lt(i) + ht*ht)
508 dcc = dc - dc0
509 dtt = dt - dt0
510 udc= one / dc
511 udt= one / dt
512 hdc= hc * udc
513 hdt= ht * udt
514 IF(dcc1 >=zero ) THEN
515c----------------------------------LOAD
516 fc = yfac(1) * finter(ifunc(1),dcc,npf,tf,fpc)
517 fpc = fpc * yfac(1)
518 kc = fpc
519 fpc = fpc * hdc
520 epsmaxc(i) = dcc
521 slmaxc(i) = fc
522 slmaxrlc(i) = slmaxc(i)
523 epsmaxrlc(i)= epsmaxc(i)
524 emincrl(i)= zero
525 slminc(i) = zero
526 ELSE ! DOMAINE POSITIF
527 fc = kkc * dcc
528 fpc = kkc * hdc! derivee % yc
529 kc = kkc
530 epsmaxc(i) = zero
531 emincrl(i) = zero
532 epsmaxrlc(i) = zero
533 slmaxc(i) = zero
534 slminc(i) = zero
535 slmaxrlc(i) = zero
536 ENDIF ! (DCC1 >=ZERO )
537c----------------------------------
538c----------------direction trame---
539c----------------------------------
540
541 IF(dtt1 >=zero ) THEN
542 ft = yfac(2)*finter(ifunc(2),dtt,npf,tf,fpt)
543 fpt = fpt * yfac(2)
544 kt = fpt
545 fpt = fpt * hdt
546 epsmaxt(i) = dtt
547 slmaxt(i) = ft
548 slmaxrlt(i) = slmaxt(i)
549 epsmaxrlt(i)= epsmaxt(i)!
550 emintrl(i)= zero
551 slmint(i) = zero
552 ELSE ! DOMAINE POSITIF
553 ! compression fibre
554 ft = kkt * dtt
555 fpt = kkt * hdt! derivee % yT
556 kt = kkt
557 epsmaxt(i) = zero
558 emintrl(i) = zero
559 epsmaxrlt(i) = zero
560 slmaxt(i) = zero
561 slmint(i) = zero
562 slmaxrlt(i) = zero
563 ENDIF
564C
565 func = kf*y + fc * hdc - ft * hdt
566 deric = kf + fpc*hdc + fc*udc*(one - hdc*hdc)
567 . + fpt*hdt + ft*udt*(one - hdt*hdt)
568 y = y - func / deric
569C
570 IF (y > 0) THEN
571 y = min(y, ht0)
572 ELSE
573 y = max(y,-hc0)
574 ENDIF
575 ENDDO !iter
576 yc(i) = y
577 yt(i) =-y
578 fn(i) = fc*hc/dc + ft*ht/dt
579 ENDIF
580 hc = hc0 + yc(i)
581 ht = ht0 + yt(i)
582 dc = sqrt(lc(i)*lc(i) + hc*hc)
583 dt = sqrt(lt(i)*lt(i) + ht*ht)
584C------------------------------------------------------------------
585C SHEAR
586C------------------------------------------------------------------
587 trace = exp(gstr(i,1) + gstr(i,2))
588 IF (trace /= zero) THEN
589 ec2 = trace / max((ec(i) + one), em6)
590 et2 = trace / max((et(i) + one), em6)
591 rfac = nc / ec2
592 rfat = nt / et2
593 ELSE
594 rfac = zero
595 rfat = zero
596 ENDIF
597C--- contraintes membrane
598 sigc = rfac * fc * lc(i) / dc
599 sigt = rfat * ft * lt(i) / dt
600 signxx(i) = sigc
601 signyy(i) = sigt
602C--- contraintes cisaillement
603 sig0 = uvar(i,10)
604 tan_phi = depsxy(i)
605 phi = atan(tan_phi)*hundred80/pi
606 signxy(i) =yfac(3)* finter(ifunc(3),phi,npf,tf,gxy)
607 uvar(i,15) = dc
608 uvar(i,16) = dt
609 ENDDO
610C---
611 DO i=1,nel
612 uvar(i,6) = depsxy(i)
613 uvar(i,7) = yc(i)
614 uvar(i,8) = yt(i)
615c
616 uvar(i,17) = epsmaxc(i)
617 uvar(i,18) = emincrl(i)
618 uvar(i,19) = epsmaxrlc(i)
619 uvar(i,20) = slmaxc(i)
620 uvar(i,21) = epsmaxt(i)
621 uvar(i,22) = emintrl(i)
622 uvar(i,23) = epsmaxrlt(i)
623 uvar(i,24) = slmaxt(i)
624 uvar(i,25) = slminc(i)
625 uvar(i,26) = slmint(i)
626 uvar(i,27) = slmaxrlc(i)
627 uvar(i,28) = slmaxrlt(i)
628 uvar(i,30) = epsmaxs(i)
629 uvar(i,31) = eminsrl(i)
630 uvar(i,32) = slmaxs(i)
631 uvar(i,33) = slmins(i)
632 uvar(i,34) = epsmaxrls(i)
633 uvar(i,35) = slmaxrls(i)
634 ENDDO
635C----------------
636 ENDIF
637C---------------- FIN hysteresis
638c
639 DO i=1,nel
640 tens(1,i) = signxx(i)
641 tens(2,i) = signyy(i)
642 tens(3,i) = signxy(i)
643 ENDDO
644C---------------- RETOUR REPERE COQUE
645 DO i=1,nel
646 r1 = dir1(i,1)
647 s1 = dir1(i,2)
648 r2 = dir2(i,1)
649 s2 = dir2(i,2)
650 rs1= r1*s1
651 rs2= r2*s2
652 r12= r1*r1
653 r22= r2*r2
654 s12= s1*s1
655 s22= s2*s2
656 rs3 = s1*s2-r1*r2
657 r3r3= one+s1*r2+r1*s2
658 r3r3= half*r3r3
659 s3s3= one-s1*r2-r1*s2
660 s3s3= half*s3s3
661 t1 = tens(1,i)
662 t2 = tens(2,i)
663 t3 = tens(3,i)
664 tens(1,i) = r12*t1 + r22*t2 - rs3*t3
665 tens(2,i) = s12*t1 + s22*t2 + rs3*t3
666 tens(3,i) = rs1*t1 + rs2*t2 + (r3r3 - s3s3)*t3
667 ENDDO
668C-----------------------
669C FORCES ET MOMENTS
670C-----------------------
671 DO i=1,nel
672 for(i,1)=for(i,1) + thk(i)*tens(1,i)
673 for(i,2)=for(i,2) + thk(i)*tens(2,i)
674 for(i,3)=for(i,3) + thk(i)*tens(3,i)
675 ENDDO
676C
677 DO i=1,nel
678 degmb(i) = degmb(i)+
679 + for(i,1)*exx(i)+for(i,2)*eyy(i)+for(i,3)*exy(i)
680 eint(i,1) = eint(i,1) + degmb(i)*half*thk(i)*area(i)
681 ENDDO
682C-----------------------------------------------------------
683C REF-STATE ZEROSTRESS OPTION
684C-----------------------------------------------------------
685 IF(zerostress /= zero)THEN
686 DO i=1,nel
687 uvar(i,11) = signxx(i)
688 uvar(i,12) = signyy(i)
689 uvar(i,13) = signxy(i)
690 ENDDO
691c ELSE
692c DO I=1,NEL
693c DEGMB(I) = DEGMB(I)+
694c + FOR(I,1)*EXX(I)+FOR(I,2)*EYY(I)+FOR(I,3)*EXY(I)
695c EINT(I,1) = EINT(I,1) + DEGMB(I)*HALF*THK(I)*AREA(I)
696c ENDDO
697 ENDIF
698C-----------
699 RETURN
700 END
cm58_refsta
subroutine cm58_refsta(nel, nuparam, niparam, uparam, iparam, for, eint, gstr, thk, dir1, dir2, nuvar, uvar, nbfunc, ifunc, npf, tf, area, exx, eyy, exy)
Definition cm58_refsta.F:35
alpha
#define alpha
Definition eval.h:35
area
subroutine area(d1, x, x2, y, y2, eint, stif0)
Definition hm_read_prop32.F:567
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
for
for(i8=*sizetab-1;i8 >=0;i8--)
Definition mumps_common.c:91