OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps120_connect_vm.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!|| sigeps120_connect_vm ../engine/source/materials/mat/mat120/sigeps120_connect_vm.F
25!||--- called by ------------------------------------------------------
26!|| sigeps120_connect_main ../engine/source/materials/mat/mat120/sigeps120_connect_main.F
27!||====================================================================
28 SUBROUTINE sigeps120_connect_vm(
29 1 NEL ,NGL ,TIME ,TIMESTEP,UPARAM ,OFF ,
30 2 EPSD ,STIFM ,THICK ,
31 3 AREA ,DEPSZZ ,DEPSYZ ,DEPSZX ,NUPARAM ,
32 4 SIGOZZ ,SIGOYZ ,SIGOZX ,SIGNZZ ,SIGNYZ ,SIGNZX ,
33 5 PLA ,JSMS ,DMELS , UVAR ,NUVAR ,DMG)
34C-----------------------------------------------
35C I m p l i c i t T y p e s
36C-----------------------------------------------
37#include "implicit_f.inc"
38C-----------------------------------------------
39C C o m m o n B l o c k s
40C-----------------------------------------------
41#include "units_c.inc"
42#include "comlock.inc"
43#include "sms_c.inc"
44C----------------------------------------------------------------
45C D u m m y A R G U M E N T S
46C----------------------------------------------------------------
47 INTEGER :: NEL,NUPARAM,NUVAR,JSMS
48 my_real :: TIME,TIMESTEP
49 INTEGER ,DIMENSION(NEL) ,INTENT(IN) :: NGL
50 my_real ,DIMENSION(NUPARAM),INTENT(IN) :: UPARAM
51 my_real ,DIMENSION(NEL) ,INTENT(IN) :: THICK,
52 . depszz,depsyz,depszx,
53 . sigozz,sigoyz,sigozx,area
54c
55 my_real ,DIMENSION(NEL) ,INTENT(OUT) :: epsd
56 my_real ,DIMENSION(NEL) ,INTENT(OUT) ::
57 . signzz,signyz,signzx
58 my_real ,DIMENSION(NEL) ,INTENT(INOUT) :: off,pla,stifm,dmels,dmg
59 my_real ,DIMENSION(NEL,NUVAR) ,INTENT(INOUT) :: uvar
60C----------------------------------------------------------------
61C L O C A L V A R I B L E S
62C----------------------------------------------------------------
63 INTEGER :: I,II,ITER,NITER,NINDX,NINDF,ITRX,IDAM
64 INTEGER ,DIMENSION(NEL) :: INDX,INDF
65c material parameters
66 my_real :: YOUNG,NU,G,G2,BULK,YLD0,QQ,BETA,HH,C11,C12,
67 . A1F,A2F,AS,D1C,D2C,D1F,D2F,D_TRX,D_JC,DM,EXPN
68c Johnson & Cook rate-dependency
69 my_real cc,epsdref,epsdmax,frate,fcut,alpha,alphai,epsdot
70c Local variables
71 my_real facr,triax,rvoce,fvm,i1p,ip,jp,tp,fact,facd,dfac,dtinv,
72 . dlam,ddep,dfdlam,dyld_dpla,dpla_dlam,dphi_dlam,dtb,
73 . gamc,gamf,np_xx,np_yy,np_zz,np_xy,np_yz,np_zx,
74 . nf_xx,nf_yy,nf_zz,nf_xy,nf_yz,nf_zx,dpxx,dpyy,dpzz,dpxy,dpyz,dpzx
75 my_real ,DIMENSION(NEL) :: i1,j2,yld,phi,dam,fdam,jcrate,stf,
76 . sxx,syy,szz,sxy,syz,szx,dpla,jcr_log,gamma,dgamm,pla_nl,dpdt_nl
77C---------------------------------------------
78
79 young = uparam(1)
80 nu = uparam(2)
81 g = uparam(3)
82 bulk = uparam(4)
83c Parameters of yield function/hardening law and flow potential
84 yld0 = uparam(5)
85 qq = uparam(6)
86 beta = uparam(7)
87 hh = uparam(8)
88 a1f = uparam(9)
89 a2f = uparam(10)
90c A1H = UPARAM(11)
91c A2H = UPARAM(12)
92 as = uparam(13)
93c Johnson-Cook failure parameters
94 d1c = uparam(14)
95 d2c = uparam(15)
96 d1f = uparam(16)
97 d2f = uparam(17)
98 d_trx = uparam(18)
99 d_jc = uparam(19)
100 expn = uparam(20)
101c Johnson-Cook strain rate-dependency and distortional hardening
102 cc = uparam(21)
103 epsdref = uparam(22)
104 epsdmax = uparam(23)
105 fcut = uparam(24)
106c
107 itrx = nint(uparam(26))
108 idam = nint(uparam(27))
109
110c
111 !YOUNG = YOUNG/THICK
112 !G = G/THICK
113 !BULK = BULK/THICK
114 g2 = g * two
115 alpha = 0.005
116 alphai = one - alpha
117 dtinv = one / (max(timestep, em20))
118C---------------------------------------------
119 stf(1:nel) = young * area(1:nel)
120 stifm(1:nel) = stifm(1:nel) + stf(1:nel)*off(1:nel)
121C---------------------------------------------
122 dam(1:nel) = dmg(1:nel)
123 dpla(1:nel) = zero
124 fdam(1:nel) = one - dam(1:nel)
125 DO i = 1,nel
126 IF (off(i) < 0.001) off(i) = zero
127 IF (off(i) < one) off(i) = off(i)*four/five
128 IF (off(i) == one) THEN
129 signzz(i) = sigozz(i)/fdam(i) + (depszz(i) /thick(i) )*young
130 signyz(i) = sigoyz(i)/fdam(i) + (depsyz(i) /thick(i) )*g
131 signzx(i) = sigozx(i)/fdam(i) + (depszx(i) /thick(i) )*g
132 i1(i) = signzz(i)
133c
134 sxx(i) = - i1(i) * third
135 syy(i) = - i1(i) * third
136 szz(i) = signzz(i) - i1(i) * third
137 syz(i) = signyz(i)
138 szx(i) = signzx(i)
139 j2(i) = (sxx(i)**2 + syy(i)**2 + szz(i)**2)*half
140 . + syz(i)**2 + szx(i)**2
141 END IF
142 ENDDO
143
144c Johnson & Cook rate dependency (total strain rate)
145 jcrate(1:nel) = one
146 DO i=1,nel
147 IF (off(i) == one) THEN
148 epsdot = sqrt(depsyz(i)**2 + depszx(i)**2 + depszz(i)**2)
149 epsdot = epsdot * dtinv / thick(i)
150 epsd(i) = min(epsdot ,epsdmax)
151 epsd(i) = alpha *epsd(i) + alphai * uvar(i,3)
152 uvar(i,3) = epsd(i)
153 IF (epsd(i) > epsdref) THEN
154 jcr_log(i) = log(epsd(i) / epsdref)
155 jcrate(i) = one + cc * jcr_log(i)
156 END IF
157 END IF
158 ENDDO
159c------------------------------------------------------------------
160c Current yield value, yield criterion and plastic flow function
161c------------------------------------------------------------------
162 DO i = 1,nel
163 rvoce = qq*(one - exp(-beta*pla(i))) + hh*pla(i)
164 yld(i) = (yld0 + rvoce) * jcrate(i)
165 fvm = max(zero, i1(i) + half*sqr3*a1f*yld0 / a2f )
166 phi(i) = j2(i) + third*a2f * fvm**2
167 . - (yld(i)**2 + fourth*(a1f*yld0)**2/a2f)
168 ENDDO
169
170c------------------------------------------------------------------
171c Plasticity
172c------------------------------------------------------------------
173 nindx = 0
174 nindf = 0
175 DO i = 1,nel
176 IF (phi(i) > zero .and. off(i) == one) THEN
177 nindx = nindx + 1
178 indx(nindx) = i
179 END IF
180 ENDDO
181c
182c
183 DO ii = 1,nindx
184 i = indx(ii)
185 dpla(i) = zero
186 i1p = max(zero, i1(i))
187 fvm = max(zero, i1(i) + half*sqr3*a1f*yld0 / a2f)
188 niter = 4
189 DO iter = 1,niter
190 ! normal to non associated plastic flow surface = d_psi/d_sig
191 jp = third * as * i1p
192 tp = two * third * as * i1p
193c
194 np_xx = sxx(i)
195 np_yy = syy(i)
196 np_zz = szz(i) + tp
197 np_yz = syz(i)
198 np_zx = szx(i)
199c
200 ! normal to plastic yield surface = d_sig_eff/d_sig
201 ip = two_third * a2f * fvm
202 nf_xx = sxx(i)
203 nf_yy = syy(i)
204 nf_zz = szz(i) + ip
205 nf_yz = syz(i)
206 nf_zx = szx(i)
207c
208 ! derivative of yld criterion : dphi/dlam = d_phi/d_sig * dsig/dlam
209 dfdlam = - nf_zz * (young * np_zz )
210 . -(nf_yz*np_yz + nf_zx*np_zx)*two*g2
211c
212 dyld_dpla = (hh + beta*qq*exp(-beta*pla(i))) * jcrate(i)
213 dpla_dlam = two*sqrt((j2(i) + tp*as*i1(i)))
214
215 dphi_dlam = dfdlam - two*yld(i)*dyld_dpla*dpla_dlam
216c----------------------------------------------------------------
217 dlam = -phi(i) / dphi_dlam
218c----------------------------------------------------------------
219 ! Plastic strains tensor update
220 dpzz = dlam*np_zz
221 dpyz = dlam*np_yz
222 dpzx = dlam*np_zx
223 ! Elasto-plastic stresses update
224 signzz(i) = signzz(i) - young *dpzz
225 signyz(i) = signyz(i) - g2 *dpyz
226 signzx(i) = signzx(i) - g2 *dpzx
227c
228 ! Plastic strain increments update
229 ddep = dlam * dpla_dlam
230 dpla(i)= max(zero, dpla(i) + ddep)
231 pla(i) = pla(i) + ddep
232c----------------------
233c criterion update
234c----------------------
235 i1(i) = signzz(i)
236 i1p = max(zero, i1(i))
237 sxx(i) = - i1(i)*third
238 syy(i) = - i1(i)*third
239 szz(i) = signzz(i) - i1(i)*third
240 syz(i) = signyz(i)
241 szx(i) = signzx(i)
242 j2(i) =(sxx(i)**2 + syy(i)**2 + szz(i)**2 ) * half
243 . + syz(i)**2 + szx(i)**2
244 rvoce = qq*(one - exp(-beta*pla(i))) + hh*pla(i)
245 yld(i) = (yld0 + rvoce) * jcrate(i)
246 fvm = max(zero, i1(i) + half*sqr3*a1f*yld0 / a2f)
247 phi(i) = j2(i) + third*a2f * fvm**2
248 . - (yld(i)**2 + fourth*(a1f*yld0)**2/a2f)
249 END DO ! Newton iterations
250c
251 uvar(i,4) = yld(i)**2 + fourth*(a1f*yld0)**2/a2f
252 uvar(i,5) = j2(i) + third*a2f * fvm**2
253 ENDDO ! II = 1,NINDX
254c---------------------------------------------------------------------
255c Update damage
256c---------------------------------------------------------------------
257 IF (idam > 0) THEN
258 IF (idam == 1) THEN
259 dgamm(1:nel) = dpla(1:nel)
260 gamma(1:nel) = pla(1:nel)
261 ELSE
262 dgamm(1:nel) = dpla(1:nel) * fdam(1:nel)
263 uvar(1:nel,1) = uvar(1:nel,1) + dgamm(1:nel)
264 gamma(1:nel) = uvar(1:nel,1)
265 END IF
266c
267 DO ii = 1,nindx
268 i = indx(ii)
269 triax = zero
270 fact = one
271 facr = one + d_jc * jcr_log(i) ! total strain rate factor on damage
272 IF ( j2(i)/= zero)triax = third * i1(i) / sqrt(three*j2(i))
273 IF (itrx == 1 ) THEN
274 fact = exp(-d_trx*triax)
275 ELSE
276 IF (triax > zero )fact = exp(-d_trx*triax)
277 ENDIF
278 gamc = (d1c + d2c * fact) * facr
279 gamf = (d1f + d2f * fact) * facr
280 IF (gamma(i) > gamc) THEN
281 IF (dam(i) == zero) THEN
282 nindf = nindf + 1
283 indf(nindf) = i
284 END IF
285 IF (expn == one) THEN
286 dam(i) = dam(i) + dgamm(i) / (gamf - gamc)
287 ELSE
288 dfac = (gamma(i) - gamc) / (gamf - gamc)
289 dam(i) = dam(i) + expn * dfac**(expn-one) * dgamm(i) / (gamf - gamc)
290 END IF
291 IF (dam(i) >= one) THEN
292 dam(i) = one
293 off(i) = four/five
294 END IF
295 dmg(i) = dam(i)
296 END IF ! GAMMA > GAMC
297 ENDDO
298c Update damaged stress
299 DO i = 1, nel
300 facd = one - dam(i)
301 signzz(i) = signzz(i) * facd
302 signyz(i) = signyz(i) * facd
303 signzx(i) = signzx(i) * facd
304 END DO
305 END IF
306c-------------------------
307 IF (nindf > 0) THEN
308 DO ii=1,nindf
309#include "lockon.inc"
310 WRITE(iout, 1000) ngl(indf(ii))
311 WRITE(istdo,1100) ngl(indf(ii)),time
312#include "lockoff.inc"
313 ENDDO
314 END IF
315c-----------------------------------------------------------------
316 1000 FORMAT(1x,'START DAMAGE IN SOLID ELEMENT NUMBER ',i10)
317 1100 FORMAT(1x,'START DAMAGE IN SOLID ELEMENT NUMBER ',i10,1x,' AT TIME :',g11.4)
318c-----------------------------------------------------------------
319c-----------------------------------------------------
320c omega = sqrt(2k/2*dmels), dt=2/omega, 2*dmels=dt**2 * 2k / 4
321 IF (idtmins==2 .AND. jsms/=0 ) THEN
322 dtb = (dtmins/dtfacs)**2
323 DO i=1,nel
324 dmels(i)=max(dmels(i),half*dtb*stf(i)*off(i))
325 ENDDO
326 ENDIF
327c-----------------------------------------------------------------
328 RETURN
329 END
alpha
#define alpha
Definition eval.h:35
area
subroutine area(d1, x, x2, y, y2, eint, stif0)
Definition hm_read_prop32.F:567
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
sigeps120_connect_vm
subroutine sigeps120_connect_vm(nel, ngl, time, timestep, uparam, off, epsd, stifm, thick, area, depszz, depsyz, depszx, nuparam, sigozz, sigoyz, sigozx, signzz, signyz, signzx, pla, jsms, dmels, uvar, nuvar, dmg)
Definition sigeps120_connect_vm.F:34