OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
cnloc_mat104_ini.F File Reference
#include "implicit_f.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine cnloc_mat104_ini (nel, ipg, ipt, nuparam, nuvar, uparam, uvar, pla, off, thkly, offg, sigoxx, sigoyy, sigoxy, sigoyz, sigozx, thk, dmg, nptr, npts, nptt, bufly, time, varnl, failure)

Function/Subroutine Documentation

◆ cnloc_mat104_ini()

subroutine cnloc_mat104_ini ( integer nel,
integer ipg,
integer ipt,
integer nuparam,
integer nuvar,
intent(in) uparam,
intent(inout) uvar,
intent(inout) pla,
intent(inout) off,
intent(inout) thkly,
intent(inout) offg,
intent(in) sigoxx,
intent(in) sigoyy,
intent(in) sigoxy,
intent(in) sigoyz,
intent(in) sigozx,
intent(inout) thk,
intent(inout) dmg,
integer nptr,
integer npts,
integer nptt,
type(buf_lay_) bufly,
time,
intent(inout) varnl,
logical failure )

Definition at line 29 of file cnloc_mat104_ini.F.

35 USE elbufdef_mod
36 !=======================================================================
37 ! Implicit types
38 !=======================================================================
39#include "implicit_f.inc"
40 !=======================================================================
41 ! Common
42 !=======================================================================
43 !=======================================================================
44 ! Dummy arguments
45 !=======================================================================
46 INTEGER NEL,IPG,IPT,NUPARAM,NUVAR,NPTR,NPTS,NPTT,IR,IS,IT
47 my_real
48 . time
49 my_real,DIMENSION(NUPARAM), INTENT(IN) ::
50 . uparam
51 my_real,DIMENSION(NEL), INTENT(IN) ::
52 . sigoxx,sigoyy,sigoxy,sigoyz,sigozx
53 my_real ,DIMENSION(NEL), INTENT(INOUT) ::
54 . pla,off,offg,thk,varnl,dmg,thkly
55 my_real ,DIMENSION(NEL,NUVAR), INTENT(INOUT) ::
56 . uvar
57 TYPE(BUF_LAY_) :: BUFLY
58 LOGICAL :: FAILURE
59c
60 !=======================================================================
61 ! Local Variables
62 !=======================================================================
63 INTEGER I,K,II,IGURSON,ITER,NITER,NINDX,INDEX(NEL),NICE
64c
65 my_real
66 . cdr,kdr,hard,yld0,qvoce,bvoce,
67 . q1,q2,ed,an,epn,kw,fr,fc,f0,nnu
68 my_real
69 . sigvm,yld2i,omega,fcosh,fsinh,dsdrdj2,dsdrdj3,normsig,
70 . dj3dsxx,dj3dsyy,dj3dsxy,dj3dszz,dj3dsyz,dj3dszx,
71 . normxx,normyy,normxy,normzz,normyz,normzx,
72 . sdpla,dphi_dtrsig,sig_dfdsig,dphi_dsig,dphi_dfdr
73 my_real, DIMENSION(NEL) ::
74 . trsig,sxx,syy,sxy,szz,syz,szx,sigm,j2,j3,sigdr,
75 . yld,trdfds,fdr,d,triax,ft,fs,fg,fn,fsh,dezz,dlam,
76 . dlam_nl,plap_nl,dpla_nl,pla_nl,etat
77c
78 !=======================================================================
79 ! INITIALIZATION OF
80 ! DRUCKER MATERIAL LAW WITH GURSON DAMAGE
81 ! USING NON LOCAL PEERLINGS METHOD
82 !=======================================================================
83 !UVAR(5) FG DAMAGE OF CAVITIES GROWTH
84 !UVAR(6) FN DAMAGE OF CAVITIES NUCLEATION
85 !UVAR(7) FSH DAMAGE OF SHEAR
86 !UVAR(8) FT TOTAL DAMAGE (FG + FN + FSH)
87 !UVAR(9) FS EFFECTIVE DAMAGE FS = F(FT)
88 !UVAR(10) PLA_NL NON-LOCAL PLASTIC STRAIN
89 !-----------------------------------------------------------------------
90 !Warning VARNL(NEL) :
91 !input => VARNL = non-local plastic strain increment
92 !output => VARNL = cumulated local plastic strain
93 !=======================================================================
94c
95 ! Recovering model parameters
96 ! Elastic parameters
97 nnu = uparam(7)
98 ! Plastic criterion and hardening parameters [Drucker, 1948]
99 nice = nint(uparam(11)) ! Choice of the Nice method
100 cdr = uparam(12) ! Drucker coefficient
101 kdr = uparam(13) ! Drucker 1/K coefficient
102 hard = uparam(15) ! Linear hardening
103 yld0 = uparam(16) ! Initial yield stress
104 qvoce = uparam(17) ! 1st Voce parameter
105 bvoce = uparam(18) ! 2nd Voce parameter
106c
107 ! Gurson damage model parameters parameters
108 igurson = nint(uparam(30)) ! Gurson switch flag:
109 ! = 0 => no damage model
110 ! = 1 => local damage model
111 ! = 2 => non local (Forest - micromorphic) damage model
112 ! = 3 => non local (Peerlings) damage model
113 q1 = uparam(31) ! Gurson yield criterion 1st parameter
114 q2 = uparam(32) ! Gurson yield criterion 2nd parameter
115 ed = uparam(34) ! Plastic strain trigger for nucleation
116 an = uparam(35) ! Nucleation rate
117 kw = uparam(36) ! Shear coefficient (Nahshon & Hutchinson)
118 fr = uparam(37) ! Failure void volume fraction
119 fc = uparam(38) ! Critical void volume fraction
120 f0 = uparam(39) ! Initial void volume fraction
121c
122 ! Non-local micromorphic method
123 IF (igurson == 2) THEN
124 ! Recovering internal variables
125 DO i=1,nel
126 ! User inputs
127 IF (nice == 1) THEN
128 pla_nl(i) = uvar(i,7) ! Non-local plastic strain
129 ELSE
130 pla_nl(i) = uvar(i,6) ! Non-local plastic strain
131 ENDIF
132 IF (off(i) == one) THEN
133 dpla_nl(i) = max(varnl(i),zero) ! Non-local plastic strain increment
134 pla_nl(i) = pla_nl(i) + dpla_nl(i) ! Non-local cumulated plastic strain
135 ELSE
136 dpla_nl(i) = zero
137 plap_nl(i) = zero
138 ENDIF
139 IF (nice == 1) THEN
140 uvar(i,7) = pla_nl(i)
141 ELSE
142 uvar(i,6) = pla_nl(i)
143 ENDIF
144 ! Variable to regularize (cumulated plastic strain)
145 IF ((nptr==1).AND.(npts==1)) THEN
146 varnl(i) = pla(i)
147 ELSE
148 IF (ipg == nptr*npts) THEN
149 varnl(i) = zero
150 DO ir = 1,nptr
151 DO is = 1,npts
152 varnl(i) = varnl(i) + (bufly%LBUF(ir,is,ipt)%PLA(i))/(nptr*npts)
153 ENDDO
154 ENDDO
155 ENDIF
156 ENDIF
157 ENDDO
158 ! Non-local Peerlings method
159 ELSEIF (igurson == 3) THEN
160 ! Recovering internal variables
161 DO i=1,nel
162 ! User inputs
163 IF (nice == 1) THEN
164 fg(i) = uvar(i,2) ! Growth damage
165 fn(i) = uvar(i,3) ! Nucleation damage
166 fsh(i) = uvar(i,4) ! Shear damage
167 ft(i) = uvar(i,5) ! Total damage
168 fs(i) = uvar(i,6) ! Effective damage
169 pla_nl(i) = uvar(i,7) ! Non-local plastic strain
170 ELSE
171 fg(i) = uvar(i,1) ! Growth damage
172 fn(i) = uvar(i,2) ! Nucleation damage
173 fsh(i) = uvar(i,3) ! Shear damage
174 ft(i) = uvar(i,4) ! Total damage
175 fs(i) = uvar(i,5) ! Effective damage
176 pla_nl(i) = uvar(i,6) ! Non-local plastic strain
177 ENDIF
178 ! Standard inputs
179 d(i) = q1*fs(i) ! Effective normalized damage
180 dezz(i) = zero ! Initialization of the transverse strain increment
181 IF (off(i) == one) THEN
182 dpla_nl(i) = max(varnl(i),zero) ! Non-local plastic strain increment
183 pla_nl(i) = pla_nl(i) + dpla_nl(i) ! Non-local cumulated plastic strain
184 ELSE
185 dpla_nl(i) = zero
186 plap_nl(i) = zero
187 ENDIF
188 ! Computation of the yield stress
189 yld(i) = yld0 + hard*pla(i) + qvoce*(one-exp(-bvoce*pla(i)))
190c
191 !========================================================================
192 ! NON-LOCAL VARIABLES INITIALIZATION
193 !========================================================================
194c
195 ! Previous value of Drucker equivalent stress
196 trsig(i)= sigoxx(i) + sigoyy(i)
197 sigm(i) = -trsig(i) * third
198 sxx(i) = sigoxx(i) + sigm(i)
199 syy(i) = sigoyy(i) + sigm(i)
200 szz(i) = sigm(i)
201 sxy(i) = sigoxy(i)
202 syz(i) = sigoyz(i)
203 szx(i) = sigozx(i)
204 j2(i) = half*(sxx(i)**2 + syy(i)**2 + szz(i)**2 )
205 . + sxy(i)**2 + syz(i)**2 + szx(i)**2
206 j3(i) = sxx(i) * syy(i) * szz(i)
207 . + sxy(i) * syz(i) * szx(i) * two
208 . - sxx(i) * syz(i)**2
209 . - syy(i) * szx(i)**2
210 . - szz(i) * sxy(i)**2
211c
212 fdr(i) = j2(i)**3 - cdr*(j3(i)**2)
213 IF (fdr(i) > zero) THEN
214 sigdr(i) = kdr * exp((one/six)*log(fdr(i)))
215 ELSE
216 sigdr(i) = zero
217 ENDIF
218 ! Computation of the stress triaxiality and the etaT factor
219 IF (sigdr(i)>zero) THEN
220 triax(i) = (trsig(i)*third)/sigdr(i)
221 ELSE
222 triax(i) = zero
223 ENDIF
224 IF (trsig(i)<zero) THEN
225 etat(i) = zero
226 ELSE
227 etat(i) = one
228 ENDIF
229c
230 ! Normal to the previous yield surface
231 IF (yld(i)>zero) THEN
232 yld2i = one / yld(i)**2
233 dphi_dsig = two * sigdr(i) * yld2i
234 fsinh = sinh(q2*etat(i)*trsig(i)/(yld(i)*two))
235 dphi_dtrsig = q1*q2*etat(i)*fs(i)*fsinh/yld(i)
236 ELSE
237 yld2i = zero
238 dphi_dsig = zero
239 fsinh = zero
240 dphi_dtrsig = zero
241 ENDIF
242c
243 ! Computation of the Eulerian norm of the stress tensor
244 normsig = sqrt(sigoxx(i)*sigoxx(i)
245 . + sigoyy(i)*sigoyy(i)
246 . + two*sigoxy(i)*sigoxy(i)
247 . + two*sigoyz(i)*sigoyz(i)
248 . + two*sigozx(i)*sigozx(i))
249c
250 ! Computation of the norm to the yield surface
251 IF (normsig>zero) THEN
252 fdr(i) = (j2(i)/(normsig**2))**3 - cdr*((j3(i)/(normsig**3))**2)
253 dphi_dfdr = dphi_dsig*kdr*(one/six)*exp(-(five/six)*log(fdr(i)))
254 dsdrdj2 = dphi_dfdr*three*(j2(i)/(normsig**2))**2
255 dsdrdj3 = -dphi_dfdr*two*cdr*(j3(i)/(normsig**3))
256 ! dJ3/dS
257 dj3dsxx = two_third*(syy(i)*szz(i)-syz(i)**2)/(normsig**2)
258 . - third*(sxx(i)*szz(i)-szx(i)**2)/(normsig**2)
259 . - third*(sxx(i)*syy(i)-sxy(i)**2)/(normsig**2)
260 dj3dsyy = - third*(syy(i)*szz(i)-syz(i)**2)/(normsig**2)
261 . + two_third*(sxx(i)*szz(i)-szx(i)**2)/(normsig**2)
262 . - third*(sxx(i)*syy(i)-sxy(i)**2)/(normsig**2)
263 dj3dszz = - third*(syy(i)*szz(i)-syz(i)**2)/(normsig**2)
264 . - third*(sxx(i)*szz(i)-szx(i)**2)/(normsig**2)
265 . + two_third*(sxx(i)*syy(i)-sxy(i)**2)/(normsig**2)
266 dj3dsxy = two*(sxx(i)*sxy(i) + sxy(i)*syy(i) + szx(i)*syz(i))/(normsig**2)
267 dj3dsyz = two*(sxy(i)*szx(i) + syy(i)*syz(i) + syz(i)*szz(i))/(normsig**2)
268 dj3dszx = two*(sxx(i)*szx(i) + sxy(i)*syz(i) + szx(i)*szz(i))/(normsig**2)
269 ! dPhi/dSig
270 normxx = dsdrdj2*sxx(i)/normsig + dsdrdj3*dj3dsxx + dphi_dtrsig
271 normyy = dsdrdj2*syy(i)/normsig + dsdrdj3*dj3dsyy + dphi_dtrsig
272 normzz = dsdrdj2*szz(i)/normsig + dsdrdj3*dj3dszz + dphi_dtrsig
273 normxy = two*dsdrdj2*sxy(i)/normsig + dsdrdj3*dj3dsxy
274 normyz = two*dsdrdj2*syz(i)/normsig + dsdrdj3*dj3dsyz
275 normzx = two*dsdrdj2*szx(i)/normsig + dsdrdj3*dj3dszx
276 trdfds(i) = normxx + normyy + normzz
277 sig_dfdsig = sigoxx(i)*normxx + sigoyy(i)*normyy
278 . + sigoxy(i)*normxy + sigoyz(i)*normyz + sigozx(i)*normzx
279 ELSE
280 normxx = zero
281 normyy = zero
282 normzz = zero
283 normxy = zero
284 normyz = zero
285 normzx = zero
286 trdfds(i) = zero
287 sig_dfdsig = zero
288 ENDIF
289c
290 ! Computation of the non-local plastic multiplier
291 IF (sig_dfdsig>zero) THEN
292 dlam_nl(i) = ((one - ft(i))*yld(i)*dpla_nl(i))/sig_dfdsig
293 ELSE
294 dlam_nl(i) = zero
295 ENDIF
296 IF (time == zero) THEN
297 IF (sig_dfdsig>em01) THEN
298 dlam(i) = (yld(i)*pla(i))/sig_dfdsig
299 dezz(i) = -nnu*dlam(i)*(normxx + normyy) - dlam(i)*normzz
300 ELSE
301 dlam(i) = zero
302 ENDIF
303 ENDIF
304 IF ((sig_dfdsig>em01).AND.(dlam_nl(i)>zero)) THEN
305 dezz(i) = dezz(i) + nnu*dlam_nl(i)*(normxx + normyy) + dlam_nl(i)*normzz
306 ENDIF
307c
308 ! Damage growth update
309 IF ((ft(i)>zero).AND.(ft(i)<fr).AND.(trdfds(i)>zero)) THEN
310 fg(i) = fg(i) + (one-ft(i))*dlam_nl(i)*trdfds(i)
311 ENDIF
312 fg(i) = max(fg(i),zero)
313c
314 ! Nucleation damage update
315 IF ((pla_nl(i) >= ed).AND.(ft(i)<fr)) THEN
316 ! Case for positive stress triaxiality
317 IF (triax(i)>=zero) THEN
318 fn(i) = fn(i) + an*dpla_nl(i)
319 ! Case for negative stress triaxiality
320 ELSEIF ((triax(i)<zero).AND.(triax(i)>=-third)) THEN
321 fn(i) = fn(i) + an*max(one + three*triax(i),zero)*dpla_nl(i)
322 ENDIF
323 ENDIF
324 fn(i) = max(fn(i),zero)
325c
326 ! Shear damage update
327 IF ((sigdr(i) > zero).AND.(ft(i)>zero).AND.(ft(i)<fr)) THEN
328 sigvm = sqrt(max(em20,three*(j2(i)/(normsig**2))))
329 omega = one - ((twenty7 *(j3(i)/(normsig**3)))/(two*(sigvm**3)))**2
330 omega = max(omega,zero)
331 omega = min(omega,one)
332 sdpla = (sxx(i)*normxx+syy(i)*normyy+ szz(i)*normzz
333 . + sxy(i)*normxy+syz(i)*normyz+ szx(i)*normzx)
334 . * dlam_nl(i)
335 fsh(i) = fsh(i) + kw*omega*ft(i)*(sdpla/sigdr(i))
336 ENDIF
337 fsh(i) = max(fsh(i),zero)
338c
339 ! Total damage update
340 ft(i) = f0 + fg(i) + fn(i) + fsh(i)
341 ft(i) = min(ft(i),fr)
342 IF (ft(i) >= fr) THEN
343 IF (off(i)==one) off(i) = zero
344 IF (.NOT.failure) failure = .true.
345 ENDIF
346c
347 ! Effective update
348 IF (ft(i) < fc)THEN
349 fs(i) = ft(i)
350 ELSE
351 fs(i) = fc + (one/q1 - fc) * (ft(i)-fc)/(fr-fc)
352 ENDIF
353 fs(i) = min(fs(i),one/q1)
354 d(i) = q1*fs(i)
355c
356 ! Storing new values
357 ! USR Outputs
358 IF (nice == 1) THEN
359 uvar(i,2) = fg(i) ! Growth damage
360 uvar(i,3) = fn(i) ! Nucleation damage
361 uvar(i,4) = fsh(i) ! Shear damage
362 uvar(i,5) = min(ft(i),fr) ! Total damage
363 uvar(i,6) = min(fs(i),one/q1) ! Effective damage
364 uvar(i,7) = pla_nl(i) ! Non-local cumulated plastic strain
365 ELSE
366 uvar(i,1) = fg(i) ! Growth damage
367 uvar(i,2) = fn(i) ! Nucleation damage
368 uvar(i,3) = fsh(i) ! Shear damage
369 uvar(i,4) = min(ft(i),fr) ! Total damage
370 uvar(i,5) = min(fs(i),one/q1) ! Effective damage
371 uvar(i,6) = pla_nl(i) ! Non-local cumulated plastic strain
372 ENDIF
373 ! Standard outputs
374 dmg(i) = ft(i)/fr ! Normalized total damage
375 ! Variable to regularize (cumulated plastic strain)
376 IF ((nptr==1).AND.(npts==1)) THEN
377 varnl(i) = pla(i)
378 ! Computation of the thickness variation
379 thk(i) = thk(i) + dezz(i)*thk(i)*off(i)
380 ELSE
381 IF (ipg == 1) THEN
382 DO ir = 1,nptr
383 DO is = 1,npts
384 bufly%LBUF(ir,is,ipt)%THK(i) = thk(i)
385 ENDDO
386 ENDDO
387 ENDIF
388 ! Computation of the thickness variation
389 thkly(i) = thkly(i) + dezz(i)*thkly(i)*off(i)
390 IF (ipg == nptr*npts) THEN
391 varnl(i) = zero
392 thk(i) = zero
393 DO ir = 1,nptr
394 DO is = 1,npts
395 varnl(i) = varnl(i) + (bufly%LBUF(ir,is,ipt)%PLA(i))/(nptr*npts)
396 thk(i) = thk(i) + (bufly%LBUF(ir,is,ipt)%THK(i))/(nptr*npts)
397 ENDDO
398 ENDDO
399 ENDIF
400 ENDIF
401 ENDDO
402c
403 DO i = 1,nel
404 ! Checking element deletion
405 ! Case of under-integrated shells
406 IF ((nptr == 1).AND.(npts == 1)) THEN
407 IF (nptt>1) THEN
408 !Initialization for checking complete failure of the shell (all integration points)
409 IF (ipt == 1) THEN
410 offg(i) = zero
411 ENDIF
412 !If one integration points is not fully broken, the shell remains
413 IF (off(i)>zero) offg(i) = one
414 ENDIF
415 ! Case of fully integrated shells
416 ELSE
417 IF ((ipg == 1).AND.(ipt == 1)) THEN
418 !Initialization for checking complete failure of the shell (all integration points)
419 offg(i) = zero
420 ! Loop over all Gauss points (thickness + surface)
421 DO ir = 1,nptr
422 DO is = 1,npts
423 DO it = 1,nptt
424 !If one integration points is not fully broken, the shell remains
425 IF (bufly%LBUF(ir,is,it)%OFF(i)>zero) offg(i) = one
426 ENDDO
427 ENDDO
428 ENDDO
429 ENDIF
430 ENDIF
431 ENDDO
432 ENDIF
433c-----------
my_real
#define my_real
Definition cppsort.cpp:32
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21