OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps65c.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!|| sigeps65c ../engine/source/materials/mat/mat065/sigeps65c.F
25!||--- called by ------------------------------------------------------
26!|| mulawc ../engine/source/materials/mat_share/mulawc.F90
27!||--- calls -----------------------------------------------------
28!|| eqsolv_2 ../engine/source/materials/mat/mat065/sigeps65c.F
29!|| finter ../engine/source/tools/curve/finter.F
30!|| interstar ../engine/source/materials/mat/mat065/sigeps65c.F
31!||====================================================================
32 SUBROUTINE sigeps65c(
33 1 NEL ,NUPARAM,NUVAR ,MFUNC ,KFUNC ,
34 2 NPF ,NPT0 ,IPT ,IFLAG ,
35 2 TF ,TIME ,TIMESTEP,UPARAM ,RHO0 ,
36 3 AREA ,EINT ,THKLY ,ISRATE ,ASRATE ,
37 4 EPSPXX ,EPSPYY ,EPSPXY ,EPSPYZ ,EPSPZX ,
38 5 DEPSXX ,DEPSYY ,DEPSXY ,DEPSYZ ,DEPSZX ,
39 6 EPSXX ,EPSYY ,EPSXY ,EPSYZ ,EPSZX ,
40 7 SIGOXX ,SIGOYY ,SIGOXY ,SIGOYZ ,SIGOZX ,
41 8 SIGNXX ,SIGNYY ,SIGNXY ,SIGNYZ ,SIGNZX ,
42 9 SIGVXX ,SIGVYY ,SIGVXY ,SIGVYZ ,SIGVZX ,
43 A SOUNDSP,VISCMAX,THK ,PLA ,UVAR ,
44 B OFF ,NGL ,IPM ,MAT ,ETSE ,
45 C GS ,YLD ,EPSP ,DPLA_I )
46C-----------------------------------------------
47C I m p l i c i t T y p e s
48C-----------------------------------------------
49#include "implicit_f.inc"
50C---------+---------+---+---+--------------------------------------------
51C VAR | SIZE |TYD| RW| DEFINITION
52C---------+---------+---+---+--------------------------------------------
53C NEL | 1 | I | R | SIZE OF THE ELEMENT GROUP NEL
54C NUPARAM | 1 | I | R | SIZE OF THE USER PARAMETER ARRAY
55C NUVAR | 1 | I | R | NUMBER OF USER ELEMENT VARIABLES
56C---------+---------+---+---+--------------------------------------------
57C NFUNC | 1 | I | R | NUMBER FUNCTION USED FOR THIS USER LAW
58C IFUNC | NFUNC | I | R | FUNCTION INDEX
59C NPF | * | I | R | FUNCTION ARRAY
60C NPT0 | 1 | I | R | NUMBER OF LAYERS OR INTEGRATION POINTS
61C IPT | 1 | I | R | LAYER OR INTEGRATION POINT NUMBER
62C IFLAG | * | I | R | GEOMETRICAL FLAGS
63C TF | * | F | R | FUNCTION ARRAY
64C---------+---------+---+---+--------------------------------------------
65C TIME | 1 | F | R | CURRENT TIME
66C TIMESTEP| 1 | F | R | CURRENT TIME STEP
67C UPARAM | NUPARAM | F | R | USER MATERIAL PARAMETER ARRAY
68C RHO0 | NEL | F | R | INITIAL DENSITY
69C AREA | NEL | F | R | AREA
70C EINT | 2*NEL | F | R | INTERNAL ENERGY(MEMBRANE,BENDING)
71C THKLY | NEL | F | R | LAYER THICKNESS
72C EPSPXX | NEL | F | R | STRAIN RATE XX
73C EPSPYY | NEL | F | R | STRAIN RATE YY
74C ... | | | |
75C DEPSXX | NEL | F | R | STRAIN INCREMENT XX
76C DEPSYY | NEL | F | R | STRAIN INCREMENT YY
77C ... | | | |
78C EPSXX | NEL | F | R | STRAIN XX
79C EPSYY | NEL | F | R | STRAIN YY
80C ... | | | |
81C SIGOXX | NEL | F | R | OLD ELASTO PLASTIC STRESS XX
82C SIGOYY | NEL | F | R | OLD ELASTO PLASTIC STRESS YY
83C ... | | | |
84C---------+---------+---+---+--------------------------------------------
85C SIGNXX | NEL | F | W | NEW ELASTO PLASTIC STRESS XX
86C SIGNYY | NEL | F | W | NEW ELASTO PLASTIC STRESS YY
87C ... | | | |
88C SIGVXX | NEL | F | W | VISCOUS STRESS XX
89C SIGVYY | NEL | F | W | VISCOUS STRESS YY
90C ... | | | |
91C SOUNDSP | NEL | F | W | SOUND SPEED (NEEDED FOR TIME STEP)
92C VISCMAX | NEL | F | W | MAXIMUM DAMPING MODULUS(NEEDED FOR TIME STEP)
93C---------+---------+---+---+--------------------------------------------
94C THK | NEL | F |R/W| THICKNESS
95C PLA | NEL | F |R/W| PLASTIC STRAIN
96C UVAR |NEL*NUVAR| F |R/W| USER ELEMENT VARIABLE ARRAY
97C OFF | NEL | F |R/W| DELETED ELEMENT FLAG (=1. ON, =0. OFF)
98C---------+---------+---+---+--------------------------------------------
99C C o m m o n B l o c k s
100C-----------------------------------------------
101#include "param_c.inc"
102#include "com01_c.inc"
103C-----------------------------------------------
104C I N P U T A r g u m e n t s
105C-----------------------------------------------
106 INTEGER NEL, NUPARAM, NUVAR, NPT0, IPT,ISRATE,
107 . IFLAG(*),NGL(NEL),MAT(NEL),IPM(NPROPMI,*)
108 my_real ,INTENT(IN) :: ASRATE
109 my_real ,DIMENSION(NEL) ,INTENT(INOUT) :: EPSP !< local strain rate used in equations
110 my_real TIME,TIMESTEP,UPARAM(*),
111 . AREA(NEL),RHO0(NEL),EINT(NEL,2),
112 . THKLY(NEL),PLA(NEL),
113 . EPSPXX(NEL),EPSPYY(NEL),
114 . EPSPXY(NEL),EPSPYZ(NEL),EPSPZX(NEL),
115 . DEPSXX(NEL),DEPSYY(NEL),
116 . DEPSXY(NEL),DEPSYZ(NEL),DEPSZX(NEL),
117 . EPSXX(NEL) ,EPSYY(NEL) ,
118 . EPSXY(NEL) ,EPSYZ(NEL) ,EPSZX(NEL) ,
119 . sigoxx(nel),sigoyy(nel),
120 . sigoxy(nel),sigoyz(nel),sigozx(nel),
121 . gs(*)
122C-----------------------------------------------
123C O U T P U T A r g u m e n t s
124C-----------------------------------------------
125 my_real
126 . signxx(nel),signyy(nel),
127 . signxy(nel),signyz(nel),signzx(nel),
128 . sigvxx(nel),sigvyy(nel),
129 . sigvxy(nel),sigvyz(nel),sigvzx(nel),
130 . soundsp(nel),viscmax(nel),etse(nel),
131 . dpla_i(nel)
132C-----------------------------------------------
133C I N P U T O U T P U T A r g u m e n t s
134C-----------------------------------------------
135 my_real uvar(nel,nuvar), off(nel),thk(nel),yld(nel)
136C-----------------------------------------------
137C VARIABLES FOR FUNCTION INTERPOLATION
138C-----------------------------------------------
139 INTEGER NPF(*), MFUNC, KFUNC(MFUNC)
140 my_real FINTER ,TF(*)
141 EXTERNAL FINTER
142C Y = FINTER(IFUNC(J),X,NPF,TF,DYDX)
143C Y : y = f(x)
144C X : x
145C DYDX : f'(x) = dy/dx
146C IFUNC(J): FUNCTION INDEX
147C J : FIRST(J=1), SECOND(J=2) .. FUNCTION USED FOR THIS LAW
148C NPF,TF : FUNCTION PARAMETER
149C-----------------------------------------------
150c L o c a l V a r i a b l e s
151C-----------------------------------------------
152 INTEGER I,I1,I2,J,J1,J2,K,K1,N0,N1,N2,N3,IADBUF,NFUNC,FUNC,FUND,
153 . np1,np2,nrate
154 INTEGER JJ(NEL),IFUNC(NEL,100),
155 . IPOS(NEL),IAD(NEL),ILEN(NEL),
156 . IPOSC1(NEL),IADC1(NEL),ILENC1(NEL),
157 . IPOSC2(NEL),IADC2(NEL),ILENC2(NEL),
158 . IPOSD1(NEL),IADD1(NEL),ILEND1(NEL),
159 . IPOSD2(NEL),IADD2(NEL),ILEND2(NEL)
160 my_real
161 . R,DEZZ,DPLA,A11,A22,RATE1,RATE2,SIGY1,SIGY2,DYDX,DTDS,ALPHA,
162 . S1,S2,SX,T1,T2,TY,X1,X2,Y1,Y2,A,B,C,EPSEFF,SFC,SFD,EPSPL
163 my_real
164 . E(NEL),A1(NEL),A2(NEL),G(NEL),G3(NEL),HC(NEL),HD(NEL),
165 . nu(nel),nnu1(nel),nu3(nel),epsmax(nel),fac(nel),deplzz(nel),
166 . sigo(nel),svm(nel),sigyld(nel),
167 . dydxc1(nel),dydxc2(nel),yc1(nel),yc2(nel),
168 . dydxd1(nel),dydxd2(nel),yd1(nel),yd2(nel),
169 . depse(nel),epse(nel),depss(nel),epss(nel),
170 . dsigxx(nel),dsigyy(nel),dsigxy(nel),
171 . crit(nel),sigfd(nel),sigfc(nel),dsigv(nel),
172 . sds(nel),dsig2(nel),sv2(nel),dsvm(nel),
173 . epsc1(nel),epsc2(nel),epsd1(nel),epsd2(nel),
174 . dydxc(nel),yc(nel),yfac1(nel),yfac2(nel)
175C======================================================================|
176c NUVAR = 5 + NRATE*4
177c UVAR(I,1) = epsilon star*
178c UVAR(I,2) = elastic strain
179c UVAR(I,3) = Von Mises stress
180c UVAR(I,4) = PLA(I)
181c UVAR(I,5) = EPSP(I)
182c UVAR(I,5+J1) = IPOSC1(I)
183c UVAR(I,5+J2) = IPOSC2(I)
184c UVAR(I,5+NRATE+J1) =EPSC1(I)
185c UVAR(I,5+NRATE+J2) =EPSC2(I)
186c UVAR(I,5+NRATE*2+J1)=IPOSD1(I)
187c UVAR(I,5+NRATE*2+J2)=IPOSD2(I)
188c UVAR(I,5+NRATE*3+J1)=EPSD1(I)
189c UVAR(I,5+NRATE*3+J2)=EPSD2(I)
190C======================================================================|
191 DO i=1,nel
192 nfunc = ipm(10,mat(i))
193 DO j=1,nfunc
194 ifunc(i,j)=ipm(10+j,mat(i))
195 ENDDO
196 ENDDO
197 iadbuf = ipm(7,mat(1)) - 1
198 DO i=1,nel
199 e(i) = uparam(iadbuf+2)
200 g(i) = uparam(iadbuf+3)
201 nu(i) = uparam(iadbuf+4)
202 a1(i) = uparam(iadbuf+5)
203 a2(i) = uparam(iadbuf+6)
204 g3(i) = uparam(iadbuf+8)
205 nnu1(i)= nu(i) / (one - nu(i))
206 nu3(i) = one - nnu1(i)
207 epsmax(i)= uparam(iadbuf+11)
208 ENDDO
209 nrate = nint(uparam(iadbuf+1))
210 n0 = 5
211 n1 = nrate + n0
212 n2 = nrate + n1
213 n3 = nrate + n2
214C------------------------------------------
215 IF (time == zero)THEN
216 IF (isigi == 0) THEN
217 DO i=1,nel
218 DO j=1,10
219 uvar(i,j)=zero
220 ENDDO
221 uvar(i,3) = sqrt(sigoxx(i)*sigoxx(i) + sigoyy(i)*sigoyy(i)
222 . -sigoxx(i)*sigoyy(i) + sigoxy(i)*sigoxy(i)*three)
223 ENDDO
224 ENDIF
225c---- find yield values
226 DO i=1,nel
227 DO i1=1,nrate
228 func = ifunc(i,i1)
229 fund = ifunc(i,i1+nrate)
230 np1 = (npf(func+1)-npf(func))/2
231 np2 = (npf(fund+1)-npf(fund))/2
232c---
233 DO j=3,np1
234 j1=2*(j-2)
235 s1=tf(npf(func)+j1)
236 s2=tf(npf(func)+j1+2)
237 t1=tf(npf(func)+j1+1)
238 t2=tf(npf(func)+j1+3)
239 DO k=3,np2
240 k1=2*(k-2)
241 x1=tf(npf(fund)+k1)
242 x2=tf(npf(fund)+k1+2)
243 y1=tf(npf(fund)+k1+1)
244 y2=tf(npf(fund)+k1+3)
245 IF (y2>=t1 .AND. y1<=t2 .AND. x2>=s1 .AND. x1<=s2) THEN
246 dydx = (y2-y1) / (x2-x1)
247 dtds = (t2-t1) / (s2-s1)
248 IF (dydx > dtds) THEN
249 sx = (t1-y1-dtds*s1+dydx*x1) / (dydx-dtds)
250 ty = t1 + dtds*(sx - s1)
251 IF (ty>=y1 .AND. ty<=y2 .AND. sx>=x1 .AND. sx<=x2)THEN
252 iadbuf = ipm(7,mat(i)) + 13
253 uparam(iadbuf+i1+nrate*2) = ty
254 GOTO 150
255 ENDIF
256 ENDIF
257 ENDIF
258 ENDDO
259 ENDDO
260 150 CONTINUE
261c---
262 ENDDO
263 ENDDO
264 ENDIF
265C-------------------
266C STRAIN RATE
267C-------------------
268 DO i=1,nel
269 epspl = half*(abs(epspxx(i)+epspyy(i))
270 . + sqrt((epspxx(i)-epspyy(i))**2 + epspxy(i)**2))
271 epsp(i) = asrate*epspl + (one-asrate)*epsp(i)
272 ENDDO
273C-------------------
274C STRAIN
275C-------------------
276 DO i=1,nel
277 epss(i) = uvar(i,1)
278 epse(i) = uvar(i,2)
279 ENDDO
280C-------------------
281C ELASTIC STRESS
282C-------------------
283 DO i=1,nel
284 a11 = e(i) * a1(i)
285 a22 = e(i) * a2(i)
286C---
287 dsigxx(i) = a11*depsxx(i)+a22*depsyy(i)
288 dsigyy(i) = a22*depsxx(i)+a11*depsyy(i)
289 dsigxy(i) = g(i)*depsxy(i)
290C
291 signxx(i) = sigoxx(i)+dsigxx(i)
292 signyy(i) = sigoyy(i)+dsigyy(i)
293 signxy(i) = sigoxy(i)+dsigxy(i)
294 signyz(i) = sigoyz(i) + gs(i)*depsyz(i)
295 signzx(i) = sigozx(i) + gs(i)*depszx(i)
296C-------------------
297 crit(i) = sigoxx(i)*depsxx(i) + sigoyy(i)*depsyy(i)
298 . + half*(sigoxy(i)*depsxy(i))
299 dsvm(i) = dsigxx(i)*dsigxx(i) + dsigyy(i)*dsigyy(i)
300 . -dsigxx(i)*dsigyy(i) + dsigxy(i)*dsigxy(i)*three
301 sigo(i) = uvar(i,3)
302 sv2(i) = signxx(i)*signxx(i) + signyy(i)*signyy(i)
303 . -signxx(i)*signyy(i) + signxy(i)*signxy(i)*three
304 svm(i) = sqrt(sv2(i))
305 dsig2(i) = dsigxx(i)*dsigxx(i) + dsigyy(i)*dsigyy(i)
306 . -dsigxx(i)*dsigyy(i) + dsigxy(i)*dsigxy(i)*three
307 sds(i) = (sigoxx(i)*dsigxx(i) + sigoyy(i)*dsigyy(i)
308 . +sigoxy(i)*dsigxy(i)*three)*two
309 . -sigoxx(i)*dsigyy(i) - sigoyy(i)*dsigxx(i)
310C-------------------
311C SOUND SPEED, TANGENT MODULUS
312C-------------------
313c SOUNDSP(I) = SQRT(A11/RHO0(I))
314 soundsp(i) = sqrt(e(i)/rho0(i))
315 viscmax(i) = zero
316 etse(i) = one
317 ENDDO
318C-------------------
319C Interpolation
320C-------------------
321 iadbuf = ipm(7,mat(1))-1+14
322 DO i=1,nel
323 jj(i) = 1
324 DO j=2,nrate-1
325 IF (epsp(i) > uparam(iadbuf+j)) jj(i) = j
326 ENDDO
327C
328 i1 = iadbuf + jj(i)
329 j1 = jj(i)
330 rate1 = uparam(i1)
331 yfac1(i) = uparam(i1+nrate)
332 sigy1 = uparam(i1+nrate*2)*yfac1(i)
333C
334 IF (nrate > 1) THEN
335 i2 = i1 + 1
336 j2 = j1 + 1
337c
338 rate2 = uparam(i2)
339 yfac2(i) = uparam(i2+nrate)
340 fac(i) = (epsp(i) - rate1) / (rate2 - rate1)
341 sigy2 = uparam(i2+nrate*2)*yfac2(i)
342 sigyld(i) = sigy1 + fac(i)*(sigy2-sigy1)
343c loading curves
344 iposc1(i) = nint(uvar(i,n0+j1))
345 iposc2(i) = nint(uvar(i,n0+j2))
346 epsc1(i) = uvar(i,n1+j1)
347 epsc2(i) = uvar(i,n1+j2)
348 iadc1(i) = npf(ipm(10+j1,mat(i)))/2 + 1
349 iadc2(i) = npf(ipm(10+j2,mat(i)))/2 + 1
350 ilenc1(i) = npf(ipm(10+j1,mat(i))+1)/2-iadc1(i)-iposc1(i)
351 ilenc2(i) = npf(ipm(10+j2,mat(i))+1)/2-iadc2(i)-iposc2(i)
352c unloading curves
353 iposd1(i) = nint(uvar(i,n2+j1))
354 iposd2(i) = nint(uvar(i,n2+j2))
355 epsd1(i) = uvar(i,n3+j1)
356 epsd2(i) = uvar(i,n3+j2)
357 iadd1(i) = npf(ipm(10+nrate+j1,mat(i)))/2 + 1
358 iadd2(i) = npf(ipm(10+nrate+j2,mat(i)))/2 + 1
359 ilend1(i) = npf(ipm(10+nrate+j1,mat(i))+1)/2
360 . - iadd1(i) - iposd1(i)
361 ilend2(i) = npf(ipm(10+nrate+j2,mat(i))+1)/2
362 . - iadd2(i) - iposd2(i)
363 ELSE
364 sigyld(i) = sigy1
365 iposc1(i) = nint(uvar(i,n0+j1))
366 iposd1(i) = nint(uvar(i,n2+j1))
367 epsc1(i) = uvar(i,n1+j1)
368 epsd1(i) = uvar(i,n3+j1)
369
370 iadc1(i) = npf(ipm(10+j1,mat(i)))/2 + 1
371 iadd1(i) = npf(ipm(10+nrate+j1,mat(i)))/2 + 1
372 ilenc1(i) = npf(ipm(10+j1,mat(i))+1)/2-iadc1(i)-iposc1(i)
373 ilend1(i) = npf(ipm(10+nrate+j1,mat(i))+1)/2
374 . - iadd1(i) - iposd1(i)
375 ENDIF
376 ENDDO
377C-------------------
378 IF (nrate > 1) THEN
379c Sigma load = f(total strain)
380 CALL interstar(tf ,iadc1,iposc1,ilenc1,nel ,
381 . epss ,epsc1,dydxc1,e ,yc1 ,yfac1 )
382 CALL interstar(tf ,iadc2,iposc2,ilenc2,nel ,
383 . epss ,epsc2,dydxc2,e ,yc2 ,yfac2 )
384c Sigma unload = f(elastic strain)
385 CALL interstar(tf ,iadd1,iposd1,ilend1,nel ,
386 . epse ,epsd1,dydxd1,e ,yd1 ,yfac1 )
387 CALL interstar(tf ,iadd2,iposd2,ilend2,nel ,
388 . epse ,epsd2,dydxd2,e ,yd2 ,yfac2 )
389 ELSE
390 CALL interstar(tf ,iadc1,iposc1,ilenc1,nel ,
391 . epss ,epsc1,dydxc1,e ,yc1 ,yfac1 )
392 CALL interstar(tf ,iadd1,iposd1,ilend1,nel ,
393 . epse ,epsd1,dydxd1,e ,yd1 ,yfac1 )
394 ENDIF
395C-------------------
396 IF (nrate > 1) THEN
397 DO i=1,nel
398 j1 = jj(i)
399 j2 = j1 + 1
400 uvar(i,n0+j1) = iposc1(i)
401 uvar(i,n0+j2) = iposc2(i)
402 uvar(i,n1+j1) = epsc1(i)
403 uvar(i,n1+j2) = epsc2(i)
404 uvar(i,n2+j1) = iposd1(i)
405 uvar(i,n2+j2) = iposd2(i)
406 uvar(i,n3+j1) = epsd1(i)
407 uvar(i,n3+j2) = epsd2(i)
408c
409 yc1(i) = yc1(i)
410 yd1(i) = yd1(i)
411 yc2(i) = yc2(i)
412 yd2(i) = yd2(i)
413 dydxc1(i)= dydxc1(i)
414 dydxd1(i)= dydxd1(i)
415 dydxc2(i)= dydxc2(i)
416 dydxd2(i)= dydxd2(i)
417c
418 hc(i) = dydxc1(i) + (dydxc2(i)-dydxc1(i))*fac(i)
419 hd(i) = dydxd1(i) + (dydxd2(i)-dydxd1(i))*fac(i)
420 sigfc(i) = yc1(i) + (yc2(i)-yc1(i))*fac(i)
421 sigfd(i) = yd1(i) + (yd2(i)-yd1(i))*fac(i)
422c
423 ENDDO
424 ELSE
425 DO i=1,nel
426 j1 = jj(i)
427 uvar(i,n0+j1) = iposc1(i)
428 uvar(i,n1+j1) = epsc1(i)
429 uvar(i,n2+j1) = iposd1(i)
430 uvar(i,n3+j1) = epsd1(i)
431c
432 sigfc(i) = yc1(i)
433 sigfd(i) = yd1(i)
434 hc(i) = dydxc1(i)
435 hd(i) = dydxd1(i)
436 ENDDO
437 ENDIF
438C-------------------
439C PROJECTION
440C-------------------
441
442 DO i=1,nel
443 depss(i) = zero
444 depse(i) = zero
445 deplzz(i) = zero
446c----
447 IF (svm(i) > sigfc(i)) THEN
448c load
449 depss(i) = (svm(i) - sigfc(i)) / (e(i) + hc(i))
450 sfc = sigfc(i) + hc(i)*depss(i)
451 sfd = sigfd(i) + hd(i)*depss(i)
452C
453 a = dsig2(i)
454 b = sds(i)
455 c = sigo(i)*sigo(i) - sfc*sfc
456 CALL eqsolv_2(a,b,c,x1,x2)
457 alpha = max(x1,x2)
458 signxx(i) = sigoxx(i) + dsigxx(i)*alpha
459 signyy(i) = sigoyy(i) + dsigyy(i)*alpha
460 signxy(i) = sigoxy(i) + dsigxy(i)*alpha
461c
462 IF (sfd > sfc) THEN
463 depss(i) = (svm(i) - sigfc(i)) / (g3(i) + hc(i))
464 sigfc(i) = sigfc(i) + hc(i)*depss(i)
465 sigfd(i) = sigfd(i) + hd(i)*depss(i)
466c
467 a = dsig2(i)
468 b = sds(i)
469 c = sigo(i)*sigo(i) - sigfc(i)*sigfc(i)
470 CALL eqsolv_2(a,b,c,x1,x2)
471 alpha = max(x1,x2)
472 signxx(i) = sigoxx(i) + dsigxx(i)*alpha
473 signyy(i) = sigoyy(i) + dsigyy(i)*alpha
474 signxy(i) = sigoxy(i) + dsigxy(i)*alpha
475c
476 dpla_i(i) = depss(i)
477 pla(i) = pla(i) + off(i)*dpla_i(i)
478 deplzz(i) = dpla_i(i)*(signxx(i)+signyy(i))*half*nu3(i)
479 . / max(em20,sigfc(i))
480 etse(i)= hc(i)/(hc(i)+e(i))
481 ELSE
482 sigfc(i) = sfc
483 sigfd(i) = sfd
484 dpla_i(i)= zero
485 ENDIF
486c----
487 ELSEIF (svm(i) < sigfd(i)) THEN
488c unload
489 depss(i) = (svm(i) - sigfd(i)) / (e(i) + hd(i))
490 sigfd(i) = sigfd(i) + hd(i)*depss(i)
491 a = dsig2(i)
492 b = sds(i)
493 c = sigo(i)*sigo(i) - sigfd(i)*sigfd(i)
494 CALL eqsolv_2(a,b,c,x1,x2)
495 x1 = max(x1,zero)
496 x2 = max(x2,zero)
497 alpha = min(x1,x2)
498 signxx(i) = sigoxx(i) + dsigxx(i)*alpha
499 signyy(i) = sigoyy(i) + dsigyy(i)*alpha
500 signxy(i) = sigoxy(i) + dsigxy(i)*alpha
501 dpla_i(i) = zero
502c----
503 ELSE
504C cas elastique, pente E => depss = 0
505 dpla_i(i) = zero
506 ENDIF
507 ENDDO
508c----------------------------------------------
509 DO i=1,nel
510 dezz = (depsxx(i)+depsyy(i))*nnu1(i) + deplzz(i)
511 thk(i)= thk(i) - dezz*thkly(i)*off(i)
512c
513 svm(i) = sqrt(signxx(i)*signxx(i) + signyy(i)*signyy(i)
514 . -signxx(i)*signyy(i) + three*signxy(i)*signxy(i))
515 epss(i) = epss(i) + depss(i)
516 epse(i) = epss(i) - pla(i)
517c
518 uvar(i,1) = epss(i)
519 uvar(i,2) = epse(i)
520 uvar(i,3) = svm(i)
521 uvar(i,4) = pla(i)
522 uvar(i,5) = epsp(i)
523 IF (pla(i) > epsmax(i) .and. off(i) == one) off(i)=four_over_5
524c EPSEFF=sqrt((EPSXX(I)*EPSXX(I)+EPSYY(I) *EPSYY(I)+EPSXY(I)
525c . *EPSXY(I))/(un+deux*NU(I)*NU(I)))
526C
527 ENDDO
528C------------------------------------------
529 RETURN
530 END SUBROUTINE sigeps65c
531
532!||====================================================================
533!|| eqsolv_2 ../engine/source/materials/mat/mat065/sigeps65c.F
534!||--- called by ------------------------------------------------------
535!|| sigeps65 ../engine/source/materials/mat/mat065/sigeps65.F
536!|| sigeps65c ../engine/source/materials/mat/mat065/sigeps65c.F
537!||====================================================================
538 SUBROUTINE eqsolv_2(A,B,C,X1,X2)
539C-----------------------------------------------
540C I m p l i c i t T y p e s
541C-----------------------------------------------
542#include "implicit_f.inc"
543C-----------------------------------------------
544C I N P U T A r g u m e n t s
545C-----------------------------------------------
546 my_real
547 . a,b,c,x
548C-----------------------------------------------
549c L o c a l V a r i a b l e s
550C-----------------------------------------------
551 my_real
552 . delta,x1,x2
553C-----------------------------------------------
554 x1 = zero
555 x2 = zero
556 IF (a /= zero) THEN
557 delta = b*b - four*a*c
558 IF (delta > zero) THEN
559 delta = sqrt(delta)
560 x2 = (-b - delta)*half/a
561 x1 = -b/a - x2
562 ELSEIF (delta == zero) THEN
563 x1 = -b*half/a
564 x2 = x1
565 ENDIF
566 ENDIF
567c-----------
568 RETURN
569 END SUBROUTINE eqsolv_2
570
571!||====================================================================
572!|| interstar ../engine/source/materials/mat/mat065/sigeps65c.F
573!||--- called by ------------------------------------------------------
574!|| sigeps65 ../engine/source/materials/mat/mat065/sigeps65.F
575!|| sigeps65c ../engine/source/materials/mat/mat065/sigeps65c.F
576!||====================================================================
577 SUBROUTINE interstar(TF ,IAD ,IPOS ,ILEN ,NEL ,
578 . X ,EPSS_i,DYDX ,E ,Y ,YFAC )
579C-----------------------------------------------
580c I m p l i c i t T y p e s
581C-----------------------------------------------
582#include "implicit_f.inc"
583 INTEGER ILEN(*),IPOS(*),IAD(*),NEL
584 INTEGER I,J,J1,J2,ICONT
585 my_real
586 . X(*),DYDX(*),Y(*),TF(2,*),EPSS_i(*),E(*)
587 my_real
588 . DX,DY,DERI,X1,X2,Y1,Y2,EPSS_ii(NEL),YFAC(NEL)
589C=======================================================================
590 J = 0
591 100 CONTINUE
592 j = j+1
593 icont = 0
594c
595 DO i=1,nel
596 j1 = ipos(i)+iad(i)
597 j2 = j1 + 1
598 x1 = tf(1,j1)
599 x2 = tf(1,j2)
600 y2 = tf(2,j2)
601 y1 = tf(2,j1)
602 dx = x2 - x1
603 dy = y2 - y1
604 deri = dy*yfac(i) / max(dx, em20)
605 dx = dx*(e(i) - deri)/e(i)
606 epss_ii(i) = epss_i(i) + dx
607 IF (j < ilen(i) .and. x(i) > epss_ii(i)) THEN
608 icont = 1
609 ipos(i) =ipos(i) + 1
610 epss_i(i) = epss_ii(i)
611 ELSEIF (ipos(i) > 0 .and. x(i) < epss_i(i)) THEN
612 icont = 1
613 epss_ii(i) = epss_i(i)
614 x2 = x1
615 y2 = y1
616 ipos(i) = ipos(i) - 1
617 j1 = ipos(i)+iad(i)
618 x1 = tf(1,j1)
619 y1 = tf(2,j1)
620 dx = x2 - x1
621 dy = y2 - y1
622 deri = dy*yfac(i) / max(dx, em20)
623 dx = dx*(e(i) - deri)/e(i)
624 epss_i(i) = epss_ii(i) - dx
625 ENDIF
626 ENDDO
627C
628 IF (icont == 1) GOTO 100
629c-----------
630 DO i=1,nel
631 j1 =ipos(i)+iad(i)
632 j2 = j1+1
633 dydx(i)=(tf(2,j2)-tf(2,j1))*yfac(i)/(epss_ii(i)-epss_i(i))
634 y(i) = tf(2,j1)*yfac(i) + dydx(i)*(x(i)-epss_i(i))
635 ENDDO
636c-----------
637 RETURN
638 END SUBROUTINE interstar
my_real
#define my_real
Definition cppsort.cpp:32
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
eqsolv_2
subroutine eqsolv_2(a, b, c, x1, x2)
Definition sigeps65c.F:539
sigeps65c
subroutine sigeps65c(nel, nuparam, nuvar, mfunc, kfunc, npf, npt0, ipt, iflag, tf, time, timestep, uparam, rho0, area, eint, thkly, israte, asrate, epspxx, epspyy, epspxy, epspyz, epspzx, 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, ipm, mat, etse, gs, yld, epsp, dpla_i)
Definition sigeps65c.F:46
interstar
subroutine interstar(tf, iad, ipos, ilen, nel, x, epss_i, dydx, e, y, yfac)
Definition sigeps65c.F:579