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

Go to the source code of this file.

Functions/Subroutines

subroutine asso_qplas76c (nel, nuparam, nuvar, nfunc, ifunc, npf, tf, numtabl, table, time, timestep, uparam, uvar, rho, depsxx, depsyy, depsxy, depsyz, depszx, sigoxx, sigoyy, sigoxy, sigoyz, sigozx, signxx, signyy, signxy, signyz, signzx, pla, dpla, epsd, off, gs, yld, soundsp, dezz, inloc, dplanl, nvartmp, vartmp, loff)

Function/Subroutine Documentation

◆ asso_qplas76c()

subroutine asso_qplas76c ( integer, intent(in) nel,
integer, intent(in) nuparam,
integer, intent(in) nuvar,
integer, intent(in) nfunc,
integer, dimension(nfunc), intent(in) ifunc,
integer, dimension(snpc) npf,
tf,
integer, intent(in) numtabl,
type(table_4d_), dimension(numtabl), intent(in) table,
intent(in) time,
intent(in) timestep,
intent(in) uparam,
intent(inout) uvar,
intent(in) rho,
intent(in) depsxx,
intent(in) depsyy,
intent(in) depsxy,
intent(in) depsyz,
intent(in) depszx,
intent(in) sigoxx,
intent(in) sigoyy,
intent(in) sigoxy,
intent(in) sigoyz,
intent(in) sigozx,
intent(out) signxx,
intent(out) signyy,
intent(out) signxy,
intent(out) signyz,
intent(out) signzx,
intent(inout) pla,
intent(out) dpla,
intent(inout) epsd,
intent(inout) off,
intent(in) gs,
intent(inout) yld,
intent(inout) soundsp,
intent(out) dezz,
integer, intent(in) inloc,
intent(in) dplanl,
integer, intent(in) nvartmp,
integer, dimension(nel,nvartmp), intent(inout) vartmp,
intent(in) loff )

Definition at line 34 of file asso_qplas76c.F.

44C-----------------------------------------------
45C M o d u l e s
46C-----------------------------------------------
47 USE table4d_mod
48 USE table_mat_vinterp_mod
49C-----------------------------------------------
50C I m p l i c i t T y p e s
51C-----------------------------------------------
52#include "implicit_f.inc"
53C-----------------------------------------------
54C G l o b a l P a r a m e t e r s
55C-----------------------------------------------
56#include "tabsiz_c.inc"
57C-----------------------------------------------
58C I N P U T A r g u m e n t s
59C-----------------------------------------------
60 INTEGER, INTENT(IN) :: NEL,NUPARAM,NUVAR,NFUNC,NUMTABL,INLOC
61 INTEGER, INTENT(IN) :: NVARTMP
62 INTEGER,DIMENSION(NFUNC),INTENT(IN) :: IFUNC
63 my_real, INTENT(IN) :: time,timestep
64 my_real, DIMENSION(NUPARAM), INTENT(IN) :: uparam
65 my_real, DIMENSION(NEL), INTENT(IN) :: gs,rho,dplanl,
66 . depsxx,depsyy,depsxy,depsyz,depszx,
67 . sigoxx,sigoyy,sigoxy,sigoyz,sigozx
68 TYPE(TABLE_4D_), DIMENSION(NUMTABL) ,INTENT(IN) :: TABLE
69 my_real, DIMENSION(NEL), INTENT(IN) :: loff
70C-----------------------------------------------
71C O U T P U T A r g u m e n t s
72C-----------------------------------------------
73 my_real, DIMENSION(NEL), INTENT(OUT) ::
74 . signxx,signyy,signxy,signyz,signzx,dezz,dpla
75C-----------------------------------------------
76C I N P U T O U T P U T A r g u m e n t s
77C-----------------------------------------------
78 INTEGER, DIMENSION(NEL,NVARTMP), INTENT(INOUT) :: VARTMP
79 my_real, DIMENSION(NEL,NUVAR), INTENT(INOUT) :: uvar
80 my_real, DIMENSION(NEL), INTENT(INOUT) :: off,yld,pla,epsd,soundsp
81C-----------------------------------------------
82C VARIABLES FOR FUNCTION INTERPOLATION
83C-----------------------------------------------
84 INTEGER NPF(SNPC)
85 my_real tf(stf)
86C-----------------------------------------------
87C L o c a l V a r i a b l e s
88C-----------------------------------------------
89 INTEGER I,II,ITER,NITER,ICONV,ICAS,NINDX
90 INTEGER ,PARAMETER :: FUNC_TRAC = 1
91 INTEGER ,PARAMETER :: FUNC_COMP = 2
92 INTEGER ,PARAMETER :: FUNC_SHEAR = 3
93 INTEGER, DIMENSION(NEL) :: INDX,IAD,ILEN
94 my_real :: lam,dlam,df1,da0,da1,da2,
95 . ca,cb,aa,bb,cc,a1s,a1c,a1t,a2s,a2c,a2t,e,nu,nu1,nupc,xfac,
96 . yy,dx2,norm,dfdsigdlam,yld_norm,epd_min,epd_max,dpt,dps,
97 . normxx,normyy,normxy,normzz,normyz,normzx,alpha,alphi,dtinv,
98 . epdt_min,epdt_max,epdc_min,epdc_max,epds_min,epds_max,sig_dfdsig
99 my_real, DIMENSION(NEL) :: ff,p,svm,svm2,ylds,sxx,syy,sxy,szz,
100 . dpxx,dpyy,dpxy,dpzz,sigt,sigc,sigs,dft,dfc,dfs,a11_2d,a12_2d,g,
101 . a0,a1,a2,nup,epspt,epspc,epsps,epdt,epdc,epds,dydx,
102 . epdt_f,epdc_f,epds_f,dplat,dplac,dplas,dlam_nl
103 my_real, DIMENSION(NEL,2) :: xvec
104 my_real, DIMENSION(NUMTABL) :: tfac
105 my_real, DIMENSION(NFUNC) :: yfac
106 LOGICAL :: CONV(NEL)
107 my_real, PARAMETER :: sfac = 1.05d0 ! Security factor of ICONV
108c-----------------------------------------------
109c associated plasticity with quadratic yield function
110c-----------------------------------------
111c icas ifunt | ifunc | ifuncs
112c -1 1 | 1 | 1
113c 0 1 | 0 | 0
114c 1 1 | 1 | 0
115c 2 1 | 0 | 1
116c-----------------------------------------
117c UVAR(1) : EPSPT
118c UVAR(2) : EPSPC
119c UVAR(3) : EPSPS
120c UVAR(4) : EPDT
121c UVAR(5) : EPDC
122c UVAR(6) : EPDS
123c UVAR(7) : NUP
124C=======================================================================
125c
126 !=======================================================================
127 ! - INITIALISATION OF COMPUTATION ON TIME STEP
128 !=======================================================================
129 ! Recovering model parameters
130 ! -> Elastic parameters
131 e = uparam(1)
132 nu = uparam(5)
133 ! -> Plastic parameters
134 nupc = uparam(9)
135 ! -> Flags
136 iconv = nint(uparam(15))
137 icas = uparam(17)
138 xfac = uparam(18)
139 alpha = min(one, uparam(16)*timestep)
140 epdt_min = uparam(19)
141 epdt_max = uparam(20)
142 epdc_min = uparam(21)
143 epdc_max = uparam(22)
144 epds_min = uparam(23)
145 epds_max = uparam(24)
146 tfac(1) = uparam(25)
147 tfac(2) = uparam(26)
148 tfac(3) = uparam(27)
149 yfac(1) = uparam(28)
150 yfac(2) = uparam(29)
151 a11_2d(1:nel) = uparam(2) ! E / (ONE - NU*NU)
152 a12_2d(1:nel) = uparam(3) ! AA2 * NU
153 g(1:nel) = uparam(4)
154 nu1 = nu/(one - nu) ! aa1/aa2
155 alphi = one-alpha
156 dtinv = one / max(em20, timestep)
157c
158 ! Initialize plastic Poisson ratio
159 IF (time == zero) THEN
160 nup(1:nel) = nupc
161 uvar(1:nel,7) = nupc
162 ELSE
163 nup(1:nel) = uvar(1:nel,7)
164 END IF
165c
166 ! Recovering internal variables
167 epspt(1:nel) = uvar(1:nel,1)
168 epspc(1:nel) = uvar(1:nel,2)
169 epsps(1:nel) = uvar(1:nel,3)
170 epdt_f(1:nel) = uvar(1:nel,4)
171 epdc_f(1:nel) = uvar(1:nel,5)
172 epds_f(1:nel) = uvar(1:nel,6)
173 dplat(1:nel) = zero
174 dplac(1:nel) = zero
175 dplas(1:nel) = zero
176 DO i=1,nel
177 epdt(i) = min(epdt_max, max(epdt_f(i),epdt_min)) * xfac
178 epdc(i) = min(epdc_max, max(epdc_f(i),epdc_min)) * xfac
179 epds(i) = min(epds_max, max(epds_f(i),epds_min)) * xfac
180 ENDDO
181c
182 ! Computation of yield stresses
183 xvec(1:nel,1) = epspt(1:nel)
184 xvec(1:nel,2) = epdt(1:nel)
185 CALL table_mat_vinterp(table(func_trac),nel,nel,vartmp(1,1),xvec,sigt,dft)
186 sigt(1:nel) = sigt(1:nel) * tfac(1)
187 dft(1:nel) = dft(1:nel) * tfac(1)
188 IF (table(func_comp)%NOTABLE > 0) THEN
189 xvec(1:nel,1) = epspc(1:nel)
190 xvec(1:nel,2) = epdc(1:nel)
191 CALL table_mat_vinterp(table(func_comp),nel,nel,vartmp(1,3),xvec,sigc,dfc)
192 sigc(1:nel) = sigc(1:nel) * tfac(2)
193 dfc(1:nel) = dfc(1:nel) * tfac(2)
194 END IF
195 IF (table(func_shear)%NOTABLE > 0) THEN
196 xvec(1:nel,1) = epsps(1:nel)
197 xvec(1:nel,2) = epds(1:nel)
198 CALL table_mat_vinterp(table(func_shear),nel,nel,vartmp(1,5),xvec,sigs,dfs)
199 sigs(1:nel) = sigs(1:nel) * tfac(3)
200 dfs(1:nel) = dfs(1:nel) * tfac(3)
201 END IF
202 IF (icas == 0) THEN
203 sigc(1:nel) = sigt(1:nel)
204 sigs(1:nel) = sigt(1:nel)/sqr3
205 ELSEIF (icas == 1) THEN
206 DO i=1,nel
207 sigs(i) = sqrt(sigt(i)*sigc(i)/three)
208 ENDDO
209 ENDIF
210 ! Ensured convexity
211 IF (iconv == 1) THEN
212 DO i = 1,nel
213 conv(i) = .false.
214 IF (sigs(i) < sfac*sqrt(sigc(i)*sigt(i)/three)) THEN
215 sigs(i) = sfac*sqrt(sigc(i)*sigt(i)/three)
216 conv(i) = .true.
217 ENDIF
218 ENDDO
219 ENDIF
220 DO i=1,nel
221 aa = one/sigc(i)/sigt(i)
222 a0(i) = three*sigs(i)**2
223 a1(i) = nine*sigs(i)**2*(sigc(i) - sigt(i))*aa
224 a2(i) = nine*(sigc(i)*sigt(i) - three*sigs(i)**2)*aa
225 END DO
226c
227 !========================================================================
228 ! - COMPUTATION OF TRIAL VALUES
229 !========================================================================
230 DO i=1,nel
231 ! Computation of the trial stress tensor
232 signxx(i) = sigoxx(i) + a11_2d(i)*depsxx(i) + a12_2d(i)*depsyy(i)
233 signyy(i) = sigoyy(i) + a11_2d(i)*depsyy(i) + a12_2d(i)*depsxx(i)
234 signxy(i) = sigoxy(i) + depsxy(i)*g(i)
235 signyz(i) = sigoyz(i) + depsyz(i)*gs(i)
236 signzx(i) = sigozx(i) + depszx(i)*gs(i)
237 p(i) = -(signxx(i) + signyy(i)) * third
238 ! Computation of the deviatoric trial stress tensor
239 sxx(i) = signxx(i) + p(i)
240 syy(i) = signyy(i) + p(i)
241 szz(i) = p(i)
242 dezz(i) = -nu1 * (depsxx(i) + depsyy(i))
243 soundsp(i) = sqrt(a11_2d(i)/rho(i))
244 ENDDO
245c
246 !========================================================================
247 ! - COMPUTATION OF YIELD FONCTION
248 !========================================================================
249 ! Alpha coefficient for non associated function
250 niter = 4
251 nindx = 0
252 DO i=1,nel
253 svm2(i) = three_half*(sxx(i)**2 + syy(i)**2 + szz(i)**2) + three*signxy(i)**2
254 svm(i) = sqrt(svm2(i))
255 ylds(i) = svm2(i) - a0(i) - a1(i)*p(i) - a2(i)*p(i)*p(i)
256 IF (ylds(i) > 0 .AND. off(i) == one) THEN
257 nindx = nindx + 1
258 indx(nindx) = i
259 ENDIF
260 ENDDO
261c
262 !====================================================================
263 ! - PLASTIC RETURN MAPPING WITH CUTTING PLANE METHOD
264 !====================================================================
265 IF (nindx > 0) THEN
266c
267 ! Loop over the iterations
268 DO iter = 1,niter
269c
270 ! Loop over yielding elements
271 DO ii = 1,nindx
272 i = indx(ii)
273c
274 ! Note : in this part, the purpose is to compute for each iteration
275 ! a plastic multiplier allowing to update internal variables to satisfy
276 ! the consistency condition using the cutting plane algorithm
277 ! Its expression at each iteration is : DLAMBDA = - PHI/DPHI_DLAMBDA
278 ! -> PHI : current value of yield function (known)
279 ! -> DPHI_DLAMBDA : derivative of PHI with respect to DLAMBDA by taking
280 ! into account of internal variables kinetic :
281 ! plasticity ...
282c
283 ! df/dsig
284 cb = a1(i) + two*a2(i)*p(i)
285 norm = one/max(em20, sqrt(six*svm(i)**2 + third*cb**2)) ! norm df/dsig
286 normxx = three * sxx(i) + cb /three ! DF/DSIG
287 normyy = three * syy(i) + cb /three
288 normzz = three * szz(i) + cb /three
289 normxy = two *three * signxy(i)
290c
291 ! DF/DSIG * DSIG/DDLAM
292 dfdsigdlam = normxx * (a11_2d(i)*normxx + a12_2d(i)*normyy)*norm
293 . + normyy * (a11_2d(i)*normyy + a12_2d(i)*normxx)*norm
294 . + normxy * normxy * g(i)*norm
295c
296 yld_norm = svm(i)*norm
297 bb = three/(one + nup(i))
298 dft(i) = dft(i) * yld_norm * bb
299 IF (table(func_comp)%NOTABLE > 0) dfc(i) = dfc(i) * yld_norm * bb
300 IF (table(func_shear)%NOTABLE > 0) dfs(i) = dfs(i) * yld_norm * sqr3/two
301 IF (icas == 0) THEN
302 dfc(i) = dft(i)
303 dfs(i) = (one/sqr3)*dft(i)
304 ELSEIF (icas == 1) THEN
305 dfs(i) = (one/sqr3)*(one/(two*sqrt(sigt(i)*sigc(i))))*
306 . (dfc(i)*sigt(i) + sigc(i)*dft(i))
307 ENDIF
308 IF (iconv == 1) THEN
309 IF (conv(i)) THEN
310 dfs(i) = sfac*(one/sqr3)*(one/(two*sqrt(sigt(i)*sigc(i))))*
311 . (dfc(i)*sigt(i) + sigc(i)*dft(i))
312 ENDIF
313 ENDIF
314c
315 ! derivatives dAi/dlam
316 cc = sigs(i)/sigc(i)/sigt(i)
317C
318 a1s = eighteen*(sigc(i) - sigt(i))*cc
319 a1c = +nine*(sigs(i)/sigc(i))**2
320 a1t = -nine*(sigs(i)/sigt(i))**2
321C
322 a2s = -cinquante4*cc
323 a2c = twenty7*cc*sigs(i)/sigc(i)
324 a2t = twenty7*cc*sigs(i)/sigt(i)
325c
326 da0 = six*sigs(i)*dfs(i)
327 da1 = a1s*dfs(i) + a1t*dft(i) + a1c*dfc(i)
328 da2 = a2s*dfs(i) + a2t*dft(i) + a2c*dfc(i)
329c
330 ff(i) = dfdsigdlam + da0 + p(i)*da1 + p(i)**2 * da2
331 ff(i) = sign(max(abs(ff(i)),em20) ,ff(i))
332c
333 dlam = ylds(i)/ff(i)
334 dpla(i) = max(zero, dpla(i) + two*dlam*yld_norm )
335 pla(i) = pla(i) + two*dlam*yld_norm
336 dpt = dlam * yld_norm*bb
337 dps = dlam * yld_norm*sqr3
338 dplat(i) = dplat(i) + dpt
339 dplac(i) = dplat(i)
340 dplas(i) = dplas(i) + dps
341c
342 ! Plastic strains tensor update
343 dpxx(i) = dlam * normxx * norm
344 dpyy(i) = dlam * normyy * norm
345 dpzz(i) = dlam * normzz * norm
346 dpxy(i) = dlam * normxy * norm
347c
348 ! Elasto-plastic stresses update
349 signxx(i) = signxx(i) - (a11_2d(i)*dpxx(i) + a12_2d(i)*dpyy(i))
350 signyy(i) = signyy(i) - (a11_2d(i)*dpyy(i) + a12_2d(i)*dpxx(i))
351 signxy(i) = signxy(i) - dpxy(i)*g(i)
352c
353 ! compute EPSPC(I), EPSPT(I), EPSPS(I)
354 epspt(i) = epspt(i) + dpt
355 epspc(i) = epspc(i) + dpt
356 epsps(i) = epsps(i) + dps
357 ENDDO
358c
359 ! Update Yld values and criterion with new plastic strain and strain rate
360 xvec(1:nel,1) = epspt(1:nel)
361 xvec(1:nel,2) = epdt(1:nel)
362 CALL table_mat_vinterp(table(func_trac),nel,nel,vartmp(1,1),xvec,sigt,dft)
363 sigt(1:nel) = sigt(1:nel) * tfac(1)
364 dft(1:nel) = dft(1:nel) * tfac(1)
365 IF (table(func_comp)%NOTABLE > 0) THEN
366 xvec(1:nel,1) = epspc(1:nel)
367 xvec(1:nel,2) = epdc(1:nel)
368 CALL table_mat_vinterp(table(func_comp),nel,nel,vartmp(1,3),xvec,sigc,dfc)
369 sigc(1:nel) = sigc(1:nel) * tfac(2)
370 dfc(1:nel) = dfc(1:nel) * tfac(2)
371 END IF
372 IF (table(func_shear)%NOTABLE > 0) THEN
373 xvec(1:nel,1) = epsps(1:nel)
374 xvec(1:nel,2) = epds(1:nel)
375 CALL table_mat_vinterp(table(func_shear),nel,nel,vartmp(1,5),xvec,sigs,dfs)
376 sigs(1:nel) = sigs(1:nel) * tfac(3)
377 dfs(1:nel) = dfs(1:nel) * tfac(3)
378 END IF
379 IF (icas == 0) THEN
380 DO ii = 1,nindx
381 i = indx(ii)
382 sigc(i) = sigt(i)
383 sigs(i) = sigt(i)/sqr3
384 ENDDO
385 ELSEIF (icas == 1) THEN
386 DO ii = 1,nindx
387 i = indx(ii)
388 sigs(i) = sqrt(sigt(i)*sigc(i)/three)
389 ENDDO
390 ENDIF
391 IF (iconv == 1) THEN
392 DO ii = 1,nindx
393 i = indx(ii)
394 conv(i) = .false.
395 IF (sigs(i) < sfac*sqrt(sigc(i)*sigt(i)/three)) THEN
396 sigs(i) = sfac*sqrt(sigc(i)*sigt(i)/three)
397 conv(i) = .true.
398 ENDIF
399 END DO
400 ENDIF
401 DO ii = 1,nindx
402 i = indx(ii)
403 aa = one/sigc(i)/sigt(i)
404 a0(i) = three*sigs(i)**2
405 a1(i) = nine*sigs(i)**2*(sigc(i) - sigt(i))*aa
406 a2(i) = nine*(sigc(i)*sigt(i) - three*sigs(i)**2)*aa
407 END DO
408c
409 ! Update yield function value
410 DO ii = 1,nindx
411 i = indx(ii)
412 p(i) = -third*(signxx(i) + signyy(i) )
413 sxx(i) = signxx(i) + p(i)
414 syy(i) = signyy(i) + p(i)
415 szz(i) = p(i)
416 svm2(i)= three_half*(sxx(i)**2 + syy(i)**2 + szz(i)**2) + three*signxy(i)**2
417 svm(i) = sqrt(svm2(i))
418 ylds(i) = svm2(i) - a0(i) - a1(i)*p(i) - a2(i)*p(i)*p(i)
419 IF (inloc == 0) THEN
420 dezz(i) = dezz(i) + nu1*(dpxx(i) + dpyy(i)) + dpzz(i)
421 ENDIF
422 ENDDO
423 ENDDO ! End Newton iterations
424 END IF ! Plasticity
425c
426 !====================================================================
427 ! - UPDATE PLASTIC POISSON RATIO
428 !====================================================================
429 IF (ifunc(1) > 0) THEN
430 iad(1:nel) = npf(ifunc(1)) / 2 + 1
431 ilen(1:nel) = npf(ifunc(1)+1) / 2 - iad(1:nel) - vartmp(1:nel,8)
432!
433 CALL vinter(tf,iad,vartmp(1:nel,8),ilen,nel,pla,dydx,nup)
434!
435 uvar(1:nel,7) = yfac(1) * nup(1:nel)
436 uvar(1:nel,7) = max(zero, min(nup(1:nel), half))
437 END IF
438c
439 !====================================================================
440 ! - NON-LOCAL THICKNESS VARIATION
441 !====================================================================
442 IF (inloc > 0) THEN
443 DO i = 1,nel
444 IF (loff(i) == one) THEN
445 cb = a1(i) + two*a2(i)*p(i)
446 norm = one/max(em20, sqrt(six*svm(i)**2 + third*cb**2)) ! norm df/dsig
447 normxx = three*sxx(i) + cb/three ! DF/DSIG
448 normyy = three*syy(i) + cb/three
449 normzz = three*szz(i) + cb/three
450 yld_norm = svm(i)*norm
451 IF (yld_norm /= zero) THEN
452 dlam_nl(i) = (one/(two*yld_norm))*max(dplanl(i),zero)
453 dezz(i) = dezz(i) + nu1*(dlam_nl(i)*normxx)*norm
454 . + nu1*(dlam_nl(i)*normyy)*norm
455 . + dlam_nl(i)*normzz*norm
456 ENDIF
457 ENDIF
458 ENDDO
459 ENDIF
460c
461 !====================================================================
462 ! - STORING NEW VALUES
463 !====================================================================
464 uvar(1:nel,1) = epspt(1:nel)
465 uvar(1:nel,2) = epspc(1:nel)
466 uvar(1:nel,3) = epsps(1:nel)
467 uvar(1:nel,4) = alpha*dplat(1:nel)*dtinv + alphi*epdt_f(1:nel)
468 uvar(1:nel,5) = alpha*dplac(1:nel)*dtinv + alphi*epdc_f(1:nel)
469 uvar(1:nel,6) = alpha*dplas(1:nel)*dtinv + alphi*epds_f(1:nel)
470 epsd(1:nel) = alpha*dpla(1:nel)*dtinv + alphi*epsd(1:nel)
471c
my_real
#define my_real
Definition cppsort.cpp:32
norm
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB
alpha
#define alpha
Definition eval.h:35
func_comp
subroutine func_comp(table, ntable, nptmax, npt_trac, npt_shear, npt_comp, x_comp, y_comp, nup)
Definition law76_upd.F:344
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
table4d_mod
Definition table4d_mod.F:52
table_mat_vinterp_mod
Definition table_mat_vinterp.F:32
table_mat_vinterp_mod::table_mat_vinterp
subroutine table_mat_vinterp(table, dimx, nel, ipos, xx, yy, dydx)
Definition table_mat_vinterp.F:45
vinter
subroutine vinter(tf, iad, ipos, ilen, nel, x, dydx, y)
Definition vinter.F:72