OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
m12law.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!|| m12law ../engine/source/materials/mat/mat012/m12law.F
25!||--- called by ------------------------------------------------------
26!|| mmain ../engine/source/materials/mat_share/mmain.F90
27!||--- calls -----------------------------------------------------
28!|| eosmain ../common_source/eos/eosmain.f
29!|| m14ama ../engine/source/materials/mat/mat014/m14ama.F
30!|| m14ftg ../engine/source/materials/mat/mat014/m14ftg.f
31!|| m14gtf ../engine/source/materials/mat/mat014/m14gtf.F
32!||--- uses -----------------------------------------------------
33!|| eosmain_mod ../common_source/eos/eosmain.f
34!|| matparam_def_mod ../common_source/modules/mat_elem/matparam_def_mod.F90
35!||====================================================================
36 SUBROUTINE m12law(
37 1 PM, OFF, SIG, EINT,
38 2 PLA, SIGF, EPSF, DAM,
39 3 EPE, EPC, A, VOL,
40 4 RX, RY, RZ, SX,
41 5 SY, SZ, MAT, VNEW,
42 6 DVOL, SSP, D1, D2,
43 7 D3, D4, D5, D6,
44 8 SOLD1, SOLD2, SOLD3, SOLD4,
45 9 SOLD5, SOLD6, SIGY, DEFP,
46 A NGL, SEQ_OUTPUT, NEL, EOSTYP,
47 B RHO0, AMU, AMU2, ESPE,
48 C DF, PSH, PNEW, DPDM,
49 D DPDE, RHO, TEMP,
50 E BUFMAT, NPC, TF, TSAIWU,
51 F VAREOS, NVAREOS, JCVT, JSPH,
52 G MAT_PARAM, NVARTMP_EOS,VARTMP_EOS)
53C-----------------------------------------------
54C M o d u l e s
55C-----------------------------------------------
56 USE matparam_def_mod , ONLY : matparam_struct_
57 USE eosmain_mod , ONLY : eosmain
58C-----------------------------------------------
59C I m p l i c i t T y p e s
60C-----------------------------------------------
61#include "implicit_f.inc"
62#include "comlock.inc"
63C-----------------------------------------------
64C G l o b a l P a r a m e t e r s
65C-----------------------------------------------
66#include "mvsiz_p.inc"
67C-----------------------------------------------
68C C o m m o n B l o c k s
69C-----------------------------------------------
70#include "com01_c.inc"
71#include "param_c.inc"
72#include "com08_c.inc"
73#include "scr17_c.inc"
74#include "units_c.inc"
75 COMMON /tablesizf/ stf,snpc
76 INTEGER STF,SNPC
77C-----------------------------------------------
78C D u m m y A r g u m e n t s
79C-----------------------------------------------
80 TYPE(MATPARAM_STRUCT_), INTENT(IN) :: MAT_PARAM !material data structure
81 INTEGER, INTENT(IN) :: JCVT
82 INTEGER, INTENT(IN) :: JSPH
83 INTEGER MAT(MVSIZ),NGL(MVSIZ),NEL,EOSTYP
84 INTEGER,INTENT(IN) :: NVAREOS
85 my_real
86 . PM(NPROPM,*), OFF(*), SIG(NEL,6), EINT(*), PLA(*), SIGF(*),
87 . RX(*), RY(*), RZ(*), SX(*), SY(*), SZ(*),
88 . EPSF(*), DAM(NEL,5), EPE(NEL,3), EPC(NEL,3), A(MVSIZ,6), VOL(*),
89 . VNEW(MVSIZ), DVOL(MVSIZ), SSP(MVSIZ),
90 . D1(MVSIZ), D2(MVSIZ), D3(MVSIZ),D4(MVSIZ),D5(MVSIZ),D6(MVSIZ),
91 . sold1(mvsiz), sold2(mvsiz), sold3(mvsiz), sold4(mvsiz),
92 . sold5(mvsiz), sold6(mvsiz),sigy(*),defp(*),seq_output(*),
93 . tsaiwu(mvsiz),vareos(nvareos*nel)
94 my_real
95 . rho0(nel), amu(nel) , amu2(nel), espe(nel), df(nel) ,
96 . psh(nel) , pnew(nel) , dpdm(nel), dpde(nel), rho(nel),
97 . temp(nel), bufmat(*),muold(nel)
98 INTEGER,INTENT(IN)::NPC(SNPC)
99 my_real,INTENT(IN)::tf(stf)
100 INTEGER,INTENT(IN) :: NVARTMP_EOS
101 INTEGER,INTENT(INOUt) :: VARTMP_EOS(NEL,NVARTMP_EOS)
102C-----------------------------------------------
103C L o c a l V a r i a b l e s
104C-----------------------------------------------
105 INTEGER KD1(MVSIZ), KD2(MVSIZ), KD3(MVSIZ),
106 . KD4(MVSIZ),ICC(MVSIZ),
107 . I, IDAM, KDX, MX, NINDX, INDEX(MVSIZ), J, ICC_1
108 my_real BID(1) !bfrac argument not used in this context
109 my_real
110 . ax(mvsiz), ay(mvsiz), az(mvsiz), bx(mvsiz), by(mvsiz),
111 . bz(mvsiz), cx(mvsiz), cy(mvsiz),cz(mvsiz),
112 . wvec(mvsiz),s1(mvsiz), s2(mvsiz), s3(mvsiz),
113 . t1(mvsiz), t2(mvsiz), t3(mvsiz),t4(mvsiz),t5(mvsiz),t6(mvsiz),
114 . e1(mvsiz), e2(mvsiz), e3(mvsiz),e4(mvsiz),e5(mvsiz),e6(mvsiz),
115 . hard(mvsiz), sigmy(mvsiz), alpha(mvsiz), efib(mvsiz),
116 . epsft(mvsiz), epsfc(mvsiz), sigeq(mvsiz), p(mvsiz),
117 . d11(mvsiz), d12(mvsiz), d13(mvsiz), d22(mvsiz),
118 . d23(mvsiz), d33(mvsiz), g12(mvsiz), g23(mvsiz), g31(mvsiz),
119 . c11(mvsiz), c22(mvsiz), c33(mvsiz), c12(mvsiz), c23(mvsiz),
120 . c13(mvsiz), f11(mvsiz), f22(mvsiz), f44(mvsiz), f55(mvsiz),
121 . f12(mvsiz), f23(mvsiz), f1(mvsiz), f2(mvsiz), f4(mvsiz),
122 . f5(mvsiz), delta(mvsiz), so1(mvsiz),
123 . so2(mvsiz), so3(mvsiz), so4(mvsiz), so5(mvsiz), so6(mvsiz),
124 . ds1(mvsiz), ds2(mvsiz), ds3(mvsiz), ds4(mvsiz), ds5(mvsiz),
125 . ds6(mvsiz), dp1(mvsiz), dp2(mvsiz), dp3(mvsiz), dp4(mvsiz),
126 . dp5(mvsiz), dp6(mvsiz), lamda(mvsiz), coef(mvsiz), plas(mvsiz),
127 . cn(mvsiz), cb(mvsiz), cnn(mvsiz), sigt1(mvsiz), sigt2(mvsiz),
128 . sigt3(mvsiz),cc(mvsiz),epdr(mvsiz), f3(mvsiz), f33(mvsiz),
129 . f6(mvsiz),f66(mvsiz),f13(mvsiz)
130 my_real
131 . deve1(mvsiz), deve2(mvsiz), deve3(mvsiz), dav(mvsiz), pold(mvsiz),
132 . ener(mvsiz)
133 my_real
134 . fib, ca, sigmx, epsp, dt5 ,
135 . sigym,ca_1, cb_1, cn_1, cc_1,
136 . epdr_1,alpha_1,f1_1,f2_1,f3_1,
137 . f4_1,f5_1,f6_1,f11_1,f22_1,
138 . f33_1,f44_1,f55_1,f66_1,f12_1,
139 . f23_1,f13_1,c11_1,c22_1,c33_1,
140 . c12_1,c23_1,c13_1,d11_1,d12_1,
141 . d13_1,d22_1,d23_1,d33_1,g12_1,
142 . g23_1,g31_1,wplaref,dwpla
143C=======================================================================
144
145 IF(ncycle==0 .AND. irun==1) THEN
146 DO i=1,nel
147 kd1(i) = 0
148 kd2(i) = 0
149 kd3(i) = 0
150 kd4(i) = 0
151 ENDDO
152 ELSE
153 DO i=1,nel
154 idam=int(dam(i,5))-10000
155 kd1(i) = idam/1000
156 kdx = idam - kd1(i)*1000
157 kd2(i) = kdx/100
158 kdx = kdx - kd2(i)*100
159 kd3(i) = kdx/10
160 kd4(i) = kdx - kd3(i)*10
161
162 ENDDO
163 ENDIF
164C------------------
165C FIBER CONTENT
166C------------------
167 fib=zero
168 mx =mat(1)
169 alpha_1 =pm(39,mx)
170 DO i=1,nel
171 alpha(i)=alpha_1
172 fib=fib+alpha_1
173 ENDDO
174C
175C--------------------------------------------
176C STRESS TRANSFORMATION (GLOBAL -> FIBER)
177C--------------------------------------------
178 CALL m14ama(
179 1 pm, a, rx, ry,
180 2 rz, sx, sy, sz,
181 3 ax, ay, az, bx,
182 4 by, bz, cx, cy,
183 5 cz, nel, jcvt, jsph)
184 CALL m14gtf(sig,ax ,ay ,az ,bx ,by,
185 2 bz ,cx ,cy ,cz ,d1 ,d2,
186 3 d3 ,d4 ,d5 ,d6 ,t1 ,t2,
187 4 t3 ,t4 ,t5 ,t6 ,e1 ,e2,
188 5 e3 ,e4 ,e5 ,e6 ,nel)
189C---------------------------------
190C MATERIAL PROPERTIES
191C---------------------------------
192C
193 nindx=0
194 mx = mat(1)
195 ca_1 = pm(25,mx)
196 cb_1 = pm(26,mx)
197 cn_1 = pm(27,mx)
198 cc_1 = pm(50,mx)
199 epdr_1 = pm(51,mx)
200 icc_1 = nint(pm(52,mx))
201 DO i=1,nel
202 ca = ca_1
203 cb(i) = cb_1
204 cn(i) = cn_1
205 cc(i) = cc_1
206 epdr(i) = epdr_1
207 icc(i) = icc_1
208 epsp = max(abs(e1(i)),abs(e2(i)),abs(e3(i)),
209 . half*abs(e4(i)),half*abs(e5(i)),half*abs(e6(i)))
210 IF(epsp>epdr(i).AND.cc(i)/=zero) THEN
211 epsp = one + cc(i) * log(epsp/epdr(i))
212 ELSE
213 epsp = one
214 ENDIF
215 IF(icc(i)==1)THEN
216 sigmx =pm(28,mx)*epsp
217 ELSEIF(icc(i)==2)THEN
218 sigmx =pm(28,mx)
219 ELSEIF(icc(i)==3)THEN
220 sigmx =pm(28,mx)*epsp
221 ELSEIF(icc(i)==4)THEN
222 sigmx =pm(28,mx)
223 ENDIF
224 cb(i) =cb(i)*epsp
225 ca =ca*epsp
226 sigmy(i)= min(sigmx,ca+cb(i)*pla(i)**cn(i))
227 IF(sigmy(i)==sigmx .AND. off(i)==one) THEN
228 off(i)=zep99
229 kd4(i)=2
230C
231 nindx=nindx+1
232 index(nindx)=i
233 ENDIF
234 sigt1(i)= pm(81,mx)
235 sigt2(i)= pm(82,mx)
236 sigt3(i)= pm(83,mx)
237 ssp(i) = pm(49,mx)
238 delta(i)= pm(79,mx)
239 ENDDO
240C
241 IF(nindx/=0)THEN
242 DO j=1,nindx
243 i=index(j)
244#include "lockon.inc"
245 WRITE(iout,1000) ngl(i)
246#include "lockoff.inc"
247 END DO
248 END IF
249C
250 DO i=1,nel
251 e1(i)=e1(i)*dt1
252 e2(i)=e2(i)*dt1
253 e3(i)=e3(i)*dt1
254 e4(i)=e4(i)*dt1
255 e5(i)=e5(i)*dt1
256 e6(i)=e6(i)*dt1
257 ENDDO
258C
259 DO i=1,nel
260 epe(i,1)=epe(i,1)+e1(i)
261 epe(i,2)=epe(i,2)+e2(i)
262 epe(i,3)=epe(i,3)+e3(i)
263 ENDDO
264C
265 mx =mat(1)
266 f1_1 =pm(53,mx)
267 f2_1 =pm(54,mx)
268 f3_1 =pm(55,mx)
269 f4_1 =pm(56,mx)
270 f5_1 =pm(57,mx)
271 f6_1 =pm(58,mx)
272 f11_1 =pm(59,mx)
273 f22_1 =pm(60,mx)
274 f33_1 =pm(61,mx)
275 f44_1 =pm(62,mx)
276 f55_1 =pm(63,mx)
277 f66_1 =pm(64,mx)
278 f12_1 =pm(65,mx)
279 f23_1 =pm(66,mx)
280 f13_1 =pm(67,mx)
281C
282 c11_1 =pm(73,mx)
283 c22_1 =pm(74,mx)
284 c33_1 =pm(75,mx)
285 c12_1 =pm(76,mx)
286 c23_1 =pm(77,mx)
287 c13_1 =pm(78,mx)
288C
289 d11_1 =pm(40,mx)
290 d12_1 =pm(41,mx)
291 d13_1 =pm(42,mx)
292 d22_1 =pm(43,mx)
293 d23_1 =pm(44,mx)
294 d33_1 =pm(45,mx)
295 g12_1 =pm(46,mx)
296 g23_1 =pm(47,mx)
297 g31_1 =pm(48,mx)
298C
299 DO i=1,nel
300C
301 f1(i) = f1_1
302 f2(i) = f2_1
303 f3(i) = f3_1
304 f4(i) = f4_1
305 f5(i) = f5_1
306 f6(i) = f6_1
307 f11(i) = f11_1
308 f22(i) = f22_1
309 f33(i) = f33_1
310 f44(i) = f44_1
311 f55(i) = f55_1
312 f66(i) = f66_1
313 f12(i) = f12_1
314 f23(i) = f23_1
315 f13(i) = f13_1
316C
317 c11(i) = c11_1
318 c22(i) = c22_1
319 c33(i) = c33_1
320 c12(i) = c12_1
321 c23(i) = c23_1
322 c13(i) = c13_1
323C
324 d11(i) = d11_1
325 d12(i) = d12_1
326 d13(i) = d13_1
327 d22(i) = d22_1
328 d23(i) = d23_1
329 d33(i) = d33_1
330 g12(i) = g12_1
331 g23(i) = g23_1
332 g31(i) = g31_1
333 ENDDO
334C-------------------------
335C NEW ELASTIC STRESSES
336C-------------------------
337 DO i=1,nel
338 so1(i)=t1(i)
339 so2(i)=t2(i)
340 so3(i)=t3(i)
341 so4(i)=t4(i)
342 so5(i)=t5(i)
343 so6(i)=t6(i)
344 ENDDO
345C
346 IF (eostyp == 18) THEN ! linear EOS by default
347 DO i=1,nel
348 t1(i)=t1(i)+d11(i)*e1(i)+d12(i)*e2(i)+d13(i)*e3(i)
349 t2(i)=t2(i)+d12(i)*e1(i)+d22(i)*e2(i)+d23(i)*e3(i)
350 t3(i)=t3(i)+d13(i)*e1(i)+d23(i)*e2(i)+d33(i)*e3(i)
351 t4(i)=t4(i)+g12(i)*e4(i)
352 t5(i)=t5(i)+g23(i)*e5(i)
353 t6(i)=t6(i)+g31(i)*e6(i)
354 ENDDO
355 ELSE
356 DO i=1,nel
357 dav(i) = third*(e1(i)+e2(i)+e3(i))
358 pold(i) = third*(t1(i)+t2(i)+t3(i))
359 ener(i) = eint(i)
360 ENDDO
361 DO i=1,nel
362 deve1(i) = e1(i) - dav(i)
363 deve2(i) = e2(i) - dav(i)
364 deve3(i) = e3(i) - dav(i)
365 ENDDO
366 DO i=1,nel
367 t1(i)=t1(i)+d11(i)*deve1(i)+d12(i)*deve2(i)+d13(i)*deve3(i)-pold(i)
368 t2(i)=t2(i)+d12(i)*deve1(i)+d22(i)*deve2(i)+d23(i)*deve3(i)-pold(i)
369 t3(i)=t3(i)+d13(i)*deve1(i)+d23(i)*deve2(i)+d33(i)*deve3(i)-pold(i)
370 ENDDO
371C
372 muold(1:nel) =( rho(1:nel)/rho0)-one
373 sig(1:nel,1:6) = zero
374 CALL eosmain(0 ,nel ,eostyp ,pm ,off , ener ,
375 2 rho ,rho0 ,amu ,amu2 ,espe ,
376 3 dvol ,df ,vnew ,mat ,psh ,
377 4 pnew ,dpdm ,dpde ,temp ,
378 5 bufmat ,sig ,muold ,12 ,
379 6 npc ,tf ,vareos ,nvareos, mat_param, bid, nvartmp_eos,vartmp_eos)
380 CALL eosmain(1 ,nel ,eostyp ,pm ,off , ener ,
381 2 rho ,rho0 ,amu ,amu2 ,espe ,
382 3 dvol ,df ,vnew ,mat ,psh ,
383 4 pnew ,dpdm ,dpde ,temp ,
384 5 bufmat ,sig ,muold ,12 ,
385 6 npc ,tf ,vareos ,nvareos, mat_param, bid, nvartmp_eos,vartmp_eos)
386 DO i=1,nel
387 t1(i)=t1(i)-pnew(i)
388 t2(i)=t2(i)-pnew(i)
389 t3(i)=t3(i)-pnew(i)
390 ENDDO
391 DO i=1,nel
392 ssp(i) = sqrt(abs(dpdm(i))/rho0(i))
393 ENDDO
394
395 ENDIF
396C------------------
397C FIBERS STRESS
398C------------------
399 IF (fib>0) THEN
400 mx = mat(1)
401 DO i=1,nel
402 efib(i) =pm(38,mx)
403 epsft(i)=pm(84,mx)
404 epsfc(i)=pm(85,mx)
405 epsf(i) = epsf(i)+ e1(i)
406 sigf(i) = efib(i)*epsf(i)
407 ENDDO
408C
409C DEACTIVATED (DAM4 AND DAM5 TEMPORARILY USED FOR GENERAL DAMAGE INFO)
410C
411 ENDIF
412C-------------------------
413C GENERAL FAILURE TEST
414C-------------------------
415 DO i=1,nel
416 IF(off(i) < em01) off(i)=zero
417 IF(off(i) < one) off(i)=off(i)*four_over_5
418 ENDDO
419C-------------------
420C TENSION DAMAGE
421C-------------------
422 DO i=1,nel
423C------------------------
424C DIRECTION 1 (FIBER)
425C------------------------
426 wvec(i)=(one-dam(i,1))*sigt1(i)
427 IF(t1(i) > wvec(i)) THEN
428C
429 IF(epc(i,1) == zero) THEN
430 epc(i,1)= max(epe(i,1),zero)
431 ELSE
432 epc(i,1)= max(epc(i,1)+e1(i),zero)
433 ENDIF
434C
435 t1(i)=wvec(i)
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 IF(kd1(i)==0) THEN
439 kd1(i)=1
440 ENDIF
441 dam(i,1)= min((dam(i,1)+delta(i)),one)
442 IF(dam(i,1)>=one .AND. kd1(i)/=2) THEN
443 kd1(i)=2
444 ENDIF
445 ENDIF
446C
447 IF(e1(i) < zero .AND. dam(i,1) > zero)
448 . epc(i,1)= max(epc(i,1)+e1(i),zero)
449C----------------
450C DIRECTION 2
451C----------------
452 wvec(i)=(one-dam(i,2))*sigt2(i)
453 IF(t2(i) > wvec(i)) THEN
454C
455 IF(epc(i,2)==zero) THEN
456 epc(i,2)= max(epe(i,2),zero)
457 ELSE
458 epc(i,2)= max(epc(i,2)+e2(i),zero)
459 ENDIF
460C
461 t1(i)=t1(i)-d12(i)*e2(i)*dam(i,2)
462 t2(i)=wvec(i)
463 t3(i)=t3(i)-d23(i)*e2(i)*dam(i,2)
464 IF(kd2(i) == 0) THEN
465 kd2(i)=1
466 ENDIF
467 dam(i,2)= min((dam(i,2)+delta(i)),one)
468 IF(dam(i,2) >= one .AND. kd2(i) /= 2) THEN
469 kd2(i)=2
470 ENDIF
471 ENDIF
472C
473 IF(e2(i) < zero .AND. dam(i,2) > zero)
474 . epc(i,2)= max(epc(i,2)+e2(i),zero)
475C----------------
476C DIRECTION 3
477C----------------
478 wvec(i)=(one-dam(i,3))*sigt3(i)
479 IF(t3(i) > wvec(i)) THEN
480C
481 IF(epc(i,3)==zero) THEN
482 epc(i,3)= max(epe(i,3),zero)
483 ELSE
484 epc(i,3)= max(epc(i,3)+e3(i),zero)
485 ENDIF
486C
487 t1(i)=t1(i)-d13(i)*e3(i)*dam(i,3)
488 t2(i)=t2(i)-d23(i)*e3(i)*dam(i,3)
489 t3(i)=wvec(i)
490 IF(kd3(i)==0) THEN
491 kd3(i)=1
492 ENDIF
493 dam(i,3)= min((dam(i,3)+delta(i)),one)
494 IF(dam(i,3) >= one .AND. kd3(i) /= 2) THEN
495 kd3(i)=2
496 ENDIF
497 ENDIF
498C
499 IF(e3(i)<zero .AND. dam(i,3)>zero)
500 . epc(i,3)= max(epc(i,3)+e3(i),zero)
501C
502 ENDDO
503C-----------------------------------
504C CRACK OPEN --> NO COMPRESSION
505C-----------------------------------
506 DO i=1,nel
507 IF(t1(i) < zero.AND. epc(i,1) > zero) THEN
508 t1(i)=zero
509 t2(i)=t2(i)-d12(i)*e1(i)*dam(i,1)
510 t3(i)=t3(i)-d13(i)*e1(i)*dam(i,1)
511 ENDIF
512 IF(t2(i) < zero.AND. epc(i,2) > zero) THEN
513 t1(i)=t1(i)-d12(i)*e2(i)*dam(i,2)
514 t2(i)=zero
515 t3(i)=t3(i)-d23(i)*e2(i)*dam(i,2)
516 ENDIF
517 IF(t3(i) < zero.AND.epc(i,3) > zero) THEN
518 t3(i)=zero
519 t1(i)=t1(i)-d13(i)*e3(i)*dam(i,3)
520 t2(i)=t2(i)-d23(i)*e3(i)*dam(i,3)
521 ENDIF
522 ENDDO
523C---------------------
524C PLASTICITY START
525C---------------------
526 DO i=1,nel
527 wvec(i)= f1(i)*t1(i) + f2(i)*t2(i) + f3(i)*t3(i)
528 . + f11(i)*t1(i)*t1(i) + f22(i)*t2(i)*t2(i) + f33(i)*t3(i)*t3(i)
529 . + f55(i)*t5(i)*t5(i) + f44(i)*t4(i)*t4(i) + f66(i)*t6(i)*t6(i)
530 . + two*f12(i)*t1(i)*t2(i) + two*f13(i)*t1(i)*t3(i)
531 . + two*f23(i)*t2(i)*t3(i)
532 dam(i,4)=wvec(i)
533 tsaiwu(i)=max(min(wvec(i)/sigmy(i),one),tsaiwu(i))
534!! CRITERION / (EQUIVALENT STRESS) FOR OUTPUT - TSAI-WU -
535!! SEQ_OUTPUT(I) = WVEC(I)
536 ENDDO
537C
538 DO i=1,nel
539 coef(i)= zero
540 IF (wvec(i) > sigmy(i) .and .off(i) == one) THEN
541 coef(i)=one
542 IF(kd4(i)==0) THEN
543 kd4(i)=1
544 ENDIF
545 ENDIF
546 ENDDO
547C
548 DO i=1,nel
549 dp1(i)=f1(i) + two*f11(i)*so1(i) + two*f12(i)*so2(i)
550 . + two*f13(i)*so3(i)
551 dp2(i)=f2(i) + two*f22(i)*so2(i) + two*f12(i)*so1(i)
552 . + two*f23(i)*so3(i)
553 dp3(i)=f3(i) + two*f33(i)*so3(i) + two*f13(i)*so1(i)
554 . + two*f23(i)*so2(i)
555 dp4(i)=two*f44(i)*so4(i)
556 dp5(i)=two*f55(i)*so5(i)
557 dp6(i)=two*f66(i)*so6(i)
558 ENDDO
559
560 3101 FORMAT(/' 205 WVEC SIGMY T1-T6 ',2e11.4/6e11.4)
561 3102 FORMAT(' F1 F2 F11 F22 F12 F23 F44 F55',3e11.4/5e11.4)
562C
563 DO i=1,nel
564 ds1(i)=t1(i)-so1(i)
565 ds2(i)=t2(i)-so2(i)
566 ds3(i)=t3(i)-so3(i)
567 ds4(i)=t4(i)-so4(i)
568 ds5(i)=t5(i)-so5(i)
569 ds6(i)=t6(i)-so6(i)
570 ENDDO
571C
572 DO i=1,nel
573 lamda(i)=(dp1(i)*ds1(i)+dp2(i)*ds2(i)+dp3(i)*ds3(i)
574 . +dp4(i)*ds4(i)+dp5(i)*ds5(i)+dp6(i)*ds6(i))*coef(i)
575 ENDDO
576C ***********
577 3103 FORMAT(' 207 LAMDA DS1-DS6 ',e11.4/6e11.4)
578C
579 DO i=1,nel
580 cnn(i)=cn(i)-one
581 plas(i)= one
582 IF(pla(i) > zero) plas(i)=pla(i)**cnn(i)
583 ENDDO
584C
585 DO i=1,nel
586 IF(lamda(i)/= zero) THEN
587 lamda(i)=lamda(i)*coef(i)/
588 . (dp1(i)*(d11(i)*dp1(i)+d12(i)*dp2(i)+d13(i)*dp3(i))+
589 . dp2(i)*(d12(i)*dp1(i)+d22(i)*dp2(i)+d23(i)*dp3(i))+
590 . dp3(i)*(d13(i)*dp1(i)+d23(i)*dp2(i)+d33(i)*dp3(i))+
591 . two*dp4(i)*g12(i)*dp4(i)+
592 . two*dp5(i)*g23(i)*dp5(i)+
593 . two*dp6(i)*g31(i)*dp6(i)+
594 . (so1(i)*dp1(i)+so2(i)*dp2(i)+so3(i)*dp3(i)+
595 . two*so4(i)*dp4(i)+two*so5(i)*dp5(i)+two*so6(i)*dp6(i))
596 . *cn(i)*cb(i)*plas(i) )
597C
598 ENDIF
599 ENDDO
600c************
601 3104 FORMAT(' 208 LAMDA ',e11.4)
602C
603 DO i=1,nel
604 dp1(i)=lamda(i)*dp1(i)
605 dp2(i)=lamda(i)*dp2(i)
606 dp3(i)=lamda(i)*dp3(i)
607 dp4(i)=lamda(i)*dp4(i)
608 dp5(i)=lamda(i)*dp5(i)
609 dp6(i)=lamda(i)*dp6(i)
610 ENDDO
611 3105 FORMAT(' 209 PLA DP1-DP6',e11.4/6e11.4)
612C
613 DO i=1,nel
614 epe(i,1)=epe(i,1)-dp1(i)
615 epe(i,2)=epe(i,2)-dp2(i)
616 epe(i,3)=epe(i,3)-dp3(i)
617 ENDDO
618C
619 DO i=1,nel
620 t1(i)=t1(i)-d11(i)*dp1(i)-d12(i)*dp2(i)-d13(i)*dp3(i)
621 t2(i)=t2(i)-d12(i)*dp1(i)-d22(i)*dp2(i)-d23(i)*dp3(i)
622 t3(i)=t3(i)-d13(i)*dp1(i)-d23(i)*dp2(i)-d33(i)*dp3(i)
623 t4(i)=t4(i)-g12(i)*dp4(i)*two
624 t5(i)=t5(i)-g23(i)*dp5(i)*two
625 t6(i)=t6(i)-g31(i)*dp6(i)*two
626 ENDDO
627C
628 mx = mat(1)
629 wplaref = pm(68 ,mx)
630 DO i=1,nel
631 dwpla = half*
632 . (dp1(i)*(t1(i)+so1(i))+
633 . dp2(i)*(t2(i)+so2(i))+
634 . dp3(i)*(t3(i)+so3(i))+
635 . two*dp4(i)*(t4(i)+so4(i))+
636 . two*dp5(i)*(t5(i)+so5(i))+
637 . two*dp6(i)*(t6(i)+so6(i)))
638 dwpla = max(dwpla ,zero) / wplaref
639 pla(i) = pla(i) + dwpla
640 pla(i)= max(pla(i),zero)
641 ENDDO
642C-------------------
643C PLASTICITY END
644C-------------------
645C--------------------------------------------
646C STRESS TRANSFORMATION (FIBER -> GLOBAL)
647C--------------------------------------------
648 CALL m14ftg(sig,ax ,ay ,az ,bx ,by ,
649 2 bz ,cx ,cy ,cz ,t1 ,t2 ,
650 3 t3 ,t4 ,t5 ,t6 ,nel)
651C
652 DO i=1,nel
653 sig(i,1)=sig(i,1)*off(i)
654 sig(i,2)=sig(i,2)*off(i)
655 sig(i,3)=sig(i,3)*off(i)
656 sig(i,4)=sig(i,4)*off(i)
657 sig(i,5)=sig(i,5)*off(i)
658 sig(i,6)=sig(i,6)*off(i)
659 ENDDO
660C--------------------
661C INTERNAL ENERGY
662C--------------------
663 dt5=half*dt1
664 DO i=1,nel
665 eint(i)=eint(i)+dt5*vnew(i)*
666 . ( d1(i)*(sold1(i)+sig(i,1))
667 . + d2(i)*(sold2(i)+sig(i,2))
668 . + d3(i)*(sold3(i)+sig(i,3))
669 . + d4(i)*(sold4(i)+sig(i,4))
670 . + d5(i)*(sold5(i)+sig(i,5))
671 . + d6(i)*(sold6(i)+sig(i,6)))
672 eint(i)=eint(i)/vol(i)
673 dam(i,5)=kd1(i)*1000 + kd2(i)*100 + kd3(i)*10 + kd4(i) + 10000
674 ENDDO
675C
676 DO i=1,nel
677 sigym = max(em20,half*(one/f11(i)+one/f22(i)))
678 sigy(i)=sigmy(i)*sqrt(sigym)
679 defp(i)=pla(i)
680 ENDDO
681C------------------------------------------------------------------------
682 1000 FORMAT(1x,'RUPTURE OF SOLID ELEMENT NUMBER ',i10)
683 6001 FORMAT(i6,' FAIL(14) E(',i5,') M(',a3,') T(',6e10.3,')')
684 RETURN
685 END
alpha
#define alpha
Definition eval.h:35
m12law
subroutine m12law(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, eostyp, rho0, amu, amu2, espe, df, psh, pnew, dpdm, dpde, rho, temp, bufmat, npc, tf, tsaiwu, vareos, nvareos, jcvt, jsph, mat_param, nvartmp_eos, vartmp_eos)
Definition m12law.F:53
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
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
eosmain_mod
Definition eosmain.F:33
eosmain_mod::eosmain
subroutine eosmain(iflag, nel, eostyp, pm, off, eint, rho, rho0, mu, mu2, espe, dvol, df, vnew, mat, psh, pnew, dpdm, dpde, theta, bufmat, sig, mu_bak, mlw, npf, tf, vareos, nvareos, mat_param, bfrac, nvartmp, vartmp)
Definition eosmain.F:80