OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps120_connect_tab_vm.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_vm (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_vm()

subroutine sigeps120_connect_tab_vm ( 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_vm.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
84c material parameters
85 my_real :: young,nu,g,g2,bulk,yld0,qq,beta,hh,c11,c12,
86 . a1f,a2f,as,d1c,d2c,d1f,d2f,d_trx,d_jc,dm,expn
87c Johnson & Cook rate-dependency
88 my_real cc,epsdref,epsdmax,frate,fcut,alpha,alphai,epsdot
89c Local variables
90 my_real facr,triax,rvoce,i1p,ip,jp,tp,fact,facd,dfac,dtinv,
91 . dlam,ddep,dfdlam,dyld_dpla,dpla_dlam,dphi_dlam,dtb,
92 . gamc,gamf,np_xx,np_yy,np_zz,np_xy,np_yz,np_zx,
93 . nf_xx,nf_yy,nf_zz,nf_xy,nf_yz,nf_zx,dpxx,dpyy,dpzz,dpxy,dpyz,dpzx
94 my_real ,DIMENSION(NEL) :: i1,j2,yld,phi,dam,fdam,jcrate,stf,hardp,hardp_i,fvm,
95 . sxx,syy,szz,sxy,syz,szx,dpla,jcr_log,gamma,dgamm,pla_nl,dpdt_nl,yld_i
96C---------------------------------------------
97
98 young = uparam(1)
99 nu = uparam(2)
100 g = uparam(3)
101 bulk = uparam(4)
102c Parameters of yield function/hardening law and flow potential
103 yld0 = uparam(5)
104 qq = uparam(6)
105 beta = uparam(7)
106 hh = uparam(8)
107 a1f = uparam(9)
108 a2f = uparam(10)
109c A1H = UPARAM(11)
110c A2H = UPARAM(12)
111 as = uparam(13)
112c Johnson-Cook failure parameters
113 d1c = uparam(14)
114 d2c = uparam(15)
115 d1f = uparam(16)
116 d2f = uparam(17)
117 d_trx = uparam(18)
118 d_jc = uparam(19)
119 expn = uparam(20)
120c Johnson-Cook strain rate-dependency and distortional hardening
121 cc = uparam(21)
122 epsdref = uparam(22)
123 epsdmax = uparam(23)
124 fcut = uparam(24)
125c
126 itrx = nint(uparam(26))
127 idam = nint(uparam(27))
128c
129 xscale = uparam(30)
130 yscale = uparam(31)
131c
132 func_yld = itable(1)
133 ndim_yld = table(func_yld)%NDIM
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 fvm(i) = max(zero, i1(i) + half*sqr3*a1f*yld0 / a2f )
230 phi(i) = j2(i) + third*a2f * fvm(i)**2
231 . - (yld(i)**2 + fourth*(a1f*yld0)**2/a2f)
232 ENDDO
233c------------------------------------------------------------------
234c Check plasticity
235c------------------------------------------------------------------
236 nindx = 0
237 nindf = 0
238 DO i = 1,nel
239 IF (phi(i) > zero .and. off(i) == one) THEN
240 nindx = nindx + 1
241 indx(nindx) = i
242 END IF
243 ENDDO
244c
245 IF (nindx > 0) THEN
246
247 niter = 4
248 dpla(1:nel) = zero
249
250 DO iter = 1,niter
251c
252 DO ii = 1,nindx
253 i = indx(ii)
254
255 ! normal to non associated plastic flow surface = d_psi/d_sig
256 i1p = max(zero, i1(i))
257 jp = third * as * i1p
258 tp = two * jp
259c
260 np_xx = sxx(i)
261 np_yy = syy(i)
262 np_zz = szz(i) + tp
263 np_yz = syz(i)
264 np_zx = szx(i)
265c
266 ! normal to plastic yield surface = d_phi/d_sig
267 ip = two_third * a2f * fvm(i)
268 nf_xx = sxx(i)
269 nf_yy = syy(i)
270 nf_zz = szz(i) + ip
271 nf_yz = syz(i)
272 nf_zx = szx(i)
273c
274 ! derivative of yld criterion : dphi/dlam = d_phi/d_sig * dsig/dlam
275 dfdlam = - nf_zz * (young * np_zz)
276 . -( nf_yz*np_yz + nf_zx*np_zx)*two*g2
277c
278 dpla_dlam = two*sqrt((j2(i) + tp*as*i1(i)))
279 dphi_dlam = dfdlam - two*yld(i)*hardp(i)*dpla_dlam
280c----------------------------------------------------------------
281 dlam = -phi(i) / dphi_dlam
282c----------------------------------------------------------------
283 ! Plastic strains tensor update
284
285 dpzz = dlam*np_zz
286 dpyz = dlam*np_yz
287 dpzx = dlam*np_zx
288 ! Elasto-plastic stresses update
289
290 signzz(i) = signzz(i) - young*dpzz
291 signyz(i) = signyz(i) - g2*dpyz
292 signzx(i) = signzx(i) - g2*dpzx
293c
294 ! Plastic strain increments update
295 ddep = dlam * dpla_dlam
296 dpla(i)= max(zero, dpla(i) + ddep)
297 pla(i) = pla(i) + ddep
298 ENDDO ! II = 1,NINDX
299c-----------------------------------------
300c Update yield from tabulated data
301c-----------------------------------------
302c
303 IF (ndim_yld == 1) THEN
304 DO ii = 1, nindx
305 i = indx(ii)
306 xvec1(ii,1) = pla(i)
307 ipos1(ii,1) = vartmp(i,1)
308 ENDDO
309c
310 CALL table_vinterp(table(func_yld),nel,nel,ipos1,xvec1,yld_i,hardp_i)
311c
312 DO ii = 1, nindx
313 i = indx(ii)
314 yld(i) = yld_i(ii)*yscale
315 hardp(i) = hardp_i(ii)*yscale
316 vartmp(i,1) = ipos1(ii,1)
317 ENDDO
318
319 ELSE IF (ndim_yld == 2) THEN
320 DO ii = 1, nindx
321 i = indx(ii)
322 xvec2(ii,1) = pla(i)
323 xvec2(ii,2) = epsd(i)
324 ipos2(ii,1) = vartmp(i,1)
325 ENDDO
326c
327 CALL table_vinterp(table(func_yld),nel,nindx,ipos2,xvec2,yld_i,hardp_i)
328c
329 DO ii = 1, nindx
330 i = indx(ii)
331 yld(i) = yld_i(ii)*yscale
332 hardp(i) = hardp_i(ii)*yscale
333 vartmp(i,1) = ipos2(ii,1)
334 ENDDO
335
336 ELSE
337 DO ii = 1, nindx
338 i = indx(ii)
339 xvec3(ii,1) = pla(i)
340 xvec3(ii,2) = epsd(i)
341 xvec3(ii,3) = temp(i)
342 ipos3(ii,1) = vartmp(i,1)
343 ENDDO
344c
345 CALL table_vinterp(table(func_yld),nel,nel,ipos3,xvec3,yld_i,hardp_i)
346c
347 DO ii = 1, nindx
348 i = indx(ii)
349 yld(i) = yld_i(ii)*yscale
350 hardp(i) = hardp_i(ii)*yscale
351 vartmp(i,1) = ipos3(ii,1)
352 ENDDO
353 END IF
354c----------------------
355c Update stress invariants and criterion function value
356c----------------------
357 DO ii = 1,nindx
358 i = indx(ii)
359 i1(i) = signzz(i)
360 sxx(i) = - i1(i)*third
361 syy(i) = - i1(i)*third
362 szz(i) = signzz(i) - i1(i)*third
363 syz(i) = signyz(i)
364 szx(i) = signzx(i)
365 j2(i) =(sxx(i)**2 + syy(i)**2 + szz(i)**2 ) * half
366 . + syz(i)**2 + szx(i)**2
367 fvm(i) = max(zero, i1(i) + half*sqr3*a1f*yld0 / a2f )
368 phi(i) = j2(i) + third*a2f * fvm(i)**2
369 . - (yld(i)**2 + fourth*(a1f*yld0)**2/a2f)
370 uvar(i,4) = phi(i) / yld(i)**2
371 END DO
372c
373 END DO ! Newton iterations
374c
375 END IF ! NINDX > 0
376c---------------------------------------------------------------------
377c Update damage
378c---------------------------------------------------------------------
379 IF (idam > 0) THEN
380 IF (idam == 1) THEN
381 dgamm(1:nel) = dpla(1:nel)
382 gamma(1:nel) = pla(1:nel)
383 ELSE
384 dgamm(1:nel) = dpla(1:nel) * fdam(1:nel)
385 uvar(1:nel,1) = uvar(1:nel,1) + dgamm(1:nel)
386 gamma(1:nel) = uvar(1:nel,1)
387 END IF
388c
389 DO ii = 1,nindx
390 i = indx(ii)
391c
392 triax = zero
393 fact = one
394 facr = one + d_jc * jcr_log(i) ! total strain rate factor on damage
395 IF ( j2(i)/= zero)triax = third * i1(i) / sqrt(three*j2(i))
396 IF (itrx == 1 ) THEN
397 fact = exp(-d_trx*triax)
398 ELSE
399 IF (triax > zero )fact = exp(-d_trx*triax)
400 ENDIF
401 gamc = (d1c + d2c * fact) * facr
402 gamf = (d1f + d2f * fact) * facr
403 IF (gamma(i) > gamc) THEN
404 IF (dam(i) == zero) THEN
405 nindf = nindf + 1
406 indf(nindf) = i
407 END IF
408 IF (expn == one) THEN
409 dam(i) = dam(i) + dgamm(i) / (gamf - gamc)
410 ELSE
411 dfac = (gamma(i) - gamc) / (gamf - gamc)
412 dam(i) = dam(i) + expn * dfac**(expn-one) * dgamm(i) / (gamf - gamc)
413 END IF
414 IF (dam(i) >= one) THEN
415 dam(i) = one
416 off(i) = four_over_5
417 END IF
418 dmg(i) = dam(i)
419 END IF ! GAMMA > GAMC
420 ENDDO
421c Update damaged stress
422 DO i = 1, nel
423 facd = one - dam(i)
424 signzz(i) = signzz(i) * facd
425 signyz(i) = signyz(i) * facd
426 signzx(i) = signzx(i) * facd
427 END DO
428 END IF
429c-------------------------
430 IF (nindf > 0) THEN
431 DO ii=1,nindf
432#include "lockon.inc"
433 WRITE(iout, 1000) ngl(indf(ii))
434 WRITE(istdo,1100) ngl(indf(ii)),time
435#include "lockoff.inc"
436 ENDDO
437 END IF
438c-----------------------------------------------------------------
439 1000 FORMAT(1x,'START DAMAGE IN SOLID ELEMENT NUMBER ',i10)
440 1100 FORMAT(1x,'START DAMAGE IN SOLID ELEMENT NUMBER ',i10,1x,' AT TIME :',g11.4)
441c-----------------------------------------------------------------
442c-----------------------------------------------------
443c omega = sqrt(2k/2*dmels), dt=2/omega, 2*dmels=dt**2 * 2k / 4
444 IF (idtmins==2 .AND. jsms/=0 ) THEN
445 dtb = (dtmins/dtfacs)**2
446 DO i=1,nel
447 dmels(i)=max(dmels(i),half*dtb*stf(i)*off(i))
448 ENDDO
449 ENDIF
450c-----------------------------------------------------------------
451 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