OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps25c.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!|| sigeps25c ../engine/source/materials/mat/mat025/sigeps25c.F
25!||--- called by ------------------------------------------------------
26!|| mulawc ../engine/source/materials/mat_share/mulawc.F90
27!||--- calls -----------------------------------------------------
28!|| m25crak ../engine/source/materials/mat/mat025/m25crak.F
29!|| mat25_crasurv_c ../engine/source/materials/mat/mat025/mat25_crasurv_c.F90
30!|| mat25_tsaiwu_c ../engine/source/materials/mat/mat025/mat25_tsaiwu_c.F90
31!|| prony25c ../engine/source/materials/mat/mat025/prony25c.F
32!||--- uses -----------------------------------------------------
33!|| mat25_crasurv_c_mod ../engine/source/materials/mat/mat025/mat25_crasurv_c.F90
34!|| mat25_tsaiwu_c_mod ../engine/source/materials/mat/mat025/mat25_tsaiwu_c.F90
35!|| matparam_def_mod ../common_source/modules/mat_elem/matparam_def_mod.F90
36!||====================================================================
37 SUBROUTINE sigeps25c(MAT_PARAM,
38 1 NEL ,PM ,OFF ,GSTR ,
39 2 DIR ,THLY ,TIME ,TIMESTEP,SHF ,
40 3 NGL ,THK0 ,EXX ,OFF_OLD ,
41 4 EYY ,EXY ,EXZ ,EYZ ,KXX ,
42 5 KYY ,KXY ,ZZ ,EPSD_PG ,RHO0 ,
43 6 SOUNDSP ,UPARAM ,OFFL ,EPSD ,ASRATE ,
44 7 SIGY ,ZCFAC ,NPTT ,ILAY ,
45 8 NFIS1 ,NFIS2 ,NFIS3 ,WPLAR ,
46 9 NPTTOT ,IGTYP ,SIGV ,SIGPLY ,
47 A SIGOXX ,SIGOYY ,SIGOXY ,SIGOYZ ,SIGOZX ,
48 B SIGNXX ,SIGNYY ,SIGNXY ,SIGNYZ ,SIGNZX ,
49 C SIGVXX ,SIGVYY ,SIGVXY ,SIGVYZ ,SIGVZX ,
50 D ISRATE ,UVARV ,ISHPLYXFEM,IPT ,SEQ_OUTP,
51 E PLY_EXX ,PLY_EYY ,PLY_EXY ,PLY_EXZ ,PLY_EYZ ,
52 F PLY_F ,PLA ,CRAK ,IERR ,
53 G IOFF_DUCT,IFAILURE,PLY_ID ,IPG ,TSAIWU ,
54 H IMCONV ,IOUT ,DMG ,L_DMG )
55C-----------------------------------------------
56C M o d u l e s
57C-----------------------------------------------
58 USE matparam_def_mod
59 use mat25_tsaiwu_c_mod
60 use mat25_crasurv_c_mod
61C-----------------------------------------------
62C I m p l i c i t T y p e s
63C-----------------------------------------------
64#include "implicit_f.inc"
65#include "comlock.inc"
66C-----------------------------------------------
67C G l o b a l P a r a m e t e r s
68C-----------------------------------------------
69#include "mvsiz_p.inc"
70C-----------------------------------------------
71C C o m m o n B l o c k s
72C-----------------------------------------------
73#include "param_c.inc"
74#include "com01_c.inc"
75C-----------------------------------------------
76C D u m m y A r g u m e n t s
77C-----------------------------------------------
78 INTEGER NEL,ISRATE,PLY_ID,
79 . IPT,NPTT,NPTTOT,ILAY,IGTYP,IFAILURE,ISHPLYXFEM,IPG
80 INTEGER ,INTENT(IN) :: IMCONV
81 INTEGER ,INTENT(IN) :: IOUT,L_DMG
82 INTEGER NGL(MVSIZ),IERR(NEL),IOFF_DUCT(MVSIZ)
83 INTEGER ,DIMENSION(NEL) ,INTENT(INOUT) :: NFIS1,NFIS2,NFIS3
84 my_real ,INTENT(IN) :: TIME
85 my_real ,INTENT(IN) :: TIMESTEP
86 my_real ,INTENT(IN) :: ASRATE
87 my_real ,DIMENSION(NEL) ,INTENT(IN) :: EPSD_PG
88 my_real ,DIMENSION(NEL) ,INTENT(INOUT) :: EPSD
89 my_real ,DIMENSION(NEL) ,INTENT(INOUT) :: WPLAR
90 my_real ::
91 . PM(NPROPM),OFF(*),OFF_OLD(*),GSTR(NEL,8),DIR(*),
92 . THLY(*),SHF(*),THK0(MVSIZ),EXX(MVSIZ),EYY(MVSIZ),
93 . EXY(MVSIZ),EXZ(MVSIZ),EYZ(MVSIZ),KXX(MVSIZ),KYY(MVSIZ),
94 . KXY(MVSIZ),ZZ(*),RHO0(*),SOUNDSP(*),UPARAM(*),
95 . sigy(*),zcfac(mvsiz,2),signxx(mvsiz),
96 . signyy(mvsiz),signxy(mvsiz),signyz(mvsiz),signzx(mvsiz),
97 . sigvxx(mvsiz),sigvyy(mvsiz),sigvxy(mvsiz),sigvyz(mvsiz),
98 . sigvzx(mvsiz),uvarv(nel,*),ply_f(mvsiz,5,*),
99 . ply_exx(mvsiz,*),ply_eyy(mvsiz,*),ply_exy(mvsiz,*),
100 . ply_exz(mvsiz,*),ply_eyz(mvsiz,*),
101 . sigoxx(nel),sigoyy(nel),sigoxy(nel),sigoyz(nel),
102 . sigozx(nel),pla(nel),crak(nel,2),seq_outp(nel),
103 . offl(nel),sigv(nel,5),sigply(nel,3),tsaiwu(nel)
104 TYPE(matparam_struct_) ,INTENT(IN) :: MAT_PARAM
105 my_real, DIMENSION(NEL,L_DMG), INTENT(INOUT) :: dmg
106C-----------------------------------------------
107C L o c a l V a r i a b l e s
108C-----------------------------------------------
109 INTEGER :: I,J,IFLAG,IOFF,JOFF,FAILNPT,NPRONY
110 my_real :: ZT,RATIO,KV,DTINV,EPST1,EPST2,EPSM1,EPSM2,DMAX
111 my_real ::
112 . EPS(MVSIZ,5),STRN1(MVSIZ),STRN2(MVSIZ),STRN3(MVSIZ),
113 . STRP1(MVSIZ),STRP2(MVSIZ),SIGE(MVSIZ,5),
114 . EPSPXX(MVSIZ),ETSE(MVSIZ),YLD(MVSIZ),
115 . EPSPYY(MVSIZ),EPSPXY(MVSIZ),EPSPYZ(MVSIZ),EPSPZX(MVSIZ),
116 . SIG(NEL,5),SIGPE(MVSIZ,5)
117 my_real,DIMENSION(:), ALLOCATABLE :: GV,BETA
118C=======================================================================
119! Formulation flag
120 IFLAG = mat_param%IPARAM(1) ! IFLAG=0 : Tsai-Wu , IFLAG=1 : Crasurv
121!
122 ioff = mat_param%IPARAM(2)
123
124 epst1 = mat_param%uparam(12) ! PM(60)
125 epst2 = mat_param%uparam(13) ! PM(61)
126 epsm1 = mat_param%uparam(14) ! PM(62)
127 epsm2 = mat_param%uparam(15) ! PM(63)
128 dmax = mat_param%uparam(18) ! PM(64)
129!
130 ratio = mat_param%uparam(19) ! PM(188)
131C-----------------------------------------------------------
132C Visco elastic model
133C-----------------------------------------------------------
134 IF (mat_param%IVISC == 1) THEN ! PRONY model
135 nprony = mat_param%VISC%IPARAM(1)
136 kv = mat_param%VISC%UPARAM(1)
137 ALLOCATE(gv(nprony),beta(nprony))
138 DO i=1,nprony
139 gv(i) = mat_param%VISC%UPARAM(1 + i)
140 beta(i) = mat_param%VISC%UPARAM(1 + nprony + i)
141 ENDDO
142 ELSE
143 ALLOCATE(gv(0),beta(0))
144 ENDIF
145C-----------------------------------------------------------
146 IF (ishplyxfem == 0 .OR. iplyxfem == 2) THEN
147 DO i=1,nel
148 zt=zz(i)*thk0(i)
149 eps(i,1)=exx(i)+zt*kxx(i)
150 eps(i,2)=eyy(i)+zt*kyy(i)
151 eps(i,3)=exy(i)+zt*kxy(i)
152 eps(i,4)=eyz(i)
153 eps(i,5)=exz(i)
154 strn1(i)= gstr(i,1)+zt*gstr(i,6)
155 strn2(i)= gstr(i,2)+zt*gstr(i,7)
156 strn3(i)=(gstr(i,3)+zt*gstr(i,8))*half
157 ENDDO
158 ELSE
159 DO i=1,nel
160 zt=zz(i)*thk0(i)
161 eps(i,1)=exx(i) + zt*kxx(i) + ply_exx(i,ipt)
162 eps(i,2)=eyy(i) + zt*kyy(i) + ply_eyy(i,ipt)
163 eps(i,3)=exy(i) + zt*kxy(i) + ply_exy(i,ipt)
164 eps(i,4)=eyz(i)
165 eps(i,5)=exz(i)
166 strn1(i)= gstr(i,1)+ zt*gstr(i,6)
167 strn2(i)= gstr(i,2)+ zt*gstr(i,7)
168 strn3(i)=(gstr(i,3)+ zt*gstr(i,8))*half
169 ENDDO
170 ENDIF
171C
172 IF (mat_param%IVISC == 1) THEN
173 dtinv = timestep / max(timestep**2,em20)
174 DO i=1,nel
175 epspxx(i) = eps(i,1)*dtinv
176 epspyy(i) = eps(i,2)*dtinv
177 epspxy(i) = eps(i,3)*dtinv
178 epspyz(i) = eps(i,4)*dtinv
179 epspzx(i) = eps(i,5)*dtinv
180 ENDDO
181 ENDIF
182C-----------------------
183 DO i=1,nel
184 sig(i,1)=sigoxx(i)
185 sig(i,2)=sigoyy(i)
186 sig(i,3)=sigoxy(i)
187 sig(i,4)=sigoyz(i)
188 sig(i,5)=sigozx(i)
189 ENDDO
190!-----------------------------------------------------------------------
191 IF (iflag == 0) THEN ! Tsai-Wu
192 CALL mat25_tsaiwu_c(mat_param,
193 1 nel ,off ,sig ,
194 2 pla ,dir ,crak ,
195 3 nfis1 ,nfis2 ,nfis3 ,ilay ,shf ,
196 4 ngl ,eps ,igtyp ,wplar ,strn1 ,
197 5 strn2 ,strn3 ,strp1 ,strp2 ,sige ,
198 6 epsd_pg ,epsd ,israte ,asrate ,offl ,
199 7 yld ,etse,ishplyxfem,ply_exx(1,ipt),ply_eyy(1,ipt),
200 8 ply_exy(1,ipt),sigply,sigpe,ply_id ,
201 9 signxx ,signyy ,signxy ,signyz ,signzx,
202 a ipg ,tsaiwu ,iplyxfem,time ,timestep,
203 b imconv ,mvsiz ,iout ,dmg ,l_dmg )
204
205 ELSE ! Crasurv
206
207 CALL mat25_crasurv_c(mat_param,
208 1 nel ,off ,sig ,
209 2 pla ,dir ,crak ,
210 3 nfis1 ,nfis2 ,nfis3 ,ilay ,shf ,
211 4 ngl ,eps ,wplar ,strn1 ,
212 5 strn2 ,strn3 ,strp1 ,strp2 ,sige ,
213 6 epsd_pg ,epsd ,israte ,asrate ,offl ,
214 7 yld ,etse ,ierr,ishplyxfem,ply_exx(1,ipt),ply_eyy(1,ipt),
215 8 ply_exy(1,ipt),sigply,sigpe,ply_id ,
216 9 signxx ,signyy ,signxy ,signyz ,signzx,
217 a ipg ,tsaiwu ,iplyxfem,time ,timestep,
218 b imconv ,mvsiz ,iout ,dmg ,l_dmg )
219
220 END IF
221!-----------------------------------------------------------------------
222 IF (iplyxfem == 1) THEN
223#include "vectorize.inc"
224 DO i=1,nel
225 ply_f(i,1,ipt) = thly(i)*sige(i,1)
226 ply_f(i,2,ipt) = thly(i)*sige(i,2)
227 ply_f(i,3,ipt) = thly(i)*sige(i,3)
228 ply_f(i,4,ipt) = thly(i)*sige(i,4)
229 ply_f(i,5,ipt) = thly(i)*sige(i,5)
230 ENDDO
231 ELSEIF (iplyxfem == 2) THEN
232#include "vectorize.inc"
233 DO i=1,nel
234 ply_f(i,1,ipt) = thly(i)*(sige(i,1)+sigpe(i,1))
235 ply_f(i,2,ipt) = thly(i)*(sige(i,2)+sigpe(i,2))
236 ply_f(i,3,ipt) = thly(i)*(sige(i,3)+sigpe(i,3))
237 ply_f(i,4,ipt) = thly(i)*sige(i,4)
238 ply_f(i,5,ipt) = thly(i)*sige(i,5)
239 ENDDO
240 ENDIF ! IF (IPLYXFEM == 1)
241C-----------------------------------------------------------
242C Visco elastic model
243C-----------------------------------------------------------
244 IF (mat_param%IVISC == 1 ) THEN
245 CALL prony25c(nel ,nprony ,beta ,kv ,
246 1 gv ,timestep,rho0 ,off ,dir ,
247 2 epspxx ,epspyy ,epspxy ,epspyz ,epspzx,
248 3 sigvxx ,sigvyy ,sigvxy ,sigvyz ,sigvzx,
249 4 sigv ,soundsp ,uvarv ,igtyp )
250C
251 IF (ishplyxfem > 0) THEN
252#include "vectorize.inc"
253 DO i=1,nel
254 ply_f(i,1,ipt) = ply_f(i,1,ipt) + thly(i)*sigvxx(i)
255 ply_f(i,2,ipt) = ply_f(i,2,ipt) + thly(i)*sigvyy(i)
256 ply_f(i,3,ipt) = ply_f(i,3,ipt) + thly(i)*sigvxy(i)
257 ply_f(i,4,ipt) = ply_f(i,4,ipt) + thly(i)*sigvyz(i)
258 ply_f(i,5,ipt) = ply_f(i,5,ipt) + thly(i)*sigvzx(i)
259 ENDDO
260 ENDIF ! IF (ISHPLYXFEM > 0)
261 ENDIF
262C-----------------------
263C For QEPH
264C-----------------------
265 sigy(1:nel) = sigy(1:nel) + yld(1:nel)/npttot ! NPTT = MAX(NPTT) for IGTYP=51
266 DO i=1,nel
267 zcfac(i,1) = zcfac(i,1) + etse(i) / npttot
268 zcfac(i,2) = min(etse(i),zcfac(i,2))
269 IF (offl(i)==zero) zcfac(i,2)=-abs(zcfac(i,2))
270 ENDDO
271C-----------------------
272C TENSILE RUPTURE
273C-----------------------
274 CALL m25crak(nel ,off ,dmg ,l_dmg ,dir ,ilay ,
275 . thly ,ngl ,strp1 ,strp2 ,ply_id ,igtyp ,
276 . ipg ,epst1 ,epst2 ,epsm1 ,epsm2 ,dmax )
277C-----------------------
278 DEALLOCATE(gv,beta)
279C---- -------------------
280 DO i=1,nel
281 IF (off(i) == zero .AND. nptt > 1) offl(i) = zero
282 ENDDO
283C---------------------------
284C GLOBAL FAILURE INDEX
285C---------------------------
286 IF (iflag == 0) THEN ! Tsai-Wu
287 DO i=1,nel
288 dmg(i,1) = max(dmg(i,2),dmg(i,3),dmg(i,4))
289 ENDDO
290 ELSE ! Crasurv
291 DO i=1,nel
292 dmg(i,1) = max(dmg(i,2),dmg(i,3),dmg(i,4))
293 IF (abs(dmg(i,5)) >= one) THEN
294 dmg(i,1) = max(abs(dmg(i,5))-one,dmg(i,1))
295 ENDIF
296 IF (abs(dmg(i,6)) >= one) THEN
297 dmg(i,1) = max(abs(dmg(i,6))-one,dmg(i,1))
298 ENDIF
299 IF (dmg(i,7) >= one) THEN
300 dmg(i,1) = max(dmg(i,7)-one,dmg(i,1))
301 ENDIF
302 ENDDO
303 ENDIF
304C----------------------------
305C Failure test ---special treatment inside law25 ---
306C----------------------------
307 DO i=1,nel
308 IF (off(i) == off_old(i) .and. off(i) > zero) THEN
309 IF (off(i) == one) THEN
310 IF (ratio < zero) THEN
311 failnpt = npttot - 1 ! NPTT =MAX(NPTT) for IGTYP=51
312 ELSE
313 failnpt = npttot - nint(npttot*(one-ratio))
314c FAILNPT = CEILING(NPTTOT*RATIO)
315 END IF
316 IF (ioff < 0) ioff=-(ioff+1)
317 joff=0
318C
319 IF (ioff == 0 .AND. wplar(i) >= one) joff=1
320 IF (ioff == 1 .AND. nint(wplar(i)) >= failnpt) joff=1
321 IF (ioff == 2 .AND. nfis1(i) >= failnpt) joff=1
322 IF (ioff == 3 .AND. nfis2(i) >= failnpt) joff=1
323 IF (ioff == 4 .AND. nfis1(i) >= failnpt
324 . .AND. nfis2(i) == failnpt) joff=1
325 IF (ioff == 5 .AND. nfis1(i) >= failnpt) joff=1
326 IF (ioff == 5 .AND. nfis2(i) >= failnpt) joff=1
327 IF (ioff == 6 .AND. nfis3(i) >= failnpt) joff=1
328C
329 IF (joff == 1) THEN
330 off(i) = four_over_5 ! progressive rupture
331 IF (ifailure == 0) ioff_duct(i) = 1 ! flag for progressive rupture
332 ENDIF
333 ELSE IF (ifailure == 0 .and. off(i) < one ) THEN
334 off(i) = off(i)*four_over_5
335 ENDIF
336 ENDIF
337 ENDDO
338c-----------
339 RETURN
340 END
m25crak
subroutine m25crak(nel, off, dmg, l_dmg, dir, ilayer, thly, ngl, strp1, strp2, ply_id, igtyp, ipg, epst1, epst2, epsm1, epsm2, dmax)
Definition m25crak.F:31
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
prony25c
subroutine prony25c(nel, nprony, beta, kv, gv, timestep, rho0, off, dir, epspxx, epspyy, epspxy, epspyz, epspzx, sigvxx, sigvyy, sigvxy, sigvyz, sigvzx, sigv, soundsp, vari, igtyp)
Definition prony25c.F:35
sigeps25c
subroutine sigeps25c(mat_param, nel, pm, off, gstr, dir, thly, time, timestep, shf, ngl, thk0, exx, off_old, eyy, exy, exz, eyz, kxx, kyy, kxy, zz, epsd_pg, rho0, soundsp, uparam, offl, epsd, asrate, sigy, zcfac, nptt, ilay, nfis1, nfis2, nfis3, wplar, npttot, igtyp, sigv, sigply, sigoxx, sigoyy, sigoxy, sigoyz, sigozx, signxx, signyy, signxy, signyz, signzx, sigvxx, sigvyy, sigvxy, sigvyz, sigvzx, israte, uvarv, ishplyxfem, ipt, seq_outp, ply_exx, ply_eyy, ply_exy, ply_exz, ply_eyz, ply_f, pla, crak, ierr, ioff_duct, ifailure, ply_id, ipg, tsaiwu, imconv, iout, dmg, l_dmg)
Definition sigeps25c.F:55