OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps120_connect_dp.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_dp ../engine/source/materials/mat/mat120/sigeps120_connect_dp.F
25!||--- called by ------------------------------------------------------
26!|| sigeps120_connect_main ../engine/source/materials/mat/mat120/sigeps120_connect_main.F
27!||====================================================================
28 SUBROUTINE sigeps120_connect_dp(
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,A1H,A2H
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,fdp,i1p,ip,jp,tp,fact,facd,dfac,dtinv,
72 . dlam,ddep,dfdlam,dyld_dpla,dpla_dlam,dphi_dlam,dfdpla,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,a1,a2,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)
90 a1h = uparam(11)
91 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
143c Johnson & Cook rate dependency (total strain rate)
144 jcrate(1:nel) = one
145 DO i=1,nel
146 IF (off(i) == one) THEN
147 epsdot = sqrt(depsyz(i)**2 + depszx(i)**2 + depszz(i)**2)
148 epsdot = epsdot * dtinv / thick(i)
149 epsd(i) = min(epsdot ,epsdmax)
150 epsd(i) = alpha *epsd(i) + alphai * uvar(i,3)
151 uvar(i,3) = epsd(i)
152 IF (epsd(i) > epsdref) THEN
153 jcr_log(i) = log(epsd(i) / epsdref)
154 jcrate(i) = one + cc * jcr_log(i)
155 END IF
156 END IF
157 ENDDO
158c------------------------------------------------------------------
159c Current yield value, yield criterion and plastic flow fONEction
160c------------------------------------------------------------------
161 DO i = 1,nel
162 rvoce = qq*(one - exp(-beta*pla(i))) + hh*pla(i)
163 yld(i)= (yld0 + rvoce) * jcrate(i)
164 a1(i) = a1f + a1h * pla(i)
165 a2(i) = max(zero, a2f + a2h * pla(i))
166 fdp = third*(sqr3*yld0*a1(i)*i1(i) + a2(i)*max(zero,i1(i) )**2 )
167 phi(i)= j2(i) + fdp - yld(i)**2
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
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 = third*(sqr3*yld0*a1(i) + two*a2(i)*i1p)
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 dfdpla = third*(sqr3*a1h*yld0*i1(i) + a2h*i1p**2)
214
215 dpla_dlam = two*sqrt((j2(i) + tp*as*i1(i)))
216c
217 dphi_dlam = dfdlam + (dfdpla - two*yld(i)*dyld_dpla)*dpla_dlam
218c----------------------------------------------------------------
219 dlam = -phi(i) / dphi_dlam
220c----------------------------------------------------------------
221 ! Plastic strains tensor update
222 dpzz = dlam*np_zz
223 dpyz = dlam*np_yz
224 dpzx = dlam*np_zx
225 ! Elasto-plastic stresses update
226 signzz(i) = signzz(i) - young *dpzz
227 signyz(i) = signyz(i) - g2 *dpyz
228 signzx(i) = signzx(i) - g2 *dpzx
229c
230 ! Plastic strain increments update
231 ddep = dlam * dpla_dlam
232 dpla(i)= max(zero, dpla(i) + ddep)
233 pla(i) = pla(i) + ddep
234c----------------------
235c criterion update
236c----------------------
237 i1(i) = signzz(i)
238 i1p = max(zero, i1(i))
239 sxx(i) = - i1(i)*third
240 syy(i) = - i1(i)*third
241 szz(i) = signzz(i) - i1(i)*third
242 syz(i) = signyz(i)
243 szx(i) = signzx(i)
244 j2(i) =(sxx(i)**2 + syy(i)**2 + szz(i)**2 ) * half
245 . + syz(i)**2 + szx(i)**2
246 rvoce = qq*(one - exp(-beta*pla(i))) + hh*pla(i)
247 yld(i) = (yld0 + rvoce) * jcrate(i)
248 a1(i) = a1f + a1h * pla(i)
249 a2(i) = max(zero, a2f + a2h * pla(i))
250 fdp = third*(sqr3*yld0*a1(i)*i1(i) + a2(i)*i1p**2)
251 phi(i) = j2(i) + fdp - yld(i)**2
252 END DO ! Newton iterations
253c
254 !UVAR(I,4) = PHI(I) / YLD0**2
255 ENDDO ! II = 1,NINDX
256
257c---------------------------------------------------------------------
258c Update damage
259c---------------------------------------------------------------------
260 IF (idam > 0) THEN
261 IF (idam == 1) THEN
262 dgamm(1:nel) = dpla(1:nel)
263 gamma(1:nel) = pla(1:nel)
264 ELSE
265 dgamm(1:nel) = dpla(1:nel) * fdam(1:nel)
266 uvar(1:nel,1) = uvar(1:nel,1) + dgamm(1:nel)
267 gamma(1:nel) = uvar(1:nel,1)
268 END IF
269c
270 DO ii = 1,nindx
271 i = indx(ii)
272 triax = zero
273 fact = one
274 facr = one + d_jc * jcr_log(i) ! total strain rate factor on damage
275 IF ( j2(i)/= zero)triax = third * i1(i) / sqrt(three*j2(i))
276 IF (itrx == 1 ) THEN
277 fact = exp(-d_trx*triax)
278 ELSE
279 IF (triax > zero )fact = exp(-d_trx*triax)
280 ENDIF
281 gamc = (d1c + d2c * fact) * facr
282 gamf = (d1f + d2f * fact) * facr
283 IF (gamma(i) > gamc) THEN
284 IF (dam(i) == zero) THEN
285 nindf = nindf + 1
286 indf(nindf) = i
287 END IF
288 IF (expn == one) THEN
289 dam(i) = dam(i) + dgamm(i) / (gamf - gamc)
290 ELSE
291 dfac = (gamma(i) - gamc) / (gamf - gamc)
292 dam(i) = dam(i) + expn * dfac**(expn-one) * dgamm(i) / (gamf - gamc)
293 END IF
294 IF (dam(i) >= one) THEN
295 dam(i) = one
296 off(i) = four/five
297 END IF
298 dmg(i) = dam(i)
299 END IF ! GAMMA > GAMC
300 ENDDO
301c Update damaged stress
302 DO i = 1, nel
303 facd = one - dam(i)
304 signzz(i) = signzz(i) * facd
305 signyz(i) = signyz(i) * facd
306 signzx(i) = signzx(i) * facd
307 END DO
308 END IF
309c-------------------------
310 IF (nindf > 0) THEN
311 DO ii=1,nindf
312#include "lockon.inc"
313 WRITE(iout, 1000) ngl(indf(ii))
314 WRITE(istdo,1100) ngl(indf(ii)),time
315#include "lockoff.inc"
316 ENDDO
317 END IF
318c-----------------------------------------------------------------
319 1000 FORMAT(1x,'START DAMAGE IN SOLID ELEMENT NUMBER ',i10)
320 1100 FORMAT(1x,'START DAMAGE IN SOLID ELEMENT NUMBER ',i10,1x,' AT TIME :',g11.4)
321c-----------------------------------------------------------------
322c-----------------------------------------------------
323c omega = sqrt(2k/2*dmels), dt=2/omega, 2*dmels=dt**2 * 2k / 4
324 IF (idtmins==2 .AND. jsms/=0 ) THEN
325 dtb = (dtmins/dtfacs)**2
326 DO i=1,nel
327 dmels(i)=max(dmels(i),half*dtb*stf(i)*off(i))
328 ENDDO
329 ENDIF
330c-----------------------------------------------------------------
331 RETURN
332 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_dp
subroutine sigeps120_connect_dp(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_dp.F:34