OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
cnloc_mat104_ini.F
Go to the documentation of this file.
1Copyright> OpenRadioss
2Copyright> Copyright (C) 1986-2025 Altair Engineering Inc.
3Copyright>
4Copyright> This program is free software: you can redistribute it and/or modify
5Copyright> it under the terms of the GNU Affero General Public License as published by
6Copyright> the Free Software Foundation, either version 3 of the License, or
7Copyright> (at your option) any later version.
8Copyright>
9Copyright> This program is distributed in the hope that it will be useful,
10Copyright> but WITHOUT ANY WARRANTY; without even the implied warranty of
11Copyright> MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12Copyright> GNU Affero General Public License for more details.
13Copyright>
14Copyright> You should have received a copy of the GNU Affero General Public License
15Copyright> along with this program. If not, see <https://www.gnu.org/licenses/>.
16Copyright>
17Copyright>
18Copyright> Commercial Alternative: Altair Radioss Software
19Copyright>
20Copyright> As an alternative to this open-source version, Altair also offers Altair Radioss
21Copyright> software under a commercial license. Contact Altair to discuss further if the
22Copyright> commercial version may interest you: https://www.altair.com/radioss/.
23!||====================================================================
24!|| cnloc_mat104_ini ../starter/source/materials/mat/mat104/cnloc_mat104_ini.F
25!||--- called by ------------------------------------------------------
26!|| cnloc_matini ../starter/source/materials/mat_share/cnloc_matini.F
27!||--- uses -----------------------------------------------------
28!||====================================================================
29 SUBROUTINE cnloc_mat104_ini(
30 . NEL ,IPG ,IPT ,NUPARAM ,NUVAR ,UPARAM ,
31 . UVAR ,PLA ,OFF ,THKLY ,OFFG ,
32 . SIGOXX ,SIGOYY ,SIGOXY ,SIGOYZ ,SIGOZX ,
33 . THK ,DMG ,NPTR ,NPTS ,NPTT ,BUFLY ,
34 . TIME ,VARNL ,FAILURE )
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-----------
434 END
cnloc_mat104_ini
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)
Definition cnloc_mat104_ini.F:35
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21