OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps120_connect_dp.F File Reference
#include "implicit_f.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_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)

Function/Subroutine Documentation

◆ sigeps120_connect_dp()

subroutine sigeps120_connect_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,
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,
intent(inout) dmg )

Definition at line 28 of file sigeps120_connect_dp.F.

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
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
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21