OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
m25law.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!|| m25law ../engine/source/materials/mat/mat025/m25law.F
25!||--- called by ------------------------------------------------------
26!|| mmain ../engine/source/materials/mat_share/mmain.F90
27!||--- calls -----------------------------------------------------
28!|| mat25_crasurv_s ../engine/source/materials/mat/mat025/mat25_crasurv_s.F90
29!|| mat25_tsaiwu_s ../engine/source/materials/mat/mat025/mat25_tsaiwu_s.F90
30!|| mreploc ../engine/source/materials/mat_share/mreploc.F
31!|| mrotens ../engine/source/materials/mat_share/mrotens.f
32!||--- uses -----------------------------------------------------
33!|| mat25_crasurv_s_mod ../engine/source/materials/mat/mat025/mat25_crasurv_s.f90
34!|| mat25_tsaiwu_s_mod ../engine/source/materials/mat/mat025/mat25_tsaiwu_s.F90
35!|| matparam_def_mod ../common_source/modules/mat_elem/matparam_def_mod.F90
36!||====================================================================
37 SUBROUTINE m25law(MAT_PARAM,
38 1 PM, OFF, SIG,
39 2 EINT, S01, S02, S03,
40 3 S04, S05, S06, D1,
41 4 D2, D3, D4, D5,
42 5 D6, RX, RY, RZ,
43 6 SX, SY, SZ, TX,
44 7 TY, TZ, GAMA, VNEW,
45 8 SSP, VOL, EPSD, WPLA,
46 9 EPST, SIGL, TSAIWU,
47 A FLAY, NGL, NEL, NFT,
48 B ILAY, NPT, IPG,
49 C JCVT, JSPH, ISORTH, DMG,
50 D L_DMG, OUTV )
51!-----------------------------------------------
52! M o d u l e s
53!-----------------------------------------------
54 use matparam_def_mod
55 use mat25_tsaiwu_s_mod
56 use mat25_crasurv_s_mod
57C-----------------------------------------------
58C I m p l i c i t T y p e s
59C-----------------------------------------------
60#include "implicit_f.inc"
61#include "comlock.inc"
62C-----------------------------------------------
63C G l o b a l P a r a m e t e r s
64C-----------------------------------------------
65#include "mvsiz_p.inc"
66C-----------------------------------------------
67C C o m m o n B l o c k s
68C-----------------------------------------------
69#include "units_c.inc"
70#include "scr17_c.inc"
71#include "param_c.inc"
72#include "com08_c.inc"
73#include "impl1_c.inc"
74C-----------------------------------------------
75C D u m m y A r g u m e n t s
76C-----------------------------------------------
77 INTEGER, INTENT(IN) :: NFT
78 INTEGER, INTENT(IN) :: NPT
79 INTEGER, INTENT(IN) :: JCVT
80 INTEGER, INTENT(IN) :: JSPH
81 INTEGER, INTENT(IN) :: ISORTH
82 INTEGER, INTENT(IN) :: L_DMG
83 INTEGER NGL(MVSIZ),ILAY,NEL,IPG
84 INTEGER, INTENT(INOUT) :: OUTV(NEL)
85 my_real ,DIMENSION(NEL) ,INTENT(INOUT) :: EPSD ! local equivalent strain rate
86 my_real
87 . PM(NPROPM),OFF(*), SIG(NEL,6), WPLA(*), GAMA(MVSIZ,6),
88 . EINT(*),RX(*) ,RY(*),RZ(*) ,SX(*),SY(*),SZ(*),TX(*) ,TY(*) ,TZ(*)
89 my_real
90 . d1(*), d2(*), d3(*), d4(*),
91 . d5(*), d6(*),
92 . s01(*),s02(*),s03(*),s04(*),s05(*),s06(*),
93 . vol(*),vnew(*),ssp(*),epst(nel,6),sigl(mvsiz,6),flay(*),
94 . tsaiwu(nel)
95 type (matparam_struct_) ,intent(in) :: mat_param
96 my_real, DIMENSION(NEL,L_DMG) :: dmg
97C-----------------------------------------------
98C L o c a l V a r i a b l e s
99C-----------------------------------------------
100 INTEGER NFIS1(MVSIZ), NFIS2(MVSIZ), NFIS3(MVSIZ), INDX(MVSIZ)
101 INTEGER I, J, JADR, J1, J2, J3, JJ, II, I2, NINDX, I1,
102 . ioff,joff,idir,iflag
103C REAL
104 my_real
105 . degmb(mvsiz),degfx(mvsiz),
106 . wplar(mvsiz),strn1(mvsiz),strn2(mvsiz),strn3(mvsiz),
107 . damcr(mvsiz,2), dmaxt(mvsiz),strp1(mvsiz) ,strp2(mvsiz),
108 . epspl(mvsiz),s1(mvsiz),s2(mvsiz),s3(mvsiz),
109 . s4(mvsiz),s5(mvsiz),s6(mvsiz),
110 . ep1(mvsiz),ep2(mvsiz),ep3(mvsiz),ep4(mvsiz),ep5(mvsiz),
111 . ep6(mvsiz),epst1(mvsiz),epst2(mvsiz),epsm1(mvsiz),epsm2(mvsiz),
112 . dmax(mvsiz),r11(mvsiz),r12(mvsiz),r13(mvsiz),r21(mvsiz),
113 . r22(mvsiz),r23(mvsiz),r31(mvsiz),r32(mvsiz),r33(mvsiz)
114 my_real
115 . asrate,eps_k2, eps_m2, sigt1, sigt2,
116 . zt, wmc, visc, dtinv ,epd,dav,dam1,dam2,dt5,
117 . epst1_1,epst2_1,epsm1_1,epsm2_1,dmax_1,ssp_1
118C=======================================================================
119C Formulation flag
120
121 iflag = mat_param%IPARAM(1) ! IFLAG=0 : Tsai-Wu , IFLAG=1 : Crasurv
122
123 epst1_1 = mat_param%UPARAM(12) ! PM(60)
124 epst2_1 = mat_param%UPARAM(13) ! PM(61)
125 epsm1_1 = mat_param%UPARAM(14) ! PM(62)
126 epsm2_1 = mat_param%UPARAM(15) ! pm(63)
127 dmax_1 = mat_param%UPARAM(18) ! PM(64)
128
129 asrate = pm(9)*dt1
130 ssp_1 = pm(27)
131C-----------------------------------------------------------
132 DO i=1,nel
133 wplar(i) = zero
134 nfis1(i) = 0
135 nfis2(i) = 0
136 nfis3(i) = 0
137 dmaxt(i) = zero
138 ENDDO
139C-----------------------------------------------------------
140C STRAIN RATE FILTERING (EQUIVALENT EPSD)
141C-----------------------------------------------------------
142 asrate = asrate / (asrate + one)
143 DO i=1,nel
144 dav = -third*(d1(i)+d2(i)+d3(i))
145 epd = half*((d1(i)+dav )**2 + (d2(i)+dav)**2 + (d3(i)+dav)**2)
146 . + fourth*(d4(i)**2 + d5(i)**2 + d6(i)**2)
147 epd = sqrt(three*epd)/three_half
148 epsd(i) = asrate*epd + (one - asrate)*epsd(i)
149 ENDDO
150C--------------------------------------------
151C STRESS TRANSFORMATION (GLOBAL -> FIBER)
152C--------------------------------------------
153 DO i=1,nel
154 ep1(i) = d1(i)*dt1
155 ep2(i) = d2(i)*dt1
156 ep3(i) = d3(i)*dt1
157 ep4(i) = d4(i)*dt1
158 ep5(i) = d5(i)*dt1
159 ep6(i) = d6(i)*dt1
160 s1(i) = sig(i,1)
161 s2(i) = sig(i,2)
162 s3(i) = sig(i,3)
163 s4(i) = sig(i,4)
164 s5(i) = sig(i,5)
165 s6(i) = sig(i,6)
166 ENDDO
167!
168 IF (isorth > 0 .AND. jcvt == 0) THEN
169 CALL mreploc(
170 1 gama, r11, r12, r13,
171 2 r21, r22, r23, r31,
172 3 r32, r33, rx, ry,
173 4 rz, sx, sy, sz,
174 5 tx, ty, tz, nel, jsph)
175 DO i=1,nel
176 ep4(i) = half*ep4(i)
177 ep5(i) = half*ep5(i)
178 ep6(i) = half*ep6(i)
179 ENDDO
180 CALL mrotens(1,nel,ep1,ep2,ep3,ep4,ep5,ep6,
181 . r11,r12,r13,
182 . r21,r22,r23,
183 . r31,r32,r33)
184 DO i=1,nel
185 ep4(i) = two*ep4(i)
186 ep5(i) = two*ep5(i)
187 ep6(i) = two*ep6(i)
188 ENDDO
189C stress
190 CALL mrotens(1,nel,s1,s2,s3,s4,s5,s6,
191 . r11,r12,r13,
192 . r21,r22,r23,
193 . r31,r32,r33)
194 ENDIF
195C-----------------------------------------------------------
196C CONTRAINTES PLASTIQUEMENT ADMISSIBLES
197C-----------------------------------------------------------
198 IF (iflag == 0) THEN ! Tsai-Wu
199 CALL mat25_tsaiwu_s(mat_param,
200 1 nel ,ngl ,off ,flay ,
201 2 s1 ,s2 ,s3 ,s4 ,s5 ,s6 ,
202 3 ep1 ,ep2 ,ep3 ,ep4 , ep5 ,ep6 ,
203 4 epst ,nfis1,nfis2 ,nfis3,
204 5 wplar ,epsd ,wpla ,sigl ,ilay ,ipg ,
205 6 tsaiwu ,tt ,imconv,mvsiz,iout ,dmg ,
206 7 l_dmg ,outv)
207
208 ELSE ! Crasurv
209
210 CALL mat25_crasurv_s(mat_param ,
211 1 nel ,ngl, off ,flay ,
212 2 s1 ,s2 ,s3 ,s4 ,s5 ,s6 ,
213 3 ep1 ,ep2 ,ep3 ,ep4 , ep5 ,ep6,
214 4 epst ,nfis1,nfis2,nfis3,
215 5 wplar,epsd, wpla, sigl, ilay,
216 6 ipg ,tsaiwu,tt,imconv ,mvsiz ,iout,
217 7 dmg ,l_dmg ,outv)
218
219 END IF
220!-----------------------------------------------------------------------
221 IF (isorth > 0 .AND. jcvt == 0) THEN
222 CALL mrotens(1,nel,s1 ,s2 ,s3 ,s4 ,s5 ,s6 ,
223 . r11,r21,r31,r12,r22,r32,r13,r23,r33 )
224 ENDIF
225
226 DO i=1,nel
227 sig(i,1) = s1(i)*off(i)
228 sig(i,2) = s2(i)*off(i)
229 sig(i,3) = s3(i)*off(i)
230 sig(i,4) = s4(i)*off(i)
231 sig(i,5) = s5(i)*off(i)
232 sig(i,6) = s6(i)*off(i)
233 ENDDO
234C-----------------------
235C TENSILE RUPTURE
236C-----------------------
237 DO i=1,nel
238 epst1(i) =epst1_1
239 epst2(i) =epst2_1
240 epsm1(i) =epsm1_1
241 epsm2(i) =epsm2_1
242 dmax(i) =dmax_1
243 ssp(i) = ssp_1
244 ENDDO
245C.....STRAINS IN ORTHOTROPIC DIRECTIONS
246C
247C.....GATHER DIRECTION 1
248 nindx=0
249 DO i=1,nel
250 IF(epst(i,1) >= epst1(i)
251 + .AND. dmg(i,2) == zero .AND. off(i) == one)THEN
252 nindx=nindx+1
253 indx(nindx)=i
254 ENDIF
255 ENDDO
256C.....1.FIRST FAILURE DIRECTION 1
257 IF(nindx>0)THEN
258 idir=1
259 DO j=1,nindx
260 i=indx(j)
261 dam1=(epst(i,1)-epst1(i))/(epsm1(i)-epst1(i))
262 dam2= dam1*epsm1(i)/epst(i,1)
263 dmg(i,2)= min(dam2,dmax(i))
264 IF(dmg(i,2)==dmax(i).AND.imconv==1)THEN
265#include "lockon.inc"
266 WRITE(iout, '(A,I1,A,I10,A,I3,A,I3,A,1PE11.4)')
267 + ' FAILURE-',idir,', ELEMENT #',ngl(i),
268 + ', LAYER #',ilay,', INTEGRATION POINT #',ipg,
269 + ', TIME=',tt
270#include "lockoff.inc"
271 ENDIF
272 ENDDO
273 ENDIF
274C.....GATHER DIRECTION 2
275 nindx=0
276 DO i=1,nel
277 IF(epst(i,2) >= epst2(i)
278 + .AND. dmg(i,3) == zero .AND. off(i) == one) THEN
279 nindx=nindx+1
280 indx(nindx)=i
281 ENDIF
282 ENDDO
283C.....1.FIRST FAILURE DIRECTION 2
284 IF(nindx > 0)THEN
285 idir=2
286 DO j=1,nindx
287 i=indx(j)
288 dam1=(epst(i,2)-epst2(i))/(epsm2(i)-epst2(i))
289 dam2= dam1*epsm2(i)/epst(i,2)
290 dmg(i,3)= min(dam2,dmax(i))
291 IF(dmg(i,3) == dmax(i) .AND. imconv == 1)THEN
292#include "lockon.inc"
293 WRITE(iout, '(A,I1,A,I10,A,I3,A,I3,A,1PE11.4)')
294 + ' FAILURE-',idir,', ELEMENT #',ngl(i),
295 + ', LAYER #',ilay,', INTEGRATION POINT #',ipg,
296 + ', TIME=',tt
297#include "lockoff.inc"
298 ENDIF
299 ENDDO
300 ENDIF
301C---------------------------
302C GLOBAL FAILURE INDEX
303C---------------------------
304 IF (iflag == 0) THEN ! Tsai-Wu
305 DO i=1,nel
306 dmg(i,1) = max(dmg(i,2),dmg(i,3),dmg(i,4))
307 ENDDO
308 ELSE ! Crasurv
309 DO i=1,nel
310 dmg(i,1) = max(dmg(i,2),dmg(i,3),dmg(i,4))
311 IF (abs(dmg(i,5)) >= one) THEN
312 dmg(i,1) = max(abs(dmg(i,5))-one,dmg(i,1))
313 ENDIF
314 IF (abs(dmg(i,6)) >= one) THEN
315 dmg(i,1) = max(abs(dmg(i,6))-one,dmg(i,1))
316 ENDIF
317 IF (dmg(i,7) >= one) THEN
318 dmg(i,1) = max(dmg(i,7)-one,dmg(i,1))
319 ENDIF
320 ENDDO
321 ENDIF
322C----------------------------
323C TEST DE RUPTURE
324C----------------------------
325 DO i=1,nel
326 IF(off(i) < em01) off(i)=zero
327 IF(off(i) < one ) off(i)=off(i)*four_over_5
328 ENDDO
329
330 nindx=0
331 DO i=1,nel
332 IF(off(i)==one) THEN
333 ioff = mat_param%iparam(2)
334 IF(ioff < 0) ioff=-(ioff+1)
335 joff=0
336c IF(IOFF == 0 .AND. WPLAR(I) >= ONE) JOFF=1
337c IF(IOFF == 1 .AND. NINT(WPLAR(I)) >= NPT) JOFF=1
338c IF(IOFF == 2 .AND. NFIS1(I) == NPT) JOFF=1
339c IF(IOFF == 3 .AND. NFIS2(I) == NPT) JOFF=1
340c IF(IOFF == 4 .AND. NFIS1(I) == NPT
341c . .AND. NFIS2(I) == NPT) JOFF=1
342c IF(IOFF == 5 .AND. NFIS1(I) == NPT) JOFF=1
343c IF(IOFF == 5 .AND. NFIS2(I) == NPT) JOFF=1
344c IF(IOFF == 6 .AND. NFIS3(I) == NPT) JOFF=1
345 IF(joff == one) THEN
346 off(i)=off(i)*four_over_5
347 ii=i+nft
348 nindx=nindx+1
349 indx(nindx)=i
350 IF(imconv==1)THEN
351#include "lockon.inc"
352 WRITE(iout,1000) ngl(i)
353 WRITE(istdo,1100) ngl(i),tt
354#include "lockoff.inc"
355 ENDIF
356 ENDIF
357 ENDIF
358 ENDDO
359 dt5=half*dt1
360 DO i=1,nel
361 eint(i)=eint(i) + dt5*vnew(i)*
362 . ( d1(i)*(s01(i)+sig(i,1))
363 . + d2(i)*(s02(i)+sig(i,2))
364 . + d3(i)*(s03(i)+sig(i,3))
365 . + d4(i)*(s04(i)+sig(i,4))
366 . + d5(i)*(s05(i)+sig(i,5))
367 . + d6(i)*(s06(i)+sig(i,6)))
368 eint(i)=eint(i)/vol(i)
369 ENDDO
370c----------------------------------------------------------------------
371 1000 FORMAT(1x,'-- RUPTURE OF SOLID ELEMENT NUMBER ',i10)
372 1100 FORMAT(1x,'-- RUPTURE OF SOLID ELEMENT :',i10,' AT TIME :',g11.4)
373c-----------
374 RETURN
375 END
m25law
subroutine m25law(mat_param, pm, off, sig, eint, s01, s02, s03, s04, s05, s06, d1, d2, d3, d4, d5, d6, rx, ry, rz, sx, sy, sz, tx, ty, tz, gama, vnew, ssp, vol, epsd, wpla, epst, sigl, tsaiwu, flay, ngl, nel, nft, ilay, npt, ipg, jcvt, jsph, isorth, dmg, l_dmg, outv)
Definition m25law.F:51
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
mreploc
subroutine mreploc(ang, r11, r12, r13, r21, r22, r23, r31, r32, r33, rx, ry, rz, sx, sy, sz, tx, ty, tz, nel, jsph)
Definition mreploc.F:39
mrotens
subroutine mrotens(lft, llt, e1, e2, e3, e4, dyz, e6, r11, r12, r13, r21, r22, r23, r31, r32, r33)
Definition mrotens.F:36