OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
m14law.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!|| m14law ../engine/source/materials/mat/mat014/m14law.F
25!||--- called by ------------------------------------------------------
26!|| mmain ../engine/source/materials/mat_share/mmain.F90
27!||--- calls -----------------------------------------------------
28!|| m14ama ../engine/source/materials/mat/mat014/m14ama.F
29!|| m14ftg ../engine/source/materials/mat/mat014/m14ftg.F
30!|| m14gtf ../engine/source/materials/mat/mat014/m14gtf.F
31!||====================================================================
32 SUBROUTINE m14law(
33 1 PM, OFF, SIG, EINT,
34 2 PLA, SIGF, EPSF, DAM,
35 3 EPE, EPC, A, VOL,
36 4 RX, RY, RZ, SX,
37 5 SY, SZ, MAT, VNEW,
38 6 DVOL, SSP, D1, D2,
39 7 D3, D4, D5, D6,
40 8 SOLD1, SOLD2, SOLD3, SOLD4,
41 9 SOLD5, SOLD6, SIGY, DEFP,
42 A NGL, SEQ_OUTPUT,NEL, TSAIWU,
43 B JCVT, JSPH)
44C-----------------------------------------------
45C I m p l i c i t T y p e s
46C-----------------------------------------------
47#include "implicit_f.inc"
48#include "comlock.inc"
49C-----------------------------------------------
50C G l o b a l P a r a m e t e r s
51C-----------------------------------------------
52#include "mvsiz_p.inc"
53C-----------------------------------------------
54C C o m m o n B l o c k s
55C-----------------------------------------------
56#include "com01_c.inc"
57#include "param_c.inc"
58#include "com08_c.inc"
59#include "scr17_c.inc"
60#include "units_c.inc"
61C-----------------------------------------------
62C D u m m y A r g u m e n t s
63C-----------------------------------------------
64 INTEGER, INTENT(IN) :: JCVT
65 INTEGER, INTENT(IN) :: JSPH
66 INTEGER MAT(MVSIZ),NGL(MVSIZ),NEL
67C REAL
68 my_real
69 . PM(NPROPM,*), OFF(*), SIG(NEL,6), EINT(*), PLA(*), SIGF(*),
70 . RX(*), RY(*), RZ(*), SX(*), SY(*), SZ(*),
71 . EPSF(*), DAM(NEL,5), EPE(NEL,3), EPC(NEL,3), A(MVSIZ,6), VOL(*),
72 . VNEW(MVSIZ), DVOL(MVSIZ), SSP(MVSIZ),
73 . D1(MVSIZ), D2(MVSIZ), D3(MVSIZ),D4(MVSIZ),D5(MVSIZ),D6(MVSIZ),
74 . SOLD1(MVSIZ), SOLD2(MVSIZ), SOLD3(MVSIZ), SOLD4(MVSIZ),
75 . sold5(mvsiz), sold6(mvsiz),sigy(*),defp(*),seq_output(*),
76 . tsaiwu(mvsiz)
77C-----------------------------------------------
78C L o c a l V a r i a b l e s
79C-----------------------------------------------
80 INTEGER KD1(MVSIZ), KD2(MVSIZ), KD3(MVSIZ),
81 . KD4(MVSIZ),ICC(MVSIZ),ICC_1,
82 . I, IDAM, KDX, MX, NINDX, INDEX(MVSIZ), J
83C REAL
84 my_real
85 . AX(MVSIZ), AY(MVSIZ), AZ(MVSIZ), BX(MVSIZ), BY(MVSIZ),
86 . BZ(MVSIZ), CX(MVSIZ), CY(MVSIZ),CZ(MVSIZ),
87 . WVEC(MVSIZ),
88 . T1(MVSIZ), T2(MVSIZ), T3(MVSIZ),T4(MVSIZ),T5(MVSIZ),T6(MVSIZ),
89 . E1(MVSIZ), E2(MVSIZ), E3(MVSIZ),E4(MVSIZ),E5(MVSIZ),E6(MVSIZ),
90 . SIGMY(MVSIZ), ALPHA(MVSIZ), EFIB(MVSIZ),
91 . epsft(mvsiz), epsfc(mvsiz),
92 . d11(mvsiz), d12(mvsiz), d13(mvsiz), d22(mvsiz),
93 . d23(mvsiz), d33(mvsiz), g12(mvsiz), g23(mvsiz), g31(mvsiz),
94 . c11(mvsiz), c22(mvsiz), c33(mvsiz), c12(mvsiz), c23(mvsiz),
95 . c13(mvsiz), f11(mvsiz), f22(mvsiz), f44(mvsiz), f55(mvsiz),
96 . f12(mvsiz), f23(mvsiz), f1(mvsiz), f2(mvsiz), f4(mvsiz),
97 . f5(mvsiz), delta(mvsiz), so1(mvsiz),
98 . so2(mvsiz), so3(mvsiz), so4(mvsiz), so5(mvsiz), so6(mvsiz),
99 . ds1(mvsiz), ds2(mvsiz), ds3(mvsiz), ds4(mvsiz), ds5(mvsiz),
100 . ds6(mvsiz), dp1(mvsiz), dp2(mvsiz), dp3(mvsiz), dp4(mvsiz),
101 . dp5(mvsiz), dp6(mvsiz), lamda(mvsiz), coef(mvsiz), plas(mvsiz),
102 . cn(mvsiz), cb(mvsiz), cnn(mvsiz), sigt1(mvsiz), sigt2(mvsiz),
103 . sigt3(mvsiz),cc(mvsiz),epdr(mvsiz)
104 my_real
105 . fib, ca, sigmx, epsp, dt5 ,sigym,
106 . cb_1,cn_1,cc_1,epdr_1,
107 . f1_1,f2_1,f4_1,f5_1,f11_1,
108 . f22_1,f44_1,f55_1,f12_1,f23_1,
109 . c11_1,c22_1,c33_1,c12_1,c23_1,
110 . c13_1,d11_1,d12_1,d13_1,d22_1,
111 . d23_1,d33_1,g12_1,g23_1,g31_1,
112 . alpha_1,wplaref,dwpla
113C=======================================================================
114 IF (ncycle==0 .AND. irun==1) THEN
115 DO i=1,nel
116 kd1(i) = 0
117 kd2(i) = 0
118 kd3(i) = 0
119 kd4(i) = 0
120 ENDDO
121 ELSE
122 DO i=1,nel
123 idam = int(dam(i,5))-10000
124 kd1(i) = idam/1000
125 kdx = idam - kd1(i)*1000
126 kd2(i) = kdx/100
127 kdx = kdx - kd2(i)*100
128 kd3(i) = kdx/10
129 kd4(i) = kdx - kd3(i)*10
130 ENDDO
131 ENDIF
132C------------------
133C FIBER CONTENT
134C------------------
135 fib = zero
136 mx = mat(1)
137 alpha_1 = pm(39,mx)
138!
139 DO i=1,nel
140 alpha(i)=alpha_1
141 fib=fib+alpha(i)
142 ENDDO
143C
144C--------------------------------------------
145C STRESS TRANSFORMATION (GLOBAL -> FIBER)
146C--------------------------------------------
147 CALL m14ama(
148 1 pm, a, rx, ry,
149 2 rz, sx, sy, sz,
150 3 ax, ay, az, bx,
151 4 by, bz, cx, cy,
152 5 cz, nel, jcvt, jsph)
153 CALL m14gtf(sig,ax ,ay ,az ,bx ,by,
154 2 bz ,cx ,cy ,cz ,d1 ,d2,
155 3 d3 ,d4 ,d5 ,d6 ,t1 ,t2,
156 4 t3 ,t4 ,t5 ,t6 ,e1 ,e2,
157 5 e3 ,e4 ,e5 ,e6 ,nel)
158C---------------------------------
159C MATERIAL PROPERTIES
160C---------------------------------
161 nindx=0
162 mx =mat(1)
163 cb_1 =pm(26,mx)
164 cn_1 =pm(27,mx)
165 cc_1 =pm(50,mx)
166 epdr_1 =pm(51,mx)
167 icc_1 =nint(pm(52,mx))
168!
169 DO i=1,nel
170 ca =pm(25,mx)
171 cb(i) =cb_1
172 cn(i) =cn_1
173 cc(i) =cc_1
174 epdr(i) =epdr_1
175 icc(i) =icc_1
176 epsp = max(abs(e1(i)),abs(e2(i)),abs(e3(i)),abs(e4(i)),abs(e5(i)))
177 IF (epsp>epdr(i) .AND. cc(i)/=zero) THEN
178 epsp=one + cc(i) * log(epdr(i)/epsp)
179 ELSE
180 epsp=one
181 ENDIF
182 IF (icc(i)==1) THEN
183 sigmx =pm(28,mx)*epsp
184 ELSEIF (icc(i)==2) THEN
185 sigmx =pm(28,mx)
186 ELSEIF (icc(i)==3) THEN
187 sigmx =pm(28,mx)*epsp
188C WPLAMX(I) = WPLAMX(I)*EPSP
189 ELSEIF (icc(i)==4) THEN
190 sigmx =pm(28,mx)
191C WPLAMX(I) = WPLAMX(I)*EPSP
192 ENDIF
193 cb(i) =cb(i)*epsp
194 ca =ca*epsp
195 sigmy(i)= min(sigmx,ca+cb(i)*pla(i)**cn(i))
196 IF (sigmy(i)==sigmx .AND. off(i)==one) THEN
197 off(i)=zep99
198 kd4(i)=2
199 nindx=nindx+1
200 index(nindx)=i
201 ENDIF
202 sigt1(i)=pm(33,mx)
203 sigt2(i)=pm(34,mx)
204 sigt3(i)=pm(35,mx)
205 ssp(i) =pm(49,mx)
206 delta(i)=pm(78,mx)
207 ENDDO
208C
209 IF (nindx/=0) THEN
210 DO j=1,nindx
211 i=index(j)
212#include "lockon.inc"
213 WRITE(iout,1000) ngl(i)
214#include "lockoff.inc"
215 END DO
216 END IF
217C
218 DO i=1,nel
219 e1(i)=e1(i)*dt1
220 e2(i)=e2(i)*dt1
221 e3(i)=e3(i)*dt1
222 e4(i)=e4(i)*dt1
223 e5(i)=e5(i)*dt1
224 e6(i)=e6(i)*dt1
225 ENDDO
226C
227 DO i=1,nel
228 epe(i,1)=epe(i,1)+e1(i)
229 epe(i,2)=epe(i,2)+e2(i)
230 epe(i,3)=epe(i,3)+e3(i)
231 ENDDO
232C
233 mx =mat(1)
234 f1_1 =pm(59,mx)
235 f2_1 =pm(60,mx)
236 f4_1 =pm(61,mx)
237 f5_1 =pm(62,mx)
238 f11_1 =pm(63,mx)
239 f22_1 =pm(64,mx)
240 f44_1 =pm(65,mx)
241 f55_1 =pm(66,mx)
242 f12_1 =pm(67,mx)
243 f23_1 =pm(68,mx)
244C
245 c11_1 =pm(69,mx)
246 c22_1 =pm(73,mx)
247 c33_1 =pm(74,mx)
248 c12_1 =pm(75,mx)
249 c23_1 =pm(76,mx)
250 c13_1 =pm(77,mx)
251C
252 d11_1 =pm(40,mx)
253 d12_1 =pm(41,mx)
254 d13_1 =pm(42,mx)
255 d22_1 =pm(43,mx)
256 d23_1 =pm(44,mx)
257 d33_1 =pm(45,mx)
258 g12_1 =pm(46,mx)
259 g23_1 =pm(47,mx)
260 g31_1 =pm(48,mx)
261C
262 DO i=1,nel
263 f1(i) = f1_1
264 f2(i) = f2_1
265 f4(i) = f4_1
266 f5(i) = f5_1
267 f11(i) = f11_1
268 f22(i) = f22_1
269 f44(i) = f44_1
270 f55(i) = f55_1
271 f12(i) = f12_1
272 f23(i) = f23_1
273C
274 c11(i) = c11_1
275 c22(i) = c22_1
276 c33(i) = c33_1
277 c12(i) = c12_1
278 c23(i) = c23_1
279 c13(i) = c13_1
280C
281 d11(i) = d11_1
282 d12(i) = d12_1
283 d13(i) = d13_1
284 d22(i) = d22_1
285 d23(i) = d23_1
286 d33(i) = d33_1
287 g12(i) = g12_1
288 g23(i) = g23_1
289 g31(i) = g31_1
290 ENDDO
291C-------------------------
292C NEW ELASTIC STRESSES
293C-------------------------
294 DO i=1,nel
295 so1(i)=t1(i)
296 so2(i)=t2(i)
297 so3(i)=t3(i)
298 so4(i)=t4(i)
299 so5(i)=t5(i)
300 so6(i)=t6(i)
301 ENDDO
302C
303 DO i=1,nel
304 t1(i)=t1(i)+d11(i)*e1(i)+d12(i)*e2(i)+d13(i)*e3(i)
305 t2(i)=t2(i)+d12(i)*e1(i)+d22(i)*e2(i)+d23(i)*e3(i)
306 t3(i)=t3(i)+d13(i)*e1(i)+d23(i)*e2(i)+d33(i)*e3(i)
307 t4(i)=t4(i)+g12(i)*e4(i)
308 t5(i)=t5(i)+g23(i)*e5(i)
309 t6(i)=t6(i)+g31(i)*e6(i)
310 ENDDO
311C------------------
312C FIBERS STRESS
313C------------------
314 IF (fib>0) THEN
315 mx = mat(1)
316 DO i=1,nel
317 efib(i) = pm(36,mx)
318 epsft(i)= pm(37,mx)
319 epsfc(i)= pm(38,mx)
320 epsf(i) = epsf(i)+ e1(i)
321 sigf(i) = efib(i)*epsf(i)
322 ENDDO
323C
324C DEACTIVATED (DAM4 AND DAM5 TEMPORARILY USED FOR GENERAL DAMAGE INFO)
325C CALL M14FIB(PM,SIGF,EPSF,DAM)
326 ENDIF
327C-------------------------
328C GENERAL FAILURE TEST
329C-------------------------
330 DO i=1,nel
331 IF (off(i)<em01) off(i)=zero
332 IF (off(i)<one) off(i)=off(i)*four_over_5
333 ENDDO
334C-------------------
335C TENSION DAMAGE
336C-------------------
337 DO i=1,nel
338C------------------------
339C DIRECTION 1 (FIBER)
340C------------------------
341 wvec(i)=(one-dam(i,1))*sigt1(i)
342 IF (t1(i)>wvec(i)) THEN
343C
344 IF (epc(i,1)==zero) THEN
345 epc(i,1) = max(epe(i,1),zero)
346 ELSE
347 epc(i,1) = max(epc(i,1)+e1(i),zero)
348 ENDIF
349C
350 t1(i)=wvec(i)
351 t2(i)=t2(i)-d12(i)*e1(i)*dam(i,1)
352 t3(i)=t3(i)-d13(i)*e1(i)*dam(i,1)
353 IF (kd1(i)==0) THEN
354 kd1(i)=1
355C* WRITE(6,6001) NCYCLE,I+NFT,'1/1',
356C* 1 T1(I),T2(I),T3(I),T4(I),T5(I),T6(I)
357 ENDIF
358 dam(i,1)= min((dam(i,1)+delta(i)),one)
359 IF (dam(i,1)>=1. .AND. kd1(i)/=2) THEN
360 kd1(i)=2
361C* WRITE(6,6001) NCYCLE,I+NFT,'1/2',
362C* 1 T1(I),T2(I),T3(I),T4(I),T5(I),T6(I)
363 ENDIF
364 ENDIF
365C
366 IF( e1(i)<zero .AND. dam(i,1)>zero)
367 . epc(i,1)= max(epc(i,1)+e1(i),zero)
368C----------------
369C DIRECTION 2
370C----------------
371 wvec(i)=(one-dam(i,2))*sigt2(i)
372 IF (t2(i)>wvec(i)) THEN
373C
374 IF (epc(i,2)==zero) THEN
375 epc(i,2)= max(epe(i,2),zero)
376 ELSE
377 epc(i,2)= max(epc(i,2)+e2(i),zero)
378 ENDIF
379C
380 t1(i)=t1(i)-d12(i)*e2(i)*dam(i,2)
381 t2(i)=wvec(i)
382 t3(i)=t3(i)-d23(i)*e2(i)*dam(i,2)
383 IF (kd2(i)==0) THEN
384 kd2(i)=1
385C* WRITE(6,6001) NCYCLE,I+NFT,'2/1',
386C* 1 T1(I),T2(I),T3(I),T4(I),T5(I),T6(I)
387 ENDIF
388 dam(i,2)= min((dam(i,2)+delta(i)),one)
389 IF (dam(i,2)>=one.AND. kd2(i)/=2) THEN
390 kd2(i)=2
391C* WRITE(6,6001) NCYCLE,I+NFT,'2/2',
392C* 1 T1(I),T2(I),T3(I),T4(I),T5(I),T6(I)
393 ENDIF
394 ENDIF
395C
396 IF (e2(i)<zero .AND. dam(i,2)>zero)
397 . epc(i,2)= max(epc(i,2)+e2(i),zero)
398C----------------
399C DIRECTION 3
400C----------------
401 wvec(i)=(one-dam(i,3))*sigt3(i)
402 IF (t3(i)>wvec(i)) THEN
403C
404 IF (epc(i,3)==zero) THEN
405 epc(i,3)= max(epe(i,3),zero)
406 ELSE
407 epc(i,3)= max(epc(i,3)+e3(i),zero)
408 ENDIF
409C
410 t1(i)=t1(i)-d13(i)*e3(i)*dam(i,3)
411 t2(i)=t2(i)-d23(i)*e3(i)*dam(i,3)
412 t3(i)=wvec(i)
413 IF (kd3(i)==0) THEN
414 kd3(i)=1
415C* WRITE(6,6001) NCYCLE,I+NFT,'3/1',
416C* 1 T1(I),T2(I),T3(I),T4(I),T5(I),T6(I)
417 ENDIF
418 dam(i,3)= min((dam(i,3)+delta(i)),one)
419 IF (dam(i,3)>=1. .AND. kd3(i)/=2) THEN
420 kd3(i)=2
421C* WRITE(6,6001) NCYCLE,I+NFT,'3/2',
422C* 1 T1(I),T2(I),T3(I),T4(I),T5(I),T6(I)
423 ENDIF
424 ENDIF
425C
426 IF (e3(i)<zero .AND. dam(i,3)>zero)
427 . epc(i,3)= max(epc(i,3)+e3(i),zero)
428C
429 ENDDO
430C-----------------------------------
431C CRACK OPEN --> NO COMPRESSION
432C-----------------------------------
433 DO i=1,nel
434 IF (t1(i)<zero.AND.epc(i,1)>zero) THEN
435 t1(i)=zero
436 t2(i)=t2(i)-d12(i)*e1(i)*dam(i,1)
437 t3(i)=t3(i)-d13(i)*e1(i)*dam(i,1)
438 ENDIF
439 IF (t2(i)<zero.AND.epc(i,2)>zero) THEN
440 t1(i)=t1(i)-d12(i)*e2(i)*dam(i,2)
441 t2(i)=zero
442 t3(i)=t3(i)-d23(i)*e2(i)*dam(i,2)
443 ENDIF
444 IF (t3(i)<zero.AND.epc(i,3)>zero) THEN
445 t3(i)=zero
446 t1(i)=t1(i)-d13(i)*e3(i)*dam(i,3)
447 t2(i)=t2(i)-d23(i)*e3(i)*dam(i,3)
448 ENDIF
449 ENDDO
450C---------------------
451C PLASTICITY START
452C---------------------
453 DO i=1,nel
454 wvec(i)=f1(i)*t1(i)+f2(i)*(t2(i)+t3(i))+
455 . f11(i)*t1(i)*t1(i)+f22(i)*(t2(i)*t2(i)+t3(i)*t3(i))+
456 . f55(i)*t5(i)*t5(i)+f44(i)*(t4(i)*t4(i)+t6(i)*t6(i))+
457 . two*f12(i)*(t1(i)*t2(i)+t1(i)*t3(i))+2*f23(i)*t2(i)*t3(i)
458 dam(i,4)=wvec(i)
459 tsaiwu(i)=max(min(wvec(i)/sigmy(i),one),tsaiwu(i))
460!! SEQ_OUTPUT(I) = WVEC(I)
461 ENDDO
462C
463 DO i=1,nel
464 coef(i)=zero
465 IF (wvec(i)>sigmy(i).AND.off(i)==one) THEN
466 coef(i)=one
467 IF (kd4(i)==0) THEN
468 kd4(i)=1
469C* WRITE(6,6001) NCYCLE,I+NFT,'4/1',
470C* 1 T1(I),T2(I),T3(I),T4(I),T5(I),T6(I)
471 ENDIF
472 ENDIF
473 ENDDO
474C
475 DO i=1,nel
476 dp1(i)=f1(i)+two*f11(i)*so1(i)+two*f12(i)*(so2(i)+so3(i))
477 dp2(i)=f2(i)+two*f22(i)*so2(i)+two*f12(i)*so1(i)
478 . + so3(i)*f23(i)*two
479 dp3(i)=f2(i)+two*f22(i)*so3(i)+two*f12(i)*so1(i)
480 . + so2(i)*f23(i)*two
481 dp4(i)=two*f44(i)*so4(i)
482 dp5(i)=two*f55(i)*so5(i)
483 dp6(i)=two*f44(i)*so6(i)
484 ENDDO
485C
486 DO i=1,nel
487 ds1(i)=t1(i)-so1(i)
488 ds2(i)=t2(i)-so2(i)
489 ds3(i)=t3(i)-so3(i)
490 ds4(i)=t4(i)-so4(i)
491 ds5(i)=t5(i)-so5(i)
492 ds6(i)=t6(i)-so6(i)
493 ENDDO
494C
495 DO i=1,nel
496 lamda(i)=(dp1(i)*ds1(i)+dp2(i)*ds2(i)+dp3(i)*ds3(i)
497 . +dp4(i)*ds4(i)+dp5(i)*ds5(i)+dp6(i)*ds6(i))*coef(i)
498 ENDDO
499C
500 DO i=1,nel
501 cnn(i)=cn(i)-one
502 ENDDO
503C
504 DO i=1,nel
505 plas(i)=one
506 IF (pla(i)>zero) plas(i)=pla(i)**cnn(i)
507 ENDDO
508C
509 DO 208 i=1,nel
510 IF (lamda(i)==zero) GO TO 208
511 lamda(i)=lamda(i)*coef(i)/
512 . (dp1(i)*(d11(i)*dp1(i)+d12(i)*dp2(i)+d13(i)*dp3(i))+
513 . dp2(i)*(d12(i)*dp1(i)+d22(i)*dp2(i)+d23(i)*dp3(i))+
514 . dp3(i)*(d13(i)*dp1(i)+d23(i)*dp2(i)+d33(i)*dp3(i))+
515 . two*dp4(i)*g12(i)*dp4(i)+
516 . two*dp5(i)*g23(i)*dp5(i)+
517 . two*dp6(i)*g31(i)*dp6(i)+
518 . (so1(i)*dp1(i)+so2(i)*dp2(i)+so3(i)*dp3(i)+
519 . two*so4(i)*dp4(i)+2.*so5(i)*dp5(i)+2.*so6(i)*dp6(i))
520 . *cn(i)*cb(i)*plas(i))
521 208 CONTINUE
522C WRITE(31,3104) LAMDA(1)
523 3104 FORMAT(' 208 LAMDA ',e11.4)
524C
525 DO i=1,nel
526 dp1(i)=lamda(i)*dp1(i)
527 dp2(i)=lamda(i)*dp2(i)
528 dp3(i)=lamda(i)*dp3(i)
529 dp4(i)=lamda(i)*dp4(i)
530 dp5(i)=lamda(i)*dp5(i)
531 dp6(i)=lamda(i)*dp6(i)
532 ENDDO
533C
534 DO i=1,nel
535 epe(i,1)=epe(i,1)-dp1(i)
536 epe(i,2)=epe(i,2)-dp2(i)
537 epe(i,3)=epe(i,3)-dp3(i)
538 ENDDO
539C
540 DO i=1,nel
541 t1(i)=t1(i)-d11(i)*dp1(i)-d12(i)*dp2(i)-d13(i)*dp3(i)
542 t2(i)=t2(i)-d12(i)*dp1(i)-d22(i)*dp2(i)-d23(i)*dp3(i)
543 t3(i)=t3(i)-d13(i)*dp1(i)-d23(i)*dp2(i)-d33(i)*dp3(i)
544 t4(i)=t4(i)-g12(i)*dp4(i)*two
545 t5(i)=t5(i)-g23(i)*dp5(i)*two
546 t6(i)=t6(i)-g31(i)*dp6(i)*two
547 ENDDO
548C
549 mx = mat(1)
550 wplaref = pm(98 ,mx)
551 DO i=1,nel
552 dwpla = half*
553 . (dp1(i)*(t1(i)+so1(i))+
554 . dp2(i)*(t2(i)+so2(i))+
555 . dp3(i)*(t3(i)+so3(i))+
556 . two*dp4(i)*(t4(i)+so4(i))+
557 . two*dp5(i)*(t5(i)+so5(i))+
558 . two*dp6(i)*(t6(i)+so6(i)))
559 dwpla = max(dwpla ,zero) / wplaref
560 pla(i) = pla(i) + dwpla
561 pla(i)= max(pla(i),zero)
562 ENDDO
563C-------------------
564C PLASTICITY END
565C-------------------
566C--------------------------------------------
567C STRESS TRANSFORMATION (FIBER -> GLOBAL)
568C--------------------------------------------
569 CALL m14ftg(sig,ax ,ay ,az ,bx ,by ,
570 2 bz ,cx ,cy ,cz ,t1 ,t2 ,
571 3 t3 ,t4 ,t5 ,t6 ,nel)
572C
573 DO i=1,nel
574 sig(i,1)=sig(i,1)*off(i)
575 sig(i,2)=sig(i,2)*off(i)
576 sig(i,3)=sig(i,3)*off(i)
577 sig(i,4)=sig(i,4)*off(i)
578 sig(i,5)=sig(i,5)*off(i)
579 sig(i,6)=sig(i,6)*off(i)
580 ENDDO
581C--------------------
582C INTERNAL ENERGY
583C--------------------
584 dt5=half*dt1
585 DO i=1,nel
586 eint(i)=eint(i)+dt5*vnew(i)*
587 . ( d1(i)*(sold1(i)+sig(i,1))
588 . + d2(i)*(sold2(i)+sig(i,2))
589 . + d3(i)*(sold3(i)+sig(i,3))
590 . + d4(i)*(sold4(i)+sig(i,4))
591 . + d5(i)*(sold5(i)+sig(i,5))
592 . + d6(i)*(sold6(i)+sig(i,6)))
593 eint(i)=eint(i)/vol(i)
594 dam(i,5)=kd1(i)*1000 + kd2(i)*100 + kd3(i)*10 + kd4(i) + 10000
595 ENDDO
596C
597 DO i=1,nel
598 sigym = max(em20,half*(one/f11(i)+one/f22(i)))
599 sigy(i)=sigmy(i)*sqrt(sigym)
600 defp(i)=pla(i)
601 ENDDO
602C
603 1000 FORMAT(1x,'RUPTURE OF SOLID ELEMENT NUMBER ',i10)
604C---
605 RETURN
606 END
m14ama
subroutine m14ama(pm, a, rx, ry, rz, sx, sy, sz, ax, ay, az, bx, by, bz, cx, cy, cz, nel, jcvt, jsph)
Definition m14ama.F:35
m14ftg
subroutine m14ftg(sig, ax, ay, az, bx, by, bz, cx, cy, cz, t1, t2, t3, t4, t5, t6, nel)
Definition m14ftg.F:32
m14gtf
subroutine m14gtf(sig, ax, ay, az, bx, by, bz, cx, cy, cz, d1, d2, d3, d4, d5, d6, t1, t2, t3, t4, t5, t6, e1, e2, e3, e4, e5, e6, nel)
Definition m14gtf.F:34
m14law
subroutine m14law(pm, off, sig, eint, pla, sigf, epsf, dam, epe, epc, a, vol, rx, ry, rz, sx, sy, sz, mat, vnew, dvol, ssp, d1, d2, d3, d4, d5, d6, sold1, sold2, sold3, sold4, sold5, sold6, sigy, defp, ngl, seq_output, nel, tsaiwu, jcvt, jsph)
Definition m14law.F:44
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21