OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps120_connect_tab_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_tab_dp ../engine/source/materials/mat/mat120/sigeps120_connect_tab_dp.F
25!||--- called by ------------------------------------------------------
26!|| sigeps120_connect_main ../engine/source/materials/mat/mat120/sigeps120_connect_main.F
27!||--- calls -----------------------------------------------------
28!|| table_vinterp ../engine/source/tools/curve/table_tools.F
29!||--- uses -----------------------------------------------------
30!|| interface_table_mod ../engine/share/modules/table_mod.F
31!|| table_mod ../engine/share/modules/table_mod.F
32!||====================================================================
33 SUBROUTINE sigeps120_connect_tab_dp(
34 1 NEL ,NGL ,TIME ,TIMESTEP,UPARAM ,OFF ,
35 2 EPSD ,STIFM ,THICK ,JTHE ,
36 3 AREA ,DEPSZZ ,DEPSYZ ,DEPSZX ,NUPARAM ,
37 4 SIGOZZ ,SIGOYZ ,SIGOZX ,SIGNZZ ,SIGNYZ ,SIGNZX ,
38 5 PLA ,JSMS ,DMELS ,UVAR ,NUVAR ,
39 6 NUMTABL ,ITABLE ,TABLE ,NVARTMP ,VARTMP ,TEMP ,DMG )
40C-----------------------------------------------
41C M o d u l e s
42C-----------------------------------------------
43 USE table_mod
44 USE interface_table_mod
45C-----------------------------------------------
46C I m p l i c i t T y p e s
47C-----------------------------------------------
48#include "implicit_f.inc"
49C-----------------------------------------------
50C C o m m o n B l o c k s
51C-----------------------------------------------
52#include "com04_c.inc"
53#include "units_c.inc"
54#include "comlock.inc"
55#include "sms_c.inc"
56C----------------------------------------------------------------
57C D u m m y A R G U M E N T S
58C----------------------------------------------------------------
59 INTEGER :: NEL,NUPARAM,NUVAR,JSMS,NUMTABL,NVARTMP,JTHE
60 my_real :: TIME,TIMESTEP
61 INTEGER ,DIMENSION(NEL) ,INTENT(IN) :: NGL
62 INTEGER ,DIMENSION(NUMTABL),INTENT(IN) :: ITABLE
63 my_real ,DIMENSION(NUPARAM),INTENT(IN) :: UPARAM
64 my_real ,DIMENSION(NEL) ,INTENT(IN) :: THICK,TEMP,
65 . depszz,depsyz,depszx,
66 . sigozz,sigoyz,sigozx,area
67c
68 my_real ,DIMENSION(NEL) ,INTENT(OUT) :: epsd
69 my_real ,DIMENSION(NEL) ,INTENT(OUT) ::
70 . signzz,signyz,signzx
71 my_real ,DIMENSION(NEL) ,INTENT(INOUT) :: off,pla,stifm,dmels,dmg
72 my_real ,DIMENSION(NEL,NUVAR) ,INTENT(INOUT) :: uvar
73 INTEGER,DIMENSION(NEL,NVARTMP) ,INTENT(INOUT) :: VARTMP
74 TYPE(ttable) ,DIMENSION(NTABLE) ,INTENT(IN) :: TABLE
75C----------------------------------------------------------------
76C L O C A L V A R I B L E S
77C----------------------------------------------------------------
78 INTEGER :: I,II,ITER,NITER,NINDX,NINDF,ITRX,IDAM
79 my_real :: xscale,yscale
80 INTEGER :: IPOS1(NEL,1), IPOS2(NEL,2), IPOS3(NEL,3)
81 INTEGER :: NDIM_YLD,FUNC_YLD
82 my_real :: XVEC1(NEL,1), XVEC2(NEL,2), XVEC3(NEL,3)
83 INTEGER ,DIMENSION(NEL) :: INDX,INDF
84 my_real :: young,nu,g,g2,bulk,yld0,qq,beta,hh,c11,c12,
85 . a1f,a2f,as,d1c,d2c,d1f,d2f,d_trx,d_jc,dm,expn,a1h,a2h
86c Johnson & Cook rate-dependency
87 my_real cc,epsdref,epsdmax,frate,fcut,alpha,alphai,epsdot
88c Local variables
89 my_real facr,triax,rvoce,i1p,ip,jp,tp,fact,facd,dfac,dtinv,
90 . dlam,ddep,dfdlam,dyld_dpla,dpla_dlam,dphi_dlam,dfdpla,dtb,
91 . gamc,gamf,np_xx,np_yy,np_zz,np_xy,np_yz,np_zx,
92 . nf_xx,nf_yy,nf_zz,nf_xy,nf_yz,nf_zx,dpxx,dpyy,dpzz,dpxy,dpyz,dpzx
93 my_real ,DIMENSION(NEL) :: a1, a2, i1,j2,yld,phi,dam,fdam,jcrate,stf,hardp,hardp_i,fvm,
94 . sxx,syy,szz,sxy,syz,szx,dpla,jcr_log,gamma,dgamm,pla_nl,dpdt_nl,yld_i,fdp
95C---------------------------------------------
96
97 young = uparam(1)
98 nu = uparam(2)
99 g = uparam(3)
100 bulk = uparam(4)
101c Parameters of yield function/hardening law and flow potential
102 yld0 = uparam(5)
103 qq = uparam(6)
104 beta = uparam(7)
105 hh = uparam(8)
106 a1f = uparam(9)
107 a2f = uparam(10)
108 a1h = uparam(11)
109 a2h = uparam(12)
110 as = uparam(13)
111c Johnson-Cook failure parameters
112 d1c = uparam(14)
113 d2c = uparam(15)
114 d1f = uparam(16)
115 d2f = uparam(17)
116 d_trx = uparam(18)
117 d_jc = uparam(19)
118 expn = uparam(20)
119c Johnson-Cook strain rate-dependency and distortional hardening
120 cc = uparam(21)
121 epsdref = uparam(22)
122 epsdmax = uparam(23)
123 fcut = uparam(24)
124c
125 itrx = nint(uparam(26))
126 idam = nint(uparam(27))
127c
128 xscale = uparam(30)
129 yscale = uparam(31)
130c
131 func_yld = itable(1)
132 ndim_yld = table(func_yld)%NDIM
133
134c
135 !YOUNG = YOUNG/THICK
136 !G = G/THICK
137 !BULK = BULK/THICK
138 g2 = g * two
139 alpha = 0.005
140 alphai = one - alpha
141 dtinv = one / (max(timestep, em20))
142C---------------------------------------------
143 stf(1:nel) = young * area(1:nel)
144 stifm(1:nel) = stifm(1:nel) + stf(1:nel)*off(1:nel)
145C---------------------------------------------
146 dam(1:nel) = dmg(1:nel)
147 dpla(1:nel) = zero
148 fdam(1:nel) = one - dam(1:nel)
149 DO i = 1,nel
150 IF (off(i) < 0.001) off(i) = zero
151 IF (off(i) < one) off(i) = off(i)*four_over_5
152 IF (off(i) == one) THEN
153 signzz(i) = sigozz(i)/fdam(i) + (depszz(i) /thick(i) )*young
154 signyz(i) = sigoyz(i)/fdam(i) + (depsyz(i) /thick(i) )*g
155 signzx(i) = sigozx(i)/fdam(i) + (depszx(i) /thick(i) )*g
156 i1(i) = signzz(i)
157c
158 sxx(i) = - i1(i) * third
159 syy(i) = - i1(i) * third
160 szz(i) = signzz(i) - i1(i) * third
161 syz(i) = signyz(i)
162 szx(i) = signzx(i)
163 j2(i) = (sxx(i)**2 + syy(i)**2 + szz(i)**2)*half
164 . + syz(i)**2 + szx(i)**2
165 END IF
166 ENDDO
167
168c Johnson & Cook rate dependency (total strain rate)
169 jcrate(1:nel) = one
170 DO i=1,nel
171 IF (off(i) == one) THEN
172 epsdot = sqrt(depsyz(i)**2 + depszx(i)**2 + depszz(i)**2)
173 epsdot = epsdot * dtinv / thick(i)
174 epsd(i) = min(epsdot ,epsdmax)
175 epsd(i) = alpha *epsd(i) + alphai * uvar(i,3)
176 uvar(i,3) = epsd(i)
177 IF (epsd(i) > epsdref) THEN
178 jcr_log(i) = log(epsd(i) / epsdref)
179 jcrate(i) = one + cc * jcr_log(i)
180 END IF
181 END IF
182 ENDDO
183c----------------------------------------------------
184c Computation of the initial yield stress
185c----------------------------------------------------
186 IF (ndim_yld == 1) THEN
187 xvec1(1:nel,1) = pla(1:nel)
188 ipos1(1:nel,1) = vartmp(1:nel,1)
189c
190 CALL table_vinterp(table(func_yld),nel,nel,ipos1,xvec1,yld,hardp)
191c
192 yld(1:nel) = yld(1:nel) * yscale * jcrate(1:nel)
193 hardp(1:nel) = hardp(1:nel) * yscale * jcrate(1:nel)
194 vartmp(1:nel,1) = ipos1(1:nel,1)
195
196 ELSE IF (ndim_yld == 2) THEN
197 xvec2(1:nel,1) = pla(1:nel)
198 xvec2(1:nel,2) = epsd(1:nel)
199 ipos2(1:nel,1) = vartmp(1:nel,1)
200 ipos2(1:nel,2) = 1
201c
202 CALL table_vinterp(table(func_yld),nel,nel,ipos2,xvec2,yld,hardp)
203c
204 yld(1:nel) = yld(1:nel) * yscale
205 hardp(1:nel) = hardp(1:nel) * yscale
206 vartmp(1:nel,1) = ipos2(1:nel,1)
207
208 ELSE
209 xvec3(1:nel,1) = pla(1:nel)
210 xvec3(1:nel,2) = epsd(1:nel)
211 IF (jthe > 0) THEN
212 xvec3(1:nel,3) = temp(1:nel)
213 ELSE
214 xvec3(1:nel,3) = zero
215 END IF
216 ipos3(1:nel,1) = vartmp(1:nel,1)
217 ipos3(1:nel,2) = 1
218 ipos3(1:nel,3) = 1
219c
220 CALL table_vinterp(table(func_yld),nel,nel,ipos3,xvec3,yld,hardp)
221c
222 yld(1:nel) = yld(1:nel) * yscale
223 hardp(1:nel) = hardp(1:nel) * yscale
224 vartmp(1:nel,1) = ipos3(1:nel,1)
225
226 END IF
227c
228 DO i = 1,nel
229 a1(i) = a1f + a1h * pla(i)
230 a2(i) = max(zero, a2f + a2h * pla(i))
231 fdp(i) = third*(sqr3*yld0*a1(i)*i1(i) + a2(i)*max(zero,i1(i))**2)
232 phi(i) = j2(i) + fdp(i) - yld(i)**2
233 ENDDO
234c------------------------------------------------------------------
235c Check plasticity
236c------------------------------------------------------------------
237 nindx = 0
238 nindf = 0
239 DO i = 1,nel
240 IF (phi(i) > zero .and. off(i) == one) THEN
241 nindx = nindx + 1
242 indx(nindx) = i
243 END IF
244 ENDDO
245c
246 IF (nindx > 0) THEN
247
248 niter = 4
249 dpla(1:nel) = zero
250
251 DO iter = 1,niter
252c
253 DO ii = 1,nindx
254 i = indx(ii)
255
256 ! normal to non associated plastic flow surface = d_psi/d_sig
257 i1p = max(zero, i1(i))
258 jp = third * as * i1p
259 tp = two * jp
260c
261 np_xx = sxx(i)
262 np_yy = syy(i)
263 np_zz = szz(i) + tp
264 np_yz = syz(i)
265 np_zx = szx(i)
266c
267 ! normal to plastic yield surface = d_phi/d_sig
268 ip = third*(sqr3*yld0*a1(i) + two*a2(i)*i1p)
269 nf_xx = sxx(i)
270 nf_yy = syy(i)
271 nf_zz = szz(i) + ip
272 nf_yz = syz(i)
273 nf_zx = szx(i)
274c
275 ! derivative of yld criterion : dphi/dlam = d_phi/d_sig * dsig/dlam
276 dfdlam = - nf_zz * (young * np_zz )
277 . -(nf_yz*np_yz + nf_zx*np_zx)*two*g2
278
279 dfdpla = third*(sqr3*a1h*yld0*i1(i) + a2h*i1p**2)
280
281 dpla_dlam = two*sqrt((j2(i) + tp*as*i1(i)))
282c
283 dphi_dlam = dfdlam + (dfdpla - two*yld(i)*hardp(i))*dpla_dlam
284c----------------------------------------------------------------
285 dlam = -phi(i) / dphi_dlam
286c----------------------------------------------------------------
287 ! Plastic strains tensor update
288 dpzz = dlam*np_zz
289 dpyz = dlam*np_yz
290 dpzx = dlam*np_zx
291 ! Elasto-plastic stresses update
292 signzz(i) = signzz(i) - young *dpzz
293 signyz(i) = signyz(i) - g2 *dpyz
294 signzx(i) = signzx(i) - g2 *dpzx
295c
296 ! Plastic strain increments update
297 ddep = dlam * dpla_dlam
298 dpla(i)= max(zero, dpla(i) + ddep)
299 pla(i) = pla(i) + ddep
300 ENDDO ! II = 1,NINDX
301c-----------------------------------------
302c Update yield from tabulated data
303c-----------------------------------------
304c
305 IF (ndim_yld == 1) THEN
306 DO ii = 1, nindx
307 i = indx(ii)
308 xvec1(ii,1) = pla(i)
309 ipos1(ii,1) = vartmp(i,1)
310 ENDDO
311c
312 CALL table_vinterp(table(func_yld),nel,nel,ipos1,xvec1,yld_i,hardp_i)
313c
314 DO ii = 1, nindx
315 i = indx(ii)
316 yld(i) = yld_i(ii)*yscale
317 hardp(i) = hardp_i(ii)*yscale
318 vartmp(i,1) = ipos1(ii,1)
319 ENDDO
320
321 ELSE IF (ndim_yld == 2) THEN
322 DO ii = 1, nindx
323 i = indx(ii)
324 xvec2(ii,1) = pla(i)
325 xvec2(ii,2) = epsd(i)
326 ipos2(ii,1) = vartmp(i,1)
327 ENDDO
328c
329 CALL table_vinterp(table(func_yld),nel,nindx,ipos2,xvec2,yld_i,hardp_i)
330c
331 DO ii = 1, nindx
332 i = indx(ii)
333 yld(i) = yld_i(ii)*yscale
334 hardp(i) = hardp_i(ii)*yscale
335 vartmp(i,1) = ipos2(ii,1)
336 ENDDO
337
338 ELSE
339 DO ii = 1, nindx
340 i = indx(ii)
341 xvec3(ii,1) = pla(i)
342 xvec3(ii,2) = epsd(i)
343 xvec3(ii,3) = temp(i)
344 ipos3(ii,1) = vartmp(i,1)
345 ENDDO
346c
347 CALL table_vinterp(table(func_yld),nel,nel,ipos3,xvec3,yld_i,hardp_i)
348c
349 DO ii = 1, nindx
350 i = indx(ii)
351 yld(i) = yld_i(ii)*yscale
352 hardp(i) = hardp_i(ii)*yscale
353 vartmp(i,1) = ipos3(ii,1)
354 ENDDO
355 END IF
356c----------------------
357c Update stress invariants and criterion function value
358c----------------------
359 DO ii = 1,nindx
360 i = indx(ii)
361 i1(i) = signzz(i)
362 i1p = max(zero, i1(i))
363 sxx(i) = - i1(i)*third
364 syy(i) = - i1(i)*third
365 szz(i) = signzz(i) - i1(i)*third
366 syz(i) = signyz(i)
367 szx(i) = signzx(i)
368 j2(i) =(sxx(i)**2 + syy(i)**2 + szz(i)**2 ) * half
369 . + syz(i)**2 + szx(i)**2
370 a1(i) = a1f + a1h * pla(i)
371 a2(i) = max(zero, a2f + a2h * pla(i))
372 fdp(i) = third*(sqr3*yld0*a1(i)*i1(i) + a2(i)*i1p**2)
373 phi(i) = j2(i) + fdp(i) - yld(i)**2
374 uvar(i,4) = phi(i) / yld(i)**2
375 END DO
376c
377 END DO ! Newton iterations
378c
379 END IF ! NINDX > 0
380c---------------------------------------------------------------------
381c Update damage
382c---------------------------------------------------------------------
383 IF (idam > 0) THEN
384 IF (idam == 1) THEN
385 dgamm(1:nel) = dpla(1:nel)
386 gamma(1:nel) = pla(1:nel)
387 ELSE
388 dgamm(1:nel) = dpla(1:nel) * fdam(1:nel)
389 uvar(1:nel,1) = uvar(1:nel,1) + dgamm(1:nel)
390 gamma(1:nel) = uvar(1:nel,1)
391 END IF
392c
393 DO ii = 1,nindx
394 i = indx(ii)
395c
396 triax = zero
397 fact = one
398 facr = one + d_jc * jcr_log(i) ! total strain rate factor on damage
399 IF ( j2(i)/= zero)triax = third * i1(i) / sqrt(three*j2(i))
400 IF (itrx == 1 ) THEN
401 fact = exp(-d_trx*triax)
402 ELSE
403 IF (triax > zero )fact = exp(-d_trx*triax)
404 ENDIF
405 gamc = (d1c + d2c * fact) * facr
406 gamf = (d1f + d2f * fact) * facr
407 IF (gamma(i) > gamc) THEN
408 IF (dam(i) == zero) THEN
409 nindf = nindf + 1
410 indf(nindf) = i
411 END IF
412 IF (expn == one) THEN
413 dam(i) = dam(i) + dgamm(i) / (gamf - gamc)
414 ELSE
415 dfac = (gamma(i) - gamc) / (gamf - gamc)
416 dam(i) = dam(i) + expn * dfac**(expn-one) * dgamm(i) / (gamf - gamc)
417 END IF
418 IF (dam(i) >= one) THEN
419 dam(i) = one
420 off(i) = four_over_5
421 END IF
422 dmg(i) = dam(i)
423 END IF ! GAMMA > GAMC
424 ENDDO
425c Update damaged stress
426 DO i = 1, nel
427 facd = one - dam(i)
428 signzz(i) = signzz(i) * facd
429 signyz(i) = signyz(i) * facd
430 signzx(i) = signzx(i) * facd
431 END DO
432 END IF
433c-------------------------
434 IF (nindf > 0) THEN
435 DO ii=1,nindf
436#include "lockon.inc"
437 WRITE(iout, 1000) ngl(indf(ii))
438 WRITE(istdo,1100) ngl(indf(ii)),time
439#include "lockoff.inc"
440 ENDDO
441 END IF
442c-----------------------------------------------------------------
443 1000 FORMAT(1x,'START DAMAGE IN SOLID ELEMENT NUMBER ',i10)
444 1100 FORMAT(1x,'START DAMAGE IN SOLID ELEMENT NUMBER ',i10,1x,' AT TIME :',g11.4)
445c-----------------------------------------------------------------
446c-----------------------------------------------------
447c omega = sqrt(2k/2*dmels), dt=2/omega, 2*dmels=dt**2 * 2k / 4
448 IF (idtmins==2 .AND. jsms/=0 ) THEN
449 dtb = (dtmins/dtfacs)**2
450 DO i=1,nel
451 dmels(i)=max(dmels(i),half*dtb*stf(i)*off(i))
452 ENDDO
453 ENDIF
454c-----------------------------------------------------------------
455 RETURN
456 END
alpha
#define alpha
Definition eval.h:35
area
subroutine area(d1, x, x2, y, y2, eint, stif0)
Definition hm_read_prop32.F:567
interface_table_mod::table_vinterp
Definition table_mod.F:251
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
interface_table_mod
Definition table_mod.F:238
table_mod
Definition table_mod.F:112
sigeps120_connect_tab_dp
subroutine sigeps120_connect_tab_dp(nel, ngl, time, timestep, uparam, off, epsd, stifm, thick, jthe, area, depszz, depsyz, depszx, nuparam, sigozz, sigoyz, sigozx, signzz, signyz, signzx, pla, jsms, dmels, uvar, nuvar, numtabl, itable, table, nvartmp, vartmp, temp, dmg)
Definition sigeps120_connect_tab_dp.F:40
table_mod::ttable
Definition table_mod.F:125