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

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