OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps120_connect_tab_dp.F File Reference
#include "implicit_f.inc"
#include "com04_c.inc"
#include "units_c.inc"
#include "comlock.inc"
#include "sms_c.inc"
#include "lockon.inc"
#include "lockoff.inc"

Go to the source code of this file.

Functions/Subroutines

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)

Function/Subroutine Documentation

◆ sigeps120_connect_tab_dp()

subroutine sigeps120_connect_tab_dp ( integer nel,
integer, dimension(nel), intent(in) ngl,
time,
timestep,
intent(in) uparam,
intent(inout) off,
intent(out) epsd,
intent(inout) stifm,
intent(in) thick,
integer jthe,
intent(in) area,
intent(in) depszz,
intent(in) depsyz,
intent(in) depszx,
integer nuparam,
intent(in) sigozz,
intent(in) sigoyz,
intent(in) sigozx,
intent(out) signzz,
intent(out) signyz,
intent(out) signzx,
intent(inout) pla,
integer jsms,
intent(inout) dmels,
intent(inout) uvar,
integer nuvar,
integer numtabl,
integer, dimension(numtabl), intent(in) itable,
type(ttable), dimension(ntable), intent(in) table,
integer nvartmp,
integer, dimension(nel,nvartmp), intent(inout) vartmp,
intent(in) temp,
intent(inout) dmg )

Definition at line 33 of file sigeps120_connect_tab_dp.F.

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
my_real
#define my_real
Definition cppsort.cpp:32
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