Line 232: Line 232:
 
| style="text-align: center;" | <math>\int _\Omega \delta v_i r_{m_i}d\Omega  + \int _{\Gamma _t} \delta v_i (\sigma _{ij} n_j - t_i)d\Gamma + \int _{\Omega } {h_j\over 2}{\partial \delta v_i \over \partial x_j} r_{m_i} d\Omega =0 </math>
 
| style="text-align: center;" | <math>\int _\Omega \delta v_i r_{m_i}d\Omega  + \int _{\Gamma _t} \delta v_i (\sigma _{ij} n_j - t_i)d\Gamma + \int _{\Omega } {h_j\over 2}{\partial \delta v_i \over \partial x_j} r_{m_i} d\Omega =0 </math>
 
|}
 
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (15a)
+
| style="width: 5px;text-align: right;white-space: nowrap;" | (15a)
 
|}
 
|}
  
Line 242: Line 242:
 
| style="text-align: center;" | <math>\int _\Omega q r_d d\Omega + \int _\Omega \left[\sum \limits _{i=1}^{n_d}\tau _i {\partial q \over \partial x_i}r_{m_i} \right]d\Omega - \int _\Gamma \left[\sum \limits _{i=1}^{n_d} q \tau _i n_i r_{m_i}\right]d\Gamma =0</math>
 
| style="text-align: center;" | <math>\int _\Omega q r_d d\Omega + \int _\Omega \left[\sum \limits _{i=1}^{n_d}\tau _i {\partial q \over \partial x_i}r_{m_i} \right]d\Omega - \int _\Gamma \left[\sum \limits _{i=1}^{n_d} q \tau _i n_i r_{m_i}\right]d\Gamma =0</math>
 
|}
 
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (15b)
+
| style="width: 5px;text-align: right;white-space: nowrap;" | (15b)
 
|}
 
|}
  
Line 254: Line 254:
 
| style="text-align: center;" | <math>\begin{array}{r} \displaystyle \int _\Omega \left[\delta v_i\rho \left({\partial v_i \over \partial t}+v_j {\partial {v_i} \over \partial x_j}\right)+ {\partial \delta v_i \over \partial x_j}\left(\mu {\partial v_i \over \partial x_j}-\delta _{ij}p \right)\right]d\Omega -  \int _{\Omega } \delta v_i b_i d\Omega - \int _{\Gamma _t} \delta v_i t_i^p  d\Gamma +\qquad \\ \displaystyle + \int _{\Omega } {h_j\over 2}{\partial \delta v_i \over \partial x_j} r_{m_i} d\Omega =0\qquad  \end{array}</math>
 
| style="text-align: center;" | <math>\begin{array}{r} \displaystyle \int _\Omega \left[\delta v_i\rho \left({\partial v_i \over \partial t}+v_j {\partial {v_i} \over \partial x_j}\right)+ {\partial \delta v_i \over \partial x_j}\left(\mu {\partial v_i \over \partial x_j}-\delta _{ij}p \right)\right]d\Omega -  \int _{\Omega } \delta v_i b_i d\Omega - \int _{\Gamma _t} \delta v_i t_i^p  d\Gamma +\qquad \\ \displaystyle + \int _{\Omega } {h_j\over 2}{\partial \delta v_i \over \partial x_j} r_{m_i} d\Omega =0\qquad  \end{array}</math>
 
|}
 
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (16a)
+
| style="width: 5px;text-align: right;white-space: nowrap;" | (16a)
 
|}
 
|}
  
Line 264: Line 264:
 
| style="text-align: center;" | <math>\int _\Omega q {\partial v_i \over \partial x_i} d\Omega + \int _\Omega \left[\sum \limits _{i=1}^{n_d} \tau _i {\partial q \over \partial x_i} r_{m_i}\right]d\Omega =0</math>
 
| style="text-align: center;" | <math>\int _\Omega q {\partial v_i \over \partial x_i} d\Omega + \int _\Omega \left[\sum \limits _{i=1}^{n_d} \tau _i {\partial q \over \partial x_i} r_{m_i}\right]d\Omega =0</math>
 
|}
 
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (16b)
+
| style="width: 5px;text-align: right;white-space: nowrap;" | (16b)
 
|}
 
|}
  
Line 366: Line 366:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
| style="text-align: center;" | <math>\displaystyle {M}\dot{\bar {v}} + {H} \bar {v} - {G}\bar {p}+{C}\bar {c}={f}</math>
+
| style="text-align: center;" | <math>\displaystyle {\boldsymbol M}\dot{\bar {\boldsymbol v}} + {\boldsymbol H} \bar {\boldsymbol v} - {\boldsymbol G}\bar {\boldsymbol p}+{\boldsymbol C}\bar {\boldsymbol c}={\boldsymbol f}</math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (24a)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (24a)
Line 376: Line 376:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
| style="text-align: center;" | <math>\displaystyle {G}^T \bar {v} + \hat{L}\bar {p}+{Q}\bar {\boldsymbol \pi }={0}</math>
+
| style="text-align: center;" | <math>\displaystyle {\boldsymbol G}^T \bar {\boldsymbol v} + \hat{\boldsymbol L}\bar {\boldsymbol p}+{\boldsymbol Q}\bar {\boldsymbol \pi }={\boldsymbol 0}</math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (24b)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (24b)
Line 386: Line 386:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
| style="text-align: center;" | <math>\displaystyle \hat {C}\bar{v}+ {M}\bar {c}={0}</math>
+
| style="text-align: center;" | <math>\displaystyle \hat {\boldsymbol C}\bar{\boldsymbol v}+ {\boldsymbol M}\bar {\boldsymbol c}={\boldsymbol 0}</math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (24c)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (24c)
Line 396: Line 396:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
| style="text-align: center;" | <math>\displaystyle {Q}^T \bar {p} + \hat {M}\bar {\boldsymbol \pi }={0}</math>
+
| style="text-align: center;" | <math>\displaystyle {\boldsymbol Q}^T \bar {\boldsymbol p} + \hat {\boldsymbol M}\bar {\boldsymbol \pi }={\boldsymbol 0}</math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (24d)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (24d)
Line 410: Line 410:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
| style="text-align: center;" | <math>{M} \displaystyle{1\over \Delta t} (\bar{v}^{n+1}-\bar{v}^{n}) + {H}^{n+\theta }\bar{v} ^{n+\theta } - {G}\bar{p}^{n+\theta }+ {C}^{n+\theta }\bar {c}^{n+\theta }={f}^{n+\theta }</math>
+
| style="text-align: center;" | <math>{\boldsymbol M} \displaystyle{1\over \Delta t} (\bar{\boldsymbol v}^{n+1}-\bar{\boldsymbol v}^{n}) + {\boldsymbol H}^{n+\theta }\bar{\boldsymbol v} ^{n+\theta } - {\boldsymbol G}\bar{\boldsymbol p}^{n+\theta }+ {\boldsymbol C}^{n+\theta }\bar {\boldsymbol c}^{n+\theta }={\boldsymbol f}^{n+\theta }</math>
 
|}
 
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (25a)
+
| style="width: 5px;text-align: right;white-space: nowrap;" | (25a)
 
|}
 
|}
  
Line 420: Line 420:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
| style="text-align: center;" | <math>{G}^T \bar {v}^{n+\theta }+\hat {L}^{n+\theta } \bar {p}^{n+\theta }+{Q}\bar {\boldsymbol \pi }^{n+\theta }={0}</math>
+
| style="text-align: center;" | <math>{\boldsymbol G}^T \bar {\boldsymbol v}^{n+\theta }+\hat {\boldsymbol L}^{n+\theta } \bar {\boldsymbol p}^{n+\theta }+{\boldsymbol Q}\bar {\boldsymbol \pi }^{n+\theta }={\boldsymbol 0}</math>
 
|}
 
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (25b)
+
| style="width: 5px;text-align: right;white-space: nowrap;" | (25b)
 
|}
 
|}
  
Line 430: Line 430:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
| style="text-align: center;" | <math>\hat {C}^{n+\theta } \bar {v}^{n+\theta }+ {M} \bar {c}^{n+\theta }={0}</math>
+
| style="text-align: center;" | <math>\hat {\boldsymbol C}^{n+\theta } \bar {\boldsymbol v}^{n+\theta }+ {\boldsymbol M} \bar {\boldsymbol c}^{n+\theta }={\boldsymbol 0}</math>
 
|}
 
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (25c)
+
| style="width: 5px;text-align: right;white-space: nowrap;" | (25c)
 
|}
 
|}
  
Line 440: Line 440:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
| style="text-align: center;" | <math>{G}^T \bar {p}^{n+\theta }+\hat {M}^{n+\theta }\bar {\boldsymbol \pi }^{n+\theta }={0}</math>
+
| style="text-align: center;" | <math>{\boldsymbol G}^T \bar {\boldsymbol p}^{n+\theta }+\hat {\boldsymbol M}^{n+\theta }\bar {\boldsymbol \pi }^{n+\theta }={\boldsymbol 0}</math>
 
|}
 
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (25d)
+
| style="width: 5px;text-align: right;white-space: nowrap;" | (25d)
 
|}
 
|}
  
where <math display="inline">{H}^{n+\theta }={H} (\bar{v}^{n+\theta })</math>, <math display="inline">\bar{v}^{n+\theta }</math> are the velocities evaluated at time <math display="inline">n+\theta </math> and the parameter <math display="inline">\theta \in [0,1]</math>. The direct monolitic solution of Eqs.(25) is possible using an adequate iterative scheme. However, we have found more convenient to  use  a fractional step method as described in the next section.
+
where <math display="inline">{\boldsymbol H}^{n+\theta }={\boldsymbol H} (\bar{\boldsymbol v}^{n+\theta })</math>, <math display="inline">\bar{\boldsymbol v}^{n+\theta }</math> are the velocities evaluated at time <math display="inline">n+\theta </math> and the parameter <math display="inline">\theta \in [0,1]</math>. The direct monolitic solution of Eqs.(25) is possible using an adequate iterative scheme. However, we have found more convenient to  use  a fractional step method as described in the next section.
  
 
===3.1 Fractional step method===
 
===3.1 Fractional step method===
Line 456: Line 456:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
| style="text-align: center;" | <math>{M} \displaystyle{1\over \Delta t} (\tilde{v}^{n+1}-\bar{v}^{n}) + {H}^{n+\theta }\bar{v} ^{n+\theta } - \alpha {G}\bar{p}^{n} + {C}^{n+\theta }\bar {c}^{n+\theta }={f}^{n+\theta }</math>
+
| style="text-align: center;" | <math>{\boldsymbol M} \displaystyle{1\over \Delta t} (\tilde{\boldsymbol v}^{n+1}-\bar{\boldsymbol v}^{n}) + {\boldsymbol H}^{n+\theta }\bar{\boldsymbol v} ^{n+\theta } - \alpha {\boldsymbol G}\bar{\boldsymbol p}^{n} + {\boldsymbol C}^{n+\theta }\bar {\boldsymbol c}^{n+\theta }={\boldsymbol f}^{n+\theta }</math>
 
|}
 
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (26)
+
| style="width: 5px;text-align: right;white-space: nowrap;" | (26)
 
|}
 
|}
  
Line 466: Line 466:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
| style="text-align: center;" | <math>{M} \displaystyle{1\over \Delta t} (\bar{v}^{n+1}-\tilde{v}^{n+1})- {G}(\bar{p}^{n+1}-\alpha \bar {p}^{n})={0}</math>
+
| style="text-align: center;" | <math>{\boldsymbol M} \displaystyle{1\over \Delta t} (\bar{\boldsymbol v}^{n+1}-\tilde{\boldsymbol v}^{n+1})- {\boldsymbol G}(\bar{\boldsymbol p}^{n+1}-\alpha \bar {\boldsymbol p}^{n})={\boldsymbol 0}</math>
 
|}
 
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (27)
+
| style="width: 5px;text-align: right;white-space: nowrap;" | (27)
 
|}
 
|}
  
In above equation <math display="inline">\tilde{v}^{n+1}</math> is a predicted value of the velocity at time <math display="inline">n+1</math> and <math display="inline">\alpha </math> is a variable whose values of interest are zero and one. For <math display="inline">\alpha =0</math> (first order scheme) the splitting error is of order <math display="inline"> 0 (\Delta t)</math>, whereas for <math display="inline">\alpha =1</math> (second order scheme) the error is of order <math display="inline">0 (\Delta t^2)</math> [45].
+
In above equation <math display="inline">\tilde{\boldsymbol v}^{n+1}</math> is a predicted value of the velocity at time <math display="inline">n+1</math> and <math display="inline">\alpha </math> is a variable whose values of interest are zero and one. For <math display="inline">\alpha =0</math> (first order scheme) the splitting error is of order <math display="inline"> 0 (\Delta t)</math>, whereas for <math display="inline">\alpha =1</math> (second order scheme) the error is of order <math display="inline">0 (\Delta t^2)</math> [45].
  
 
Eqs.(26) and (27) are completed with the following three equations emanating from Eqs.(25b-d)
 
Eqs.(26) and (27) are completed with the following three equations emanating from Eqs.(25b-d)
Line 480: Line 480:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
| style="text-align: center;" | <math>{G}^T\bar{v}^{n+1}+\hat {L}^n \bar{p}^{n+1}+{Q}\bar {\boldsymbol \pi }^{n}= {0}</math>
+
| style="text-align: center;" | <math>{\boldsymbol G}^T\bar{\boldsymbol v}^{n+1}+\hat {\boldsymbol L}^n \bar{\boldsymbol p}^{n+1}+{\boldsymbol Q}\bar {\boldsymbol \pi }^{n}= {\boldsymbol 0}</math>
 
|}
 
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (28a)
+
| style="width: 5px;text-align: right;white-space: nowrap;" | (28a)
 
|}
 
|}
  
Line 490: Line 490:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
| style="text-align: center;" | <math>\hat {C}^{n+1} \bar{v}^{n+1} + {M} \bar {c}^{n+1}= {0}</math>
+
| style="text-align: center;" | <math>\hat {\boldsymbol C}^{n+1} \bar{\boldsymbol v}^{n+1} + {\boldsymbol M} \bar {\boldsymbol c}^{n+1}= {\boldsymbol 0}</math>
 
|}
 
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (28b)
+
| style="width: 5px;text-align: right;white-space: nowrap;" | (28b)
 
|}
 
|}
  
Line 500: Line 500:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
| style="text-align: center;" | <math>{Q}^T \bar{p}^{n+1}+ \hat {M}^{n+1} \bar {\boldsymbol \pi }^{n+1}= {0}</math>
+
| style="text-align: center;" | <math>{\boldsymbol Q}^T \bar{\boldsymbol p}^{n+1}+ \hat {\boldsymbol M}^{n+1} \bar {\boldsymbol \pi }^{n+1}= {\boldsymbol 0}</math>
 
|}
 
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (28c)
+
| style="width: 5px;text-align: right;white-space: nowrap;" | (28c)
 
|}
 
|}
  
The value of <math display="inline">\bar{v}^{n+1}</math> obtained from Eq.(27) is substituted into Eq.(28a) to give
+
The value of <math display="inline">\bar{\boldsymbol v}^{n+1}</math> obtained from Eq.(27) is substituted into Eq.(28a) to give
  
 
{| class="formulaSCP" style="width: 100%; text-align: left;"  
 
{| class="formulaSCP" style="width: 100%; text-align: left;"  
Line 512: Line 512:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
| style="text-align: center;" | <math>{G}^T\tilde{v}^{n+1}+ \Delta t {G}^T {M}^{-1} {G} (\bar{p}^{n+1} - \alpha \bar{p}^n)+ \hat {L}^n {p}^{n+1}+{Q}\bar {\boldsymbol \pi }^n={0} </math>
+
| style="text-align: center;" | <math>{\boldsymbol G}^T\tilde{\boldsymbol v}^{n+1}+ \Delta t {\boldsymbol G}^T {\boldsymbol M}^{-1} {\boldsymbol G} (\bar{\boldsymbol p}^{n+1} - \alpha \bar{\boldsymbol p}^n)+ \hat {\boldsymbol L}^n {\boldsymbol p}^{n+1}+{\boldsymbol Q}\bar {\boldsymbol \pi }^n={\boldsymbol 0} </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (29)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (29)
 
|}
 
|}
  
The product <math display="inline">{G}^T {M}^{-1} {G}</math> can be approximated by a laplacian matrix, i.e.
+
The product <math display="inline">{\boldsymbol G}^T {\boldsymbol M}^{-1} {\boldsymbol G}</math> can be approximated by a laplacian matrix, i.e.
  
 
{| class="formulaSCP" style="width: 100%; text-align: left;"  
 
{| class="formulaSCP" style="width: 100%; text-align: left;"  
Line 524: Line 524:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
| style="text-align: center;" | <math>{G}^T {M}^{-1} {G}={1\over \rho } {L}\quad \hbox{with }L^{ab}=\int _{ \Omega ^e} {\boldsymbol \nabla }^T N^a {\boldsymbol \nabla } N^b d\Omega  </math>
+
| style="text-align: center;" | <math>{\boldsymbol G}^T {\boldsymbol M}^{-1} {\boldsymbol G}={1\over \rho } {\boldsymbol L}\quad \hbox{with }L^{ab}=\int _{ \Omega ^e} {\boldsymbol \nabla }^T N^a {\boldsymbol \nabla } N^b d\Omega  </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (30)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (30)
 
|}
 
|}
  
where <math display="inline">L^{ab}</math> are the element contributions to <math display="inline">{L}</math> (see Appendix).
+
where <math display="inline">L^{ab}</math> are the element contributions to <math display="inline">{\boldsymbol L}</math> (see Appendix).
  
 
The steps of the fractional step scheme chosen here are:
 
The steps of the fractional step scheme chosen here are:
Line 535: Line 535:
 
==Step 1==
 
==Step 1==
  
The fractional nodal velocities <math display="inline">\tilde {v}^{n+1}</math> can be  explicitely computed from Eq.(26) by
+
The fractional nodal velocities <math display="inline">\tilde {\boldsymbol v}^{n+1}</math> can be  explicitely computed from Eq.(26) by
  
 
{| class="formulaSCP" style="width: 100%; text-align: left;"  
 
{| class="formulaSCP" style="width: 100%; text-align: left;"  
Line 542: Line 542:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
| style="text-align: center;" | <math>\tilde{v}^{n+1} = \bar{v}^{n} - \Delta t  {M}^{-1}_d [\tilde{H}^{n}\bar{v}^{n} - \alpha {G} \bar{p}^n + {C}^{n}\bar{c}^{n} - {f}^{n}]</math>
+
| style="text-align: center;" | <math>\tilde{\boldsymbol v}^{n+1} = \bar{\boldsymbol v}^{n} - \Delta t  {\boldsymbol M}^{-1}_d [\tilde{\boldsymbol H}^{n}\bar{\boldsymbol v}^{n} - \alpha {\boldsymbol G} \bar{\boldsymbol p}^n + {\boldsymbol C}^{n}\bar{\boldsymbol c}^{n} - {\boldsymbol f}^{n}]</math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (31)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (31)
 
|}
 
|}
  
where <math display="inline">{M}_d</math> is the diagonal form of '''M''' obtaining by lumping the row terms into the corresponding diagonal terms.
+
where <math display="inline">{\boldsymbol M}_d</math> is the diagonal form of '''M''' obtaining by lumping the row terms into the corresponding diagonal terms.
  
 
<br/>
 
<br/>
  
'''Step 2''' Compute <math display="inline">\bar{p}^{n+1}</math> from Eq.(29) as
+
'''Step 2''' Compute <math display="inline">\bar{\boldsymbol p}^{n+1}</math> from Eq.(29) as
  
 
{| class="formulaSCP" style="width: 100%; text-align: left;"  
 
{| class="formulaSCP" style="width: 100%; text-align: left;"  
Line 558: Line 558:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
| style="text-align: center;" | <math>\bar{p}^{n+1}= -[\hat{L}^n + {\Delta t \over \rho } {L}]^{-1} [{G}^T\tilde{v}^{n+1} - \alpha {\Delta t \over \rho }{L}\bar{p}^n +{Q} \bar{\boldsymbol \pi }^{n}] </math>
+
| style="text-align: center;" | <math>\bar{\boldsymbol p}^{n+1}= -[\hat{\boldsymbol L}^n + {\Delta t \over \rho } {\boldsymbol L}]^{-1} [{\boldsymbol G}^T\tilde{\boldsymbol v}^{n+1} - \alpha {\Delta t \over \rho }{\boldsymbol L}\bar{\boldsymbol p}^n +{\boldsymbol Q} \bar{\boldsymbol \pi }^{n}] </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (32)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (32)
 
|}
 
|}
  
'''Step 3''' Compute <math display="inline"> \bar{v}^{n+1}</math> explicitely from Eq.(28a) as
+
'''Step 3''' Compute <math display="inline"> \bar{\boldsymbol v}^{n+1}</math> explicitely from Eq.(28a) as
  
 
{| class="formulaSCP" style="width: 100%; text-align: left;"  
 
{| class="formulaSCP" style="width: 100%; text-align: left;"  
Line 570: Line 570:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
| style="text-align: center;" | <math>\bar {v}^{n+1}=\tilde{v}^{n+1}+ \Delta t  {M}_d^{-1} {G} (\bar {p}^{n+1}- \alpha \bar {p}^n) </math>
+
| style="text-align: center;" | <math>\bar {\boldsymbol v}^{n+1}=\tilde{\boldsymbol v}^{n+1}+ \Delta t  {\boldsymbol M}_d^{-1} {\boldsymbol G} (\bar {\boldsymbol p}^{n+1}- \alpha \bar {\boldsymbol p}^n) </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (33)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (33)
Line 577: Line 577:
 
<br/>
 
<br/>
  
'''Step 4''' Compute <math display="inline"> \bar{c}^{n+1}</math> explicitely from Eq.(28b) as
+
'''Step 4''' Compute <math display="inline"> \bar{\boldsymbol c}^{n+1}</math> explicitely from Eq.(28b) as
  
 
{| class="formulaSCP" style="width: 100%; text-align: left;"  
 
{| class="formulaSCP" style="width: 100%; text-align: left;"  
Line 584: Line 584:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
| style="text-align: center;" | <math>\bar{c}^{n+1}=- {M}_d^{-1}\hat {C}^{n+1}\bar{v}^{n+1} </math>
+
| style="text-align: center;" | <math>\bar{\boldsymbol c}^{n+1}=- {\boldsymbol M}_d^{-1}\hat {\boldsymbol C}^{n+1}\bar{\boldsymbol v}^{n+1} </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (34)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (34)
Line 596: Line 596:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
| style="text-align: center;" | <math>\bar{\boldsymbol \pi }^{n+1}=- \hat {M}_d^{-1} {Q}^T \bar {p}^{n+1} </math>
+
| style="text-align: center;" | <math>\bar{\boldsymbol \pi }^{n+1}=- \hat {\boldsymbol M}_d^{-1} {\boldsymbol Q}^T \bar {\boldsymbol p}^{n+1} </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (35)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (35)
 
|}
 
|}
  
where <math display="inline">\hat {M}_d</math> is the lumped form of <math display="inline">\hat {M}</math>. A standard diagonal lumping procedure based in summing up the terms of each row has been used.
+
where <math display="inline">\hat {\boldsymbol M}_d</math> is the lumped form of <math display="inline">\hat {\boldsymbol M}</math>. A standard diagonal lumping procedure based in summing up the terms of each row has been used.
  
 
Note that all steps can be solved explicitely except for the computation of the pressure in Eq.(32) which requires to invert the sum of two laplacian matrices. This can be effectively performed using an iterative solution scheme such as the conjugate gradient method.
 
Note that all steps can be solved explicitely except for the computation of the pressure in Eq.(32) which requires to invert the sum of two laplacian matrices. This can be effectively performed using an iterative solution scheme such as the conjugate gradient method.
  
Above algorithm has improved stabilization properties versus the standard segregation methods due to the introduction of the laplacian matrix <math display="inline">\hat {L}</math> in Eq.(32). This matrix emanates from the FIC formulation.
+
Above algorithm has improved stabilization properties versus the standard segregation methods due to the introduction of the laplacian matrix <math display="inline">\hat {\boldsymbol L}</math> in Eq.(32). This matrix emanates from the FIC formulation.
  
The boundary conditions are applied as follows. No condition is applied in the computation of the fractional velocities <math display="inline">\tilde{v}^{n+1}</math> in Eq.(31). The prescribed velocities at the boundary are applied when solving for <math display="inline">\bar{v}^{n+1}</math> in  step 3. The prescribed pressures at the boundary are imposed by making <math display="inline">\bar{p}^n</math> equal to the prescribed pressure values when solving Eq.(32).
+
The boundary conditions are applied as follows. No condition is applied in the computation of the fractional velocities <math display="inline">\tilde{\boldsymbol v}^{n+1}</math> in Eq.(31). The prescribed velocities at the boundary are applied when solving for <math display="inline">\bar{\boldsymbol v}^{n+1}</math> in  step 3. The prescribed pressures at the boundary are imposed by making <math display="inline">\bar{\boldsymbol p}^n</math> equal to the prescribed pressure values when solving Eq.(32).
  
 
===3.2 Stokes flow===
 
===3.2 Stokes flow===
Line 632: Line 632:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
| style="text-align: center;" | <math>\begin{array}{l}\displaystyle {M}\dot{\bar {v}} + {K}\bar{v} - {G}\bar {p}={f}\\  \displaystyle {G}^T \bar {v} + \hat{L}\bar {p}+{Q}\bar {\boldsymbol \pi }={0}\\ \displaystyle {Q}^T \bar {p} + \hat {M}\bar {\boldsymbol \pi }={0} \end{array} </math>
+
| style="text-align: center;" | <math>\begin{array}{l}\displaystyle {\boldsymbol M}\dot{\bar {\boldsymbol v}} + {\boldsymbol K}\bar{\boldsymbol v} - {\boldsymbol G}\bar {\boldsymbol p}={\boldsymbol f}\\  \displaystyle {\boldsymbol G}^T \bar {\boldsymbol v} + \hat{\boldsymbol L}\bar {\boldsymbol p}+{\boldsymbol Q}\bar {\boldsymbol \pi }={\boldsymbol 0}\\ \displaystyle {\boldsymbol Q}^T \bar {\boldsymbol p} + \hat {\boldsymbol M}\bar {\boldsymbol \pi }={\boldsymbol 0} \end{array} </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (37)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (37)
Line 646: Line 646:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
| style="text-align: center;" | <math>\left[\begin{matrix}{K}&-{G}&{0}\\ -{G}^T & - \hat{L} &-{Q}\\ {0}& -{Q}^T&-\hat {M}\\\end{matrix}\right]\left\{\begin{matrix}\bar {v}\\ \bar {p}\\ \bar {\boldsymbol \pi }\\\end{matrix}\right\}= \left\{\begin{matrix}{f}\\ {0}\\{0}\\ \end{matrix}\right\} </math>
+
| style="text-align: center;" | <math>\left[\begin{matrix}{\boldsymbol K}&-{\boldsymbol G}&{\boldsymbol 0}\\ -{\boldsymbol G}^T & - \hat{\boldsymbol L} &-{\boldsymbol Q}\\ {\boldsymbol 0}& -{\boldsymbol Q}^T&-\hat {\boldsymbol M}\\\end{matrix}\right]\left\{\begin{matrix}\bar {\boldsymbol v}\\ \bar {\boldsymbol p}\\ \bar {\boldsymbol \pi }\\\end{matrix}\right\}= \left\{\begin{matrix}{\boldsymbol f}\\ {\boldsymbol 0}\\{\boldsymbol 0}\\ \end{matrix}\right\} </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (38)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (38)
 
|}
 
|}
  
The system is symmetric and always positive definite and therefore leads to a non singular solution. This property holds for ''any interpolation function'' chosen for <math display="inline">\bar {v},\bar {p}</math> and <math display="inline">\bar {\boldsymbol \pi }</math>, therefore overcoming the Babuŝka-Brezzi (BB) restrictions [8].
+
The system is symmetric and always positive definite and therefore leads to a non singular solution. This property holds for ''any interpolation function'' chosen for <math display="inline">\bar {\boldsymbol v},\bar {\boldsymbol p}</math> and <math display="inline">\bar {\boldsymbol \pi }</math>, therefore overcoming the Babuŝka-Brezzi (BB) restrictions [8].
  
 
A reduced velocity-pressure formulation can be obtained by elliminating the <math display="inline">\bar \pi </math> variables from the last row of Eq.(38)  [37].
 
A reduced velocity-pressure formulation can be obtained by elliminating the <math display="inline">\bar \pi </math> variables from the last row of Eq.(38)  [37].
Line 745: Line 745:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
| style="text-align: center;" | <math>{h}=h_s {{v}\over {v}}+h_{c} {{\boldsymbol \nabla } v\over \vert{\boldsymbol \nabla }v\vert }</math>
+
| style="text-align: center;" | <math>{\boldsymbol h}=h_s {{\boldsymbol v}\over {v}}+h_{c} {{\boldsymbol \nabla } v\over \vert{\boldsymbol \nabla }v\vert }</math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (41a)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (41a)
 
|}
 
|}
  
where <math display="inline">v=\vert {v}\vert </math> and <math display="inline">h_s</math> and <math display="inline">h_{c}</math> are the “streamline” and “cross wind” contributions given by
+
where <math display="inline">v=\vert {\boldsymbol v}\vert </math> and <math display="inline">h_s</math> and <math display="inline">h_{c}</math> are the “streamline” and “cross wind” contributions given by
  
 
{| class="formulaSCP" style="width: 100%; text-align: left;"  
 
{| class="formulaSCP" style="width: 100%; text-align: left;"  
Line 757: Line 757:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
| style="text-align: center;" | <math>h_s=\max ({l}^T_j {v})/{v} </math>
+
| style="text-align: center;" | <math>h_s=\max ({\boldsymbol l}^T_j {\boldsymbol v})/{v} </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (41b)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (41b)
Line 767: Line 767:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
| style="text-align: center;" | <math>h_{c}=\max ({l}^T_j {\boldsymbol \nabla }v)/ \vert {\boldsymbol \nabla }v\vert \quad , \quad  j=1,n_s</math>
+
| style="text-align: center;" | <math>h_{c}=\max ({\boldsymbol l}^T_j {\boldsymbol \nabla }v)/ \vert {\boldsymbol \nabla }v\vert \quad , \quad  j=1,n_s</math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (41c)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (41c)
 
|}
 
|}
  
where <math display="inline">{l}_j</math> are the vectors defining the element sides (<math display="inline">n_s=6</math> for tetrahedra).
+
where <math display="inline">{\boldsymbol l}_j</math> are the vectors defining the element sides (<math display="inline">n_s=6</math> for tetrahedra).
  
As for the free surface equation the following value of the characteristic length vector <math display="inline">\bar{h}</math> in Eq.(40a) has been taken
+
As for the free surface equation the following value of the characteristic length vector <math display="inline">\bar{\boldsymbol h}</math> in Eq.(40a) has been taken
  
 
{| class="formulaSCP" style="width: 100%; text-align: left;"  
 
{| class="formulaSCP" style="width: 100%; text-align: left;"  
Line 781: Line 781:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
| style="text-align: center;" | <math>\bar{h} =\bar h_s {{v}\over {v}}+\bar h_c {{\boldsymbol \nabla }\beta \over \vert {\boldsymbol \nabla }\beta \vert } </math>
+
| style="text-align: center;" | <math>\bar{\boldsymbol h} =\bar h_s {{\boldsymbol v}\over {v}}+\bar h_c {{\boldsymbol \nabla }\beta \over \vert {\boldsymbol \nabla }\beta \vert } </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (42a)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (42a)
Line 793: Line 793:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
| style="text-align: center;" | <math>\bar h_c = \max [{l}_j^T {\boldsymbol \nabla }\beta ] {1\over \vert {\boldsymbol \nabla }\beta \vert } \quad ,\quad j=1,2,3</math>
+
| style="text-align: center;" | <math>\bar h_c = \max [{\boldsymbol l}_j^T {\boldsymbol \nabla }\beta ] {1\over \vert {\boldsymbol \nabla }\beta \vert } \quad ,\quad j=1,2,3</math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (42b)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (42b)
Line 1,087: Line 1,087:
 
| style="text-align: center;" | <math>\int _{\Omega }\delta u_i \rho \frac{\partial ^2{u_i}}{\partial{t}^2}\,d\Omega + \int _{\Omega }\delta {\varepsilon }_{ij}{\sigma }_{ij}\,d\Omega{-} \int _{\Omega }\delta u_i b_i \,d\Omega  - \int _{\Gamma _t}\delta {u}_it_i^p \,d\Gamma _t=0</math>
 
| style="text-align: center;" | <math>\int _{\Omega }\delta u_i \rho \frac{\partial ^2{u_i}}{\partial{t}^2}\,d\Omega + \int _{\Omega }\delta {\varepsilon }_{ij}{\sigma }_{ij}\,d\Omega{-} \int _{\Omega }\delta u_i b_i \,d\Omega  - \int _{\Gamma _t}\delta {u}_it_i^p \,d\Gamma _t=0</math>
 
|}
 
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (63a)
+
| style="width: 5px;text-align: right;white-space: nowrap;" | (63a)
 
|}
 
|}
  
Line 1,097: Line 1,097:
 
| style="text-align: center;" | <math>\int _{\Omega }q \left({\Delta p\over K} + {\partial (\Delta u_i) \over \partial x_i} \right)\,d\Omega + \int _\Omega \left[\sum \limits _{i=1}^{n_d}{\partial q \over \partial x_i} \tau _{i}\left({\partial p \over \partial x_i}+\pi _i \right)\right]d\Omega =0</math>
 
| style="text-align: center;" | <math>\int _{\Omega }q \left({\Delta p\over K} + {\partial (\Delta u_i) \over \partial x_i} \right)\,d\Omega + \int _\Omega \left[\sum \limits _{i=1}^{n_d}{\partial q \over \partial x_i} \tau _{i}\left({\partial p \over \partial x_i}+\pi _i \right)\right]d\Omega =0</math>
 
|}
 
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (63b)
+
| style="width: 5px;text-align: right;white-space: nowrap;" | (63b)
 
|}
 
|}
  
Line 1,107: Line 1,107:
 
| style="text-align: center;" | <math>\int _{\Omega }\delta \pi _i \tau _{i} \left({\partial p \over \partial x_i}+\pi _i \right)d\Omega =0\quad \hbox{no sum in }i </math>
 
| style="text-align: center;" | <math>\int _{\Omega }\delta \pi _i \tau _{i} \left({\partial p \over \partial x_i}+\pi _i \right)d\Omega =0\quad \hbox{no sum in }i </math>
 
|}
 
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (63c)
+
| style="width: 5px;text-align: right;white-space: nowrap;" | (63c)
 
|}
 
|}
  
Line 1,119: Line 1,119:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
| style="text-align: center;" | <math>{M} \ddot{\bar {u}} + {g} -  {f}={0}</math>
+
| style="text-align: center;" | <math>{\boldsymbol M} \ddot{\bar {\boldsymbol u}} + {\boldsymbol g} -  {\boldsymbol f}={\boldsymbol 0}</math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (64a)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (64a)
Line 1,129: Line 1,129:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
| style="text-align: center;" | <math>{G}^T \Delta \bar {u} +{C} \Delta \bar {p}+ {L} \bar {p} + {Q}\bar {\boldsymbol \pi }=  {0}</math>
+
| style="text-align: center;" | <math>{\boldsymbol G}^T \Delta \bar {\boldsymbol u} +{\boldsymbol C} \Delta \bar {\boldsymbol p}+ {\boldsymbol L} \bar {\boldsymbol p} + {\boldsymbol Q}\bar {\boldsymbol \pi }=  {\boldsymbol 0}</math>
 
|}
 
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (64b)
+
| style="width: 5px;text-align: right;white-space: nowrap;" | (64b)
 
|}
 
|}
  
Line 1,139: Line 1,139:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
| style="text-align: center;" | <math>{Q}^T \bar {p}+ \bar {C} \bar {\boldsymbol \pi }=  {0}</math>
+
| style="text-align: center;" | <math>{\boldsymbol Q}^T \bar {\boldsymbol p}+ \bar {\boldsymbol C} \bar {\boldsymbol \pi }=  {\boldsymbol 0}</math>
 
|}
 
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (64c)
+
| style="width: 5px;text-align: right;white-space: nowrap;" | (64c)
 
|}
 
|}
  
where <math display="inline">\ddot{\bar {u}}</math> is the nodal acceleration vector,
+
where <math display="inline">\ddot{\bar {\boldsymbol u}}</math> is the nodal acceleration vector,
  
 
{| class="formulaSCP" style="width: 100%; text-align: left;"  
 
{| class="formulaSCP" style="width: 100%; text-align: left;"  
Line 1,163: Line 1,163:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
| style="text-align: center;" | <math>{g}= \int _\Omega {B}^T {\boldsymbol \sigma } d\Omega  </math>
+
| style="text-align: center;" | <math>{\boldsymbol g}= \int _\Omega {\boldsymbol B}^T {\boldsymbol \sigma } d\Omega  </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (66)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (66)
 
|}
 
|}
  
is the internal nodal force vector and the rest of matrices and vectors are defined in the Appendix. Note that the expression of <math display="inline">{g}</math> of eq.(66) is adequate for non linear  analysis.
+
is the internal nodal force vector and the rest of matrices and vectors are defined in the Appendix. Note that the expression of <math display="inline">{\boldsymbol g}</math> of eq.(66) is adequate for non linear  analysis.
  
 
A four steps semi-explicit time integration algorithm can be derived as follows
 
A four steps semi-explicit time integration algorithm can be derived as follows
Line 1,175: Line 1,175:
  
 
'''Step step.'''step
 
'''Step step.'''step
* Compute the nodal velocities <math display="inline">\dot{\bar {u}}^{n+1/2}</math> from Eq.(64a)
+
* Compute the nodal velocities <math display="inline">\dot{\bar {\boldsymbol u}}^{n+1/2}</math> from Eq.(64a)
  
 
{| class="formulaSCP" style="width: 100%; text-align: left;"  
 
{| class="formulaSCP" style="width: 100%; text-align: left;"  
Line 1,182: Line 1,182:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
| style="text-align: center;" | <math> \dot{\bar {u}}^{n+1/2} = \dot{\bar {u}}^{n-1/2}+\Delta t {M}^{-1}_d ({f}^n - {g}^n)</math>
+
| style="text-align: center;" | <math> \dot{\bar {\boldsymbol u}}^{n+1/2} = \dot{\bar {\boldsymbol u}}^{n-1/2}+\Delta t {\boldsymbol M}^{-1}_d ({\boldsymbol f}^n - {\boldsymbol g}^n)</math>
 
|}
 
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (67a)
+
| style="width: 5px;text-align: right;white-space: nowrap;" | (67a)
 
|}
 
|}
  
* Compute the nodal displacements <math display="inline">\bar {u}^{n+1}</math>
+
* Compute the nodal displacements <math display="inline">\bar {\boldsymbol u}^{n+1}</math>
  
 
{| class="formulaSCP" style="width: 100%; text-align: left;"  
 
{| class="formulaSCP" style="width: 100%; text-align: left;"  
Line 1,194: Line 1,194:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
| style="text-align: center;" | <math>\bar {u}^{n+1} = \bar{u}^n + \Delta t \dot{\bar {u}}^{n+1/2} </math>
+
| style="text-align: center;" | <math>\bar {\boldsymbol u}^{n+1} = \bar{\boldsymbol u}^n + \Delta t \dot{\bar {\boldsymbol u}}^{n+1/2} </math>
 
|}
 
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (67b)
+
| style="width: 5px;text-align: right;white-space: nowrap;" | (67b)
 
|}
 
|}
  
* Compute the nodal pressures <math display="inline">\bar {p}^{n+1}</math> from Eq.(64b)
+
* Compute the nodal pressures <math display="inline">\bar {\boldsymbol p}^{n+1}</math> from Eq.(64b)
  
 
{| class="formulaSCP" style="width: 100%; text-align: left;"  
 
{| class="formulaSCP" style="width: 100%; text-align: left;"  
Line 1,206: Line 1,206:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
| style="text-align: center;" | <math>\bar {p}^{n+1}=-[ {C} + {L}]^{-1} [\Delta t {G}^T \dot{\bar {u}}^{n+1/2}- {C} \bar {p}^n + {Q} \bar {\boldsymbol \pi }^n]</math>
+
| style="text-align: center;" | <math>\bar {\boldsymbol p}^{n+1}=-[ {\boldsymbol C} + {\boldsymbol L}]^{-1} [\Delta t {\boldsymbol G}^T \dot{\bar {\boldsymbol u}}^{n+1/2}- {\boldsymbol C} \bar {\boldsymbol p}^n + {\boldsymbol Q} \bar {\boldsymbol \pi }^n]</math>
 
|}
 
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (67c)
+
| style="width: 5px;text-align: right;white-space: nowrap;" | (67c)
 
|}
 
|}
  
Line 1,218: Line 1,218:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
| style="text-align: center;" | <math>\bar {\boldsymbol \pi }^{n+1}=- \bar {C}_d^{-1} {Q}^T \bar {p}^{n+1}</math>
+
| style="text-align: center;" | <math>\bar {\boldsymbol \pi }^{n+1}=- \bar {\boldsymbol C}_d^{-1} {\boldsymbol Q}^T \bar {\boldsymbol p}^{n+1}</math>
 
|}
 
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (67d)
+
| style="width: 5px;text-align: right;white-space: nowrap;" | (67d)
 
|}
 
|}
  
Line 1,230: Line 1,230:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
| style="text-align: center;" | <math>{g}^n =\int _{\Omega ^e} [{B}^T {\boldsymbol \sigma }]^n\,d\Omega  </math>
+
| style="text-align: center;" | <math>{\boldsymbol g}^n =\int _{\Omega ^e} [{\boldsymbol B}^T {\boldsymbol \sigma }]^n\,d\Omega  </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (68)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (68)
Line 1,237: Line 1,237:
 
where the stresses <math display="inline">{\boldsymbol \sigma }^n</math> are obtained by consistent integration of the adequate (non linear) constitutive law.
 
where the stresses <math display="inline">{\boldsymbol \sigma }^n</math> are obtained by consistent integration of the adequate (non linear) constitutive law.
  
Note that steps 1, 2 and 4 are ''fully explicit'' as a diagonal form of matrices <math display="inline">{M}</math> and <math display="inline">\bar{C}</math> has been chosen. The solution of step 3 requires invariably the inverse of a Laplacian matrix. This can be an inexpensive process using an iterative equation solution method (e.g. a preconditioned conjugate gradient method).
+
Note that steps 1, 2 and 4 are ''fully explicit'' as a diagonal form of matrices <math display="inline">{\boldsymbol M}</math> and <math display="inline">\bar{\boldsymbol C}</math> has been chosen. The solution of step 3 requires invariably the inverse of a Laplacian matrix. This can be an inexpensive process using an iterative equation solution method (e.g. a preconditioned conjugate gradient method).
  
For the full incompressible case <math display="inline">K=\infty </math> and <math display="inline">{C}=0 </math> in all above equations.
+
For the full incompressible case <math display="inline">K=\infty </math> and <math display="inline">{\boldsymbol C}=0 </math> in all above equations.
  
 
The critical time step <math display="inline">\Delta t</math> is taken as that of the standard explicit dynamic scheme (see Section 6.5 and [9,37]).
 
The critical time step <math display="inline">\Delta t</math> is taken as that of the standard explicit dynamic scheme (see Section 6.5 and [9,37]).
Line 1,245: Line 1,245:
 
===Fully explicit algorithm===
 
===Fully explicit algorithm===
  
A fully explicit four steps algorithm can be obtained by computing <math display="inline">\bar {p}^{n+1}</math> from step 3 in eq.(67c) as follows
+
A fully explicit four steps algorithm can be obtained by computing <math display="inline">\bar {\boldsymbol p}^{n+1}</math> from step 3 in eq.(67c) as follows
  
 
{| class="formulaSCP" style="width: 100%; text-align: left;"  
 
{| class="formulaSCP" style="width: 100%; text-align: left;"  
Line 1,252: Line 1,252:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
| style="text-align: center;" | <math>\bar{p}^{n+1}= - {C}_d^{-1}[\Delta t {G}^T \dot{\bar {u}}^{n+1/2} - ({C}_d - {L}) \bar{p}^n + {Q} \bar{\boldsymbol \pi }^n] </math>
+
| style="text-align: center;" | <math>\bar{\boldsymbol p}^{n+1}= - {\boldsymbol C}_d^{-1}[\Delta t {\boldsymbol G}^T \dot{\bar {\boldsymbol u}}^{n+1/2} - ({\boldsymbol C}_d - {\boldsymbol L}) \bar{\boldsymbol p}^n + {\boldsymbol Q} \bar{\boldsymbol \pi }^n] </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (69)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (69)
 
|}
 
|}
  
Note that the explicit algorithm is not applicable in the full incompressible limit as the solution of Eq.(69) breaks down for <math display="inline">K=\infty </math> and <math display="inline">{C}=0</math>. The explicit form can however be used with success in problems where quasi-incompressible regions exist adjacent to standard “compressible” zones. Examples of this kind are shown in the next section. In both cases  the semi-explicit and fully explicit schemes gave identical results with important savings in both computer time and memory storage requirements obtained when using the explicit form.
+
Note that the explicit algorithm is not applicable in the full incompressible limit as the solution of Eq.(69) breaks down for <math display="inline">K=\infty </math> and <math display="inline">{\boldsymbol C}=0</math>. The explicit form can however be used with success in problems where quasi-incompressible regions exist adjacent to standard “compressible” zones. Examples of this kind are shown in the next section. In both cases  the semi-explicit and fully explicit schemes gave identical results with important savings in both computer time and memory storage requirements obtained when using the explicit form.
  
 
===6.4 Thermal coupled effects===
 
===6.4 Thermal coupled effects===
Line 1,352: Line 1,352:
 
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
 
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
 
|-
 
|-
|[[Image:Draft_Samper_805805629-fig1cast.png|351px|Finite element discretization of the aluminium casting.]]
+
|[[Image:Draft_Samper_805805629-test-fig1cast.png|351px|Finite element discretization of the aluminium casting.]]
 
|- style="text-align: center; font-size: 75%;"
 
|- style="text-align: center; font-size: 75%;"
 
| colspan="1" | '''Figure 1:''' Finite element discretization of the aluminium casting.
 
| colspan="1" | '''Figure 1:''' Finite element discretization of the aluminium casting.
Line 1,360: Line 1,360:
 
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
 
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
 
|-
 
|-
|[[Image:Draft_Samper_805805629-fig2cast.png|351px|Temperature die-cycling.]]
+
|[[Image:Draft_Samper_805805629-test-fig2cast.png|351px|Temperature die-cycling.]]
 
|- style="text-align: center; font-size: 75%;"
 
|- style="text-align: center; font-size: 75%;"
 
| colspan="1" | '''Figure 2:''' Temperature die-cycling.
 
| colspan="1" | '''Figure 2:''' Temperature die-cycling.
Line 1,372: Line 1,372:
 
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
 
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
 
|-
 
|-
|[[Image:Draft_Samper_805805629-MATER.png|300px|]]
+
|[[Image:Draft_Samper_805805629-test-MATER.png|300px|]]
|[[Image:Draft_Samper_805805629-MATER-01.png|300px|]]
+
|[[Image:Draft_Samper_805805629-test-MATER-01.png|300px|]]
 
|-
 
|-
|[[Image:Draft_Samper_805805629-MATER-02.png|300px|]]
+
|[[Image:Draft_Samper_805805629-test-MATER-02.png|300px|]]
|[[Image:Draft_Samper_805805629-MATER-03.png|300px|]]
+
|[[Image:Draft_Samper_805805629-test-MATER-03.png|300px|]]
 
|-
 
|-
| colspan="2"|[[Image:Draft_Samper_805805629-MATER-04.png|300px|Filling evolution: pressure die-casting simulation.]]
+
| colspan="2"|[[Image:Draft_Samper_805805629-test-MATER-04.png|300px|Filling evolution: pressure die-casting simulation.]]
 
|- style="text-align: center; font-size: 75%;"
 
|- style="text-align: center; font-size: 75%;"
 
| colspan="2" | '''Figure 3:''' Filling evolution: pressure die-casting simulation.
 
| colspan="2" | '''Figure 3:''' Filling evolution: pressure die-casting simulation.
Line 1,386: Line 1,386:
 
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
 
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
 
|-
 
|-
|[[Image:Draft_Samper_805805629-tempe_1.png|210px|]]
+
|[[Image:Draft_Samper_805805629-test-tempe_1.png|210px|]]
|[[Image:Draft_Samper_805805629-tempe_2.png|210px|]]
+
|[[Image:Draft_Samper_805805629-test-tempe_2.png|210px|]]
|[[Image:Draft_Samper_805805629-tempe_3.png|210px|Temperature evolution during cooling phase.]]
+
|[[Image:Draft_Samper_805805629-test-tempe_3.png|210px|Temperature evolution during cooling phase.]]
 
|- style="text-align: center; font-size: 75%;"
 
|- style="text-align: center; font-size: 75%;"
 
| colspan="3" | '''Figure 4:''' Temperature evolution during cooling phase.
 
| colspan="3" | '''Figure 4:''' Temperature evolution during cooling phase.
Line 1,396: Line 1,396:
 
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
 
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
 
|-
 
|-
|[[Image:Draft_Samper_805805629-solid_frac.png|360px|]]
+
|[[Image:Draft_Samper_805805629-test-solid_frac.png|360px|]]
|[[Image:Draft_Samper_805805629-solid_frac-01.png|360px|]]
+
|[[Image:Draft_Samper_805805629-test-solid_frac-01.png|360px|]]
 
|-
 
|-
|[[Image:Draft_Samper_805805629-solid_frac-02.png|360px|]]
+
|[[Image:Draft_Samper_805805629-test-solid_frac-02.png|360px|]]
|[[Image:Draft_Samper_805805629-solid_frac-03.png|360px|Liquid-fraction evolution during phase-change.]]
+
|[[Image:Draft_Samper_805805629-test-solid_frac-03.png|360px|Liquid-fraction evolution during phase-change.]]
 
|- style="text-align: center; font-size: 75%;"
 
|- style="text-align: center; font-size: 75%;"
 
| colspan="2" | '''Figure 5:''' Liquid-fraction evolution during phase-change.
 
| colspan="2" | '''Figure 5:''' Liquid-fraction evolution during phase-change.
Line 1,408: Line 1,408:
 
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
 
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
 
|-
 
|-
|[[Image:Draft_Samper_805805629-fig6cast.png|351px|Aluminium casting. Stress-trace and von Mises deviatoric stress indicator during  phase-change (plane xy).]]
+
|[[Image:Draft_Samper_805805629-test-fig6cast.png|351px|Aluminium casting. Stress-trace and von Mises deviatoric stress indicator during  phase-change (plane xy).]]
 
|- style="text-align: center; font-size: 75%;"
 
|- style="text-align: center; font-size: 75%;"
 
| colspan="1" | '''Figure 6:''' Aluminium casting. Stress-trace and von Mises deviatoric stress indicator during  phase-change (plane xy).
 
| colspan="1" | '''Figure 6:''' Aluminium casting. Stress-trace and von Mises deviatoric stress indicator during  phase-change (plane xy).
Line 1,415: Line 1,415:
 
===7.2 Side pressing of a cylinder===
 
===7.2 Side pressing of a cylinder===
  
A cylinder 100 mm long with a radius of 100 mm is subjected to sidepressing between two plane dies. It is compressed to 100 mm. The material properties are the following: <math display="inline">E =217</math> GPa, <math display="inline">\nu =0.3</math>, <math display="inline">\rho =7830</math> kg/m<math display="inline">^3</math>, <math display="inline">\sigma _0=170</math> MPa, <math display="inline">H=30</math> MPa,  friction coefficient = 0.2. The die velocity is assumed to be 2 m/s. Initial set-up is shown in Figure [[#img-7a|7a]]. A quarter of a cylinder was discretized with tri-linear hexahedra and linear tetrahedra meshes.
+
A cylinder 100 mm long with a radius of 100 mm is subjected to sidepressing between two plane dies. It is compressed to 100 mm. The material properties are the following: <math display="inline">E =217</math> GPa, <math display="inline">\nu =0.3</math>, <math display="inline">\rho =7830</math> kg/m<math display="inline">^3</math>, <math display="inline">\sigma _0=170</math> MPa, <math display="inline">H=30</math> MPa,  friction coefficient = 0.2. The die velocity is assumed to be 2 m/s. Initial set-up is shown in Figure [[#img-7|7]]. A quarter of a cylinder was discretized with tri-linear hexahedra and linear tetrahedra meshes.
  
Figure [[#img-7|7]] shows the results obtained using the hexahedral mesh and a standard mixed formulation. The results show the distribution of the effective plastic strain and pressure on the deformed shape. The sensitivity of the FIC results with the expression of the intrinsic time parameter of Eq.(71) was studied by defining <math display="inline">\tau ^{(e)} =a  {[h^{(e)}]^2\over G}</math> and solving the problem for different values of <math display="inline">a</math>. The results obtained with the FIC method are shown in Figs. cylt-fine-fast-fic01 and cylt-fine-fast-fic for <math display="inline">a =0.1</math> and <math display="inline">0.03</math>, respectively. The alternative expression for <math display="inline">\tau ^{(e)}</math> of Eq.(74) has also been studied. The results for this case are shown in Fig. cylt-fine-fast-fic-dtcr. Quite a good agreement can be seen between the FIC solutions and the reference solution with the best results for the effective plastic strain obtained for <math display="inline">\tau ^{(e)}</math> calculated according to Eq. (71) with <math display="inline">a=0.03</math>. The results for the two alternative formulae for <math display="inline">\tau ^{(e)}</math> are similar, but those obtained using Eq. (71) seem to be slightly better. In any case, the results on the pressure distribution are quite insentitive to the value of <math display="inline">\tau ^{(e)}</math>. This is also confirmed in Fig. cylt-liniab, which displays the distribution of the pressure along the line <math display="inline">ABCDEA</math> defined in Fig. cylt-liniaa. A small perturbation can be seen at the sharp edges of the deformed body.
+
Figure [[#img-8|8]] shows the results obtained using the hexahedral mesh and a standard mixed formulation. The results show the distribution of the effective plastic strain and pressure on the deformed shape. The sensitivity of the FIC results with the expression of the intrinsic time parameter of Eq.(71) was studied by defining <math display="inline">\tau ^{(e)} =a  {[h^{(e)}]^2\over G}</math> and solving the problem for different values of <math display="inline">a</math>. The results obtained with the FIC method are shown in Figs. [[#img-9|9]] and [[#img-10|10]] for <math display="inline">a =0.1</math> and <math display="inline">0.03</math>, respectively. The alternative expression for <math display="inline">\tau ^{(e)}</math> of Eq.(74) has also been studied. The results for this case are shown in Fig. [[#img-11|11]]. Quite a good agreement can be seen between the FIC solutions and the reference solution with the best results for the effective plastic strain obtained for <math display="inline">\tau ^{(e)}</math> calculated according to Eq. (71) with <math display="inline">a=0.03</math>. The results for the two alternative formulae for <math display="inline">\tau ^{(e)}</math> are similar, but those obtained using Eq. (71) seem to be slightly better. In any case, the results on the pressure distribution are quite insentitive to the value of <math display="inline">\tau ^{(e)}</math>. This is also confirmed in Fig. [[#img-12|12]]b, which displays the distribution of the pressure along the line <math display="inline">ABCDEA</math> defined in Fig. [[#img-12|12]]a. A small perturbation can be seen at the sharp edges of the deformed body.
  
 
All the calculations here have been carried out using fully explicit version of the algorithm which is more efficient than the semi-implicit one with giving practically the same results.
 
All the calculations here have been carried out using fully explicit version of the algorithm which is more efficient than the semi-implicit one with giving practically the same results.
  
<div id='img-7a'></div>
 
 
<div id='img-7'></div>
 
<div id='img-7'></div>
 
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
 
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
 
|-
 
|-
|[[Image:Draft_Samper_805805629-cyltn3-pi-v5e5-msh.png|270px|]]
+
|[[Image:Draft_Samper_805805629-test-cyltn3-pi-v5e5-msh.png|270px|]]
|[[Image:Draft_Samper_805805629-cylhc-msh-new.png|270px|Sidepressing of a cylinder: (a) initial tetrahedral mesh; (b) initial hexahedral mesh]]
+
|[[Image:Draft_Samper_805805629-test-cylhc-msh-new.png|270px|Sidepressing of a cylinder: (a) initial tetrahedral mesh; (b) initial hexahedral mesh]]
| (a) Sidepressing of a cylinder: (a) initial tetrahedral mesh; (b) initial hexahedral mesh
+
|- style="text-align: center; font-size: 75%;"
 +
| colspan="2" | '''Figure 7:''' Sidepressing of a cylinder: (a) initial tetrahedral mesh; (b) initial hexahedral mesh
 +
|}
 +
 
 +
<div id='img-8'></div>
 +
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
 +
|-
 +
|[[Image:Draft_Samper_805805629-test-cylhc-ep-new.png|228px|]]
 +
|[[Image:Draft_Samper_805805629-test-cylhc-hp-new.png|228px|Sidepressing of a cylinder, mixed formulation,  hexahedral mesh: (a) effective plastic strain; (b) pressure distribution ]]
 +
|- style="text-align: center; font-size: 75%;"
 +
| colspan="2" | '''Figure 8:''' Sidepressing of a cylinder, mixed formulation,  hexahedral mesh: (a) effective plastic strain; (b) pressure distribution
 +
|}
 +
 
 +
<div id='img-9'></div>
 +
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
 +
|-
 +
|[[Image:Draft_Samper_805805629-test-cylt-fine-fast-fic01-eps.png|228px|]]
 +
|[[Image:Draft_Samper_805805629-test-cylt-fine-fast-fic01-p.png|228px|Sidepressing of a cylinder, FIC algorithm  (α=0.1),  tetrahedra, mesh  of 22186 elements: (a) effective plastic strain; (b) pressure distribution ]]
 +
|- style="text-align: center; font-size: 75%;"
 +
| colspan="2" | '''Figure 9:''' Sidepressing of a cylinder, FIC algorithm  (<math>\alpha=0.1</math>),  tetrahedra, mesh  of 22186 elements: (a) effective plastic strain; (b) pressure distribution
 +
|}
 +
 
 +
<div id='img-10'></div>
 +
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
 +
|-
 +
|[[Image:Draft_Samper_805805629-test-cylt-fine-fast-fic-eps.png|228px|]]
 +
|[[Image:Draft_Samper_805805629-test-cylt-fine-fast-fic-p.png|228px|Sidepressing of a cylinder, FIC algorithm  (a=0.03),  tetrahedra, mesh  of 22186 elements: (a) effective plastic strain; (b) pressure distribution ]]
 +
|- style="text-align: center; font-size: 75%;"
 +
| colspan="2" | '''Figure 10:''' Sidepressing of a cylinder, FIC algorithm  (<math>a=0.03</math>),  tetrahedra, mesh  of 22186 elements: (a) effective plastic strain; (b) pressure distribution
 +
|}
 +
 
 +
<div id='img-11'></div>
 +
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
 +
|-
 +
|[[Image:Draft_Samper_805805629-test-cylt-fine07-fast-v1e12-ab2.png|228px|]]
 +
|[[Image:Draft_Samper_805805629-test-cylt-fine07-fast-v1e12-p-ab2.png|228px|Sidepressing of a cylinder, FIC algorithm  (τ<sup>(e)</sup> calculated according to Eq. (74)),  tetrahedra, mesh  of 22186 elements: (a) effective plastic strain; (b) pressure distribution ]]
 +
|- style="text-align: center; font-size: 75%;"
 +
| colspan="2" | '''Figure 11:''' Sidepressing of a cylinder, FIC algorithm  (<math>\tau ^{(e)}</math> calculated according to Eq. (74)),  tetrahedra, mesh  of 22186 elements: (a) effective plastic strain; (b) pressure distribution
 +
|}
 +
 
 +
<div id='img-12'></div>
 +
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
 +
|-
 +
|[[Image:Draft_Samper_805805629-test-cylt-linia.png|300px|]]
 +
|[[Image:Draft_Samper_805805629-test-cylt-p-mo-fic.png|350px|a) Definition of the  line for comparison of pressure distribution b) Pressure distribution along the line ABCDEA]]
 +
|- style="text-align: center; font-size: 75%;"
 +
| colspan="2" | '''Figure 12:''' a) Definition of the  line for comparison of pressure distribution b) Pressure distribution along the line <math>ABCDEA</math>
 +
|}
 +
 
 +
===7.3 Backward extrusion===
 +
 
 +
Backward extrusion of a cylinder made of steel 16MNCr5 has been analysed using an axi-symmetric formulation. This is a benchmark example of the finite element program for forming simulation MARC/Autoforge [marc]. The tooling and billet geometry are given in Fig. [[#img-13|13]]a. Initial material dimensions are the following: length 30 mm and diameter 30 mm. The punch of diameter 20 mm has a prescribed stroke of 28 mm. Material properties are as follows: Young's modulus <math display="inline">E=3.24\cdot 10^5</math> MPa, Poissson's coefficient <math display="inline">\nu=0.3</math>, material density <math display="inline">\rho=8120</math> kg/m<math display="inline">^3</math>, yield stress <math display="inline">\sigma _{Y0}=300</math> MPa and hardening modulus <math display="inline">H=50</math> MPa. Friction between the material and tools is defined by the Coulomb friction coefficient <math display="inline">\mu=0.1</math>.
 +
 
 +
<div id='img-13'></div>
 +
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
 +
|-
 +
|[[Image:Draft_Samper_805805629-test-bext-geo.png|330px|]]
 +
|[[Image:Draft_Samper_805805629-test-bext4a-quad-nurbs-sprn9-eps.png|266px|]]
 +
|-
 +
| colspan="2"|[[Image:Draft_Samper_805805629-test-bext4a_036-eps.png|279px|Backward extrusion a) geometry definition. Final deformed shape with effective plastic strain distribution; b)  solution with quadrilaterals and mixed formulation; c)  solution with triangles and the FIC algorithm ]]
 +
|- style="text-align: center; font-size: 75%;"
 +
| colspan="2" | '''Figure 13:''' Backward extrusion a) geometry definition. Final deformed shape with effective plastic strain distribution; b)  solution with quadrilaterals and mixed formulation; c)  solution with triangles and the FIC algorithm
 +
|}
 +
 
 +
The simulation of the backward extrusion process was carried out with a particularization of the FIC formulation of Section 6 for axisymmetric solids. Regeneration of the  mesh was performed when element distorsion was excessive. Figures [[#img-13|13]]b and [[#img-13|13]]c show the results in the form of the final deformed shape with the distribution of the effective plastic strain obtained using quadrilaterals with a mixed formulation, and using triangles and the FIC algorithm, respectively. The results are in a good agreement with the solution given in [marc]. This example demonstrates the efficiency of the FIC algorithm for simulation of bulk forming processes. Different stages of the forming process are shown in Fig. [[#img-14|14]].
 +
 
 +
<div id='img-14'></div>
 +
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
 +
|-
 +
|[[Image:Draft_Samper_805805629-test-bext4a_002-eps.png|330px|]]
 +
|[[Image:Draft_Samper_805805629-test-bext4a_014-eps.png|330px|]]
 +
|-
 +
| colspan="2"|[[Image:Draft_Samper_805805629-test-bext4a_022-eps.png|330px|Backward extrusion &#8211; deformed shapes with effective plastic strain distribution at different stages of forming. Solution with triangles and the FIC algorithm ]]
 +
|- style="text-align: center; font-size: 75%;"
 +
| colspan="2" | '''Figure 14:''' Backward extrusion &#8211; deformed shapes with effective plastic strain distribution at different stages of forming. Solution with triangles and the FIC algorithm
 +
|}
 +
 
 +
===7.4 Forming of a hose-clamp band===
 +
 
 +
The stabilized formulation was applied to model the manufacturing of a hose-clamp band  of steel AISI 409L (Fig. [[#img-15|15]]). The initial set-up of the tooling and band is shown in Fig. [[#img-16|16]]. A series of grooves are forged in the band by the roll passing over the band placed on the toothed punch. The band thickness is 0.7 mm and its width 8 mm. Plane strain conditions have been assumed. Material properties are as follows: Young's modulus <math display="inline">E=2.1\cdot 10^5</math> MPa, Poissson's coefficient <math display="inline">\nu=0.33</math>, material density <math display="inline">\rho=7800</math> kg/m<math display="inline">^3</math>, the true stress&#8211;true strain relationship is given by the power Ludwik&#8211;Nadai equation <math display="inline">\sigma _{Y}=623(0.36822\cdot 10^{-2}+\varepsilon _p)^{0.1362}</math> MPa, friction between the material and tools is defined by the Coulomb friction coefficient <math display="inline">\mu=0.1</math>.
 +
 
 +
<div id='img-15'></div>
 +
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
 +
|-
 +
|[[Image:Draft_Samper_805805629-test-hose-clip.png|330px|A hose clamp ]]
 +
|- style="text-align: center; font-size: 75%;"
 +
| colspan="1" | '''Figure 15:''' A hose clamp
 +
|}
 +
 
 +
<div id='img-16'></div>
 +
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
 +
|-
 +
|[[Image:Draft_Samper_805805629-test-p1222e_ini.png|210px|A hose clamp &#8211; initial set-up ]]
 +
|- style="text-align: center; font-size: 75%;"
 +
| colspan="1" | '''Figure 16:''' A hose clamp &#8211; initial set-up
 +
|}
 +
 
 +
<div id='img-17'></div>
 +
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
 +
|-
 +
|[[Image:Draft_Samper_805805629-test-p1222d_12-eps.png|270px|]]
 +
|[[Image:Draft_Samper_805805629-test-p1222e_29-eps.png|240px|A hose clamp &#8211; deformed shapes at different stages of forming with distribution of effective plastic strain ]]
 +
|- style="text-align: center; font-size: 75%;"
 +
| colspan="2" | '''Figure 17:''' A hose clamp &#8211; deformed shapes at different stages of forming with distribution of effective plastic strain
 +
|}
 +
 
 +
Simulation was carried out using linear triangular elements and  the FIC formulation. As in the previous example the  meshes where  regenerated when element distorsion was excessive. The purpose of the simulation was to check if the expected groove depth and tooth height in the band were obtained. Figures [[#img-17|17]]a and [[#img-17|17]]b show the results in the form of the  deformed shape with the distribution of the effective plastic strain at different stages of forming. The results are in good agreement with the real process. Figure  [[#img-18|18]] shows a detail  of the deformed shape with finite element discretization and  distribution of effective plastic strain. In this figure the obtained dimensions are compared with the required ones shown in brackets. Effects of elastic springback are also clearly seen.
 +
 
 +
<div id='img-18'></div>
 +
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
 +
|-
 +
|[[Image:Draft_Samper_805805629-test-clamp-detail.png|360px|Detail of the deformed shape with finite element discretization and  distribution of effective plastic strain ]]
 +
|- style="text-align: center; font-size: 75%;"
 +
| colspan="1" | '''Figure 18:''' Detail of the deformed shape with finite element discretization and  distribution of effective plastic strain
 +
|}
 +
 
 +
===Remark===
 +
 
 +
The examples presented show that the FIC/FEM formulation is an effective procedure for solving bulk metal forming problems involving full or quasi-incompressible situations. The key advantage of the FIC approach versus more standard mixed FEM formulations is that it provides a natural theoretical framework for equal order finite element interpolations for the velocity and pressure variables, both in the context of implicit and explicit solution schemes. We note the simplicity and effectiveness of the full explicit algorithm as demonstrated in the examples presented. The FIC formulation reproduces also the best feature of the so called stabilized FEM method for incompressible problems, such as the CBS scheme [20,21,25,26,58], the pressure gradient operator method [22] and the subgrid scale method [23] among others.
 +
 
 +
==8 LAGRANGIAN FLOWS. THE PARTICLE FINITE ELEMENT  METHOD==
 +
 
 +
===8.1 The Particle Finite Element Method (PFEM)===
 +
 
 +
The  Lagrangian formulation is  an excellent procedure for treating bulk forming processes involving the interaction of fluids and solids using a unified formulation. An important advantage of the Lagrangian formulation is that both the  motions of the solid  and the fluid  are defined in the same frame of reference and modelled with the some governing equations.
 +
 
 +
The Lagrangian fluid flow equations can be simply obtained by noting that the velocity of the mesh nodes and that of the fluid particles are the same. Hence the the convective terms vanish in the momentum equations, while the rest of the fluid flow equations remain unchanged. The resulting governing equations have an identical form as those of the Stokes flow problem, with the motion of the flow particles being referred now to a Lagrangian coordinate frame.
 +
 
 +
The FEM algorithms for solving the Lagrangian flow equations are very similar to those presented for incompressible solids in a previous section. Here we present a particular class of Lagrangian formulation  called the ''particle finite element method'' (PFEM) [38&#8211;41]. The PFEM treats the mesh nodes in the fluid and solid domains as  dimensionless particles which can freely move an even separate from the main fluid domain representing, for instance, the effect of fluid  drops. A finite element mesh connects the nodes defining the discretized domain where the governing equations are solved in the standard FEM fashion.
 +
 
 +
The quality of the numerical solution  depends on the discretization chosen as in the standard FEM. Adaptive mesh refinement techniques can be used to improve the solution in zones where large motions of the fluid or the structure occur.
 +
 
 +
A typical solution with the PFEM involves the following steps.
 +
 
 +
<ol>
 +
 
 +
<li>Identify the external boundaries for both the fluid and solid domains. This is  an essential step as some boundaries (such as the free surface in the fluid) may be severely distorted during the solution process including separation and re-entring of nodes. The Alpha Shape method  is used for the boundary definition (see Section 8.4).  </li>
 +
 
 +
<li>Discretize the fluid and solid domains with a finite element mesh. For the mesh generation process we use and extended Delaunay technique  [52].  </li>
 +
 
 +
<li>Solve iteratively the coupled Lagrangian equations of motion for  the fluid and the solid domains. Compute the relevant state variables in both domains at each time step: velocities, pressure and viscous stresses in the fluid and displacements, stresses and strains in the solid.  </li>
 +
 
 +
The iterative solution scheme chosen is a particularization of the fractional step algorithm of Section 3.1 for <math display="inline">\alpha =1</math>. In summary the solution steps are the following.
 +
 
 +
''Step 3.1.'' Compute the predicted velocities (viz, Eq.(31))
 +
 
 +
{| class="formulaSCP" style="width: 100%; text-align: left;"
 +
|-
 +
|
 +
{| style="text-align: left; margin:auto;width: 100%;"
 +
|-
 +
| style="text-align: center;" | <math>
 +
 
 +
\tilde{\boldsymbol v}^{n+1,i+1} = \bar{\boldsymbol v}^n - \Delta t {\boldsymbol M}^{-1}_d [{\boldsymbol K} \bar{\boldsymbol v}^{n} - {\boldsymbol G}\bar {\boldsymbol p}^{n} -{\boldsymbol f}^{n+1}] </math>
 +
|}
 +
| style="width: 5px;text-align: right;white-space: nowrap;" | (75)
 +
|}
 +
 
 +
''Step 3.2.'' Compute <math display="inline">\bar {\boldsymbol p}^{n+1,i+1}</math> from Eq.(32) for <math display="inline">\alpha =1</math>.
 +
 
 +
''Step 3.3.'' Compute explicitely <math display="inline"> \bar{\boldsymbol v}^{n+1,i+1}</math> from Eq.(33) with <math display="inline">\alpha =1</math>.
 +
 
 +
''Step 3.4.'' Compute <math display="inline"> \bar{\boldsymbol \pi }^{n+1,i+1}</math> explicitely from Eq.(35).
 +
 
 +
''Step 3.5.'' Solve for the motion of the solid. This can be performed by integrating the dynamic equations of motion in the solid. Here the algorithm of Section 6 can be used, as it applies to both “compressible” and incompressible materials.
 +
 
 +
''Step 3.6.'' Estimate a new position of the mesh nodes in terms of the time increment size as
 +
 
 +
{| class="formulaSCP" style="width: 100%; text-align: left;"
 +
|-
 +
|
 +
{| style="text-align: left; margin:auto;width: 100%;"
 +
|-
 +
| style="text-align: center;" | <math>
 +
 
 +
{\boldsymbol x}_j^{n+1,i+1} = {\boldsymbol x}_j^{n}+\bar {\boldsymbol u}_j^{n+1,i+1} \Delta t </math>
 +
|}
 +
| style="width: 5px;text-align: right;white-space: nowrap;" | (76)
 +
|}
 +
 
 +
It is important to note that all matrices in Steps 3.1&#8211;3.5 are evaluated at the last predicted configuration <math display="inline">\Omega ^{n+1,i}</math>.
 +
 
 +
In steps 3.1&#8211;3.6 superindex <math display="inline">i</math> denotes the iteration within each time step.
 +
 
 +
''Step 3.7.'' Check the convergence of the velocity and pressure fields in the fluid and the displacements, strains and stresses in the solid. If convergence is achieved froze the final position of the mesh nodes and move to the next time increment (step 4), otherwise return to step 3.1 for the next iteration.
 +
 
 +
<li>Go back to step 1 and repeat the solution process for the next time step. </li>
 +
 
 +
</ol>
 +
 
 +
Above algorithm can be found to be analogous to the standard ''updated lagrangian'' scheme typically used in non linear solid mechanics problems [8,57].
 +
 
 +
Despite the motion of the nodes within the iterative process, in general there is no need to regenerate the mesh at each iteration.  ''In the examples presented in the paper the mesh in the fluid domain has been regenerated at each time step''. A cheaper alternative is to generate a new mesh only after a prescribed number of converged time steps, or when  the nodal displacements induce significant geometrical distorsions in some elements.
 +
 
 +
The boundary conditions are applied as described in Section 3.1.
 +
 
 +
In the examples presented in the paper the time increment size has been chosen as
 +
 
 +
{| class="formulaSCP" style="width: 100%; text-align: left;"
 +
|-
 +
|
 +
{| style="text-align: left; margin:auto;width: 100%;"
 +
|-
 +
| style="text-align: center;" | <math>\Delta t =\min (\Delta t_i ) \quad \hbox{with}\quad \Delta t_i ={\vert {\boldsymbol v}\vert \over h_i^{\min }} </math>
 +
|}
 +
| style="width: 5px;text-align: right;white-space: nowrap;" | (77)
 +
|}
 +
 
 +
where <math display="inline">h_i^{\min }</math> is the distance between node <math display="inline">i</math> and the closest node in the mesh.
 +
 
 +
===8.2 Treatment of contact between fluid and solid interfaces===
 +
 
 +
The  condition of prescribed velocities or pressures at the solid boundaries in the PFEM are  applied in strong form to the boundary nodes. These nodes might belong to fixed external boundaries or to moving boundaries linked to the interacting solids. In some problems it is useful to define a layer of nodes adjacent to the external boundary in the fluid where the condition of prescribed velocity is imposed. These nodes typically remain fixed during the solution process. Contact between liquid particles and the solid boundaries is accounted for by the incompressibility condition which ''naturally prevents the liquid nodes to penetrate into the solid boundaries'' [38&#8211;41].
 +
 
 +
===8.3 Generation of a new mesh===
 +
 
 +
One of the key points for the success of the PFEM  is the fast regeneration of a finite element mesh at every time step on the basis of the position of the nodes in the space domain. In our work the mesh is typically generated using the so called extended Delaunay tesselation (EDT) [38,52]. The EDT allows one to generate non standard meshes combining elements of  arbitrary polyhedrical shapes  (triangles, quadrilaterals and other polygons in 2D and tetrahedra, hexahedra and arbitrary polyhedra in 3D) in a computing time of order <math display="inline">n</math>, where <math display="inline">n</math> is the total number of nodes in the mesh. The <math display="inline">C^\circ </math> continuous shape functions of the elements are obtained using the so called meshless finite element interpolation (MFEM) [53].
 +
 
 +
===8.4 Identification of boundary surfaces===
 +
 
 +
The PFEM requires the correct definition of the boundary domain. Sometimes, boundary nodes are clearly distinguished from internal nodes. In other cases, the total set of nodes is the only information available and the algorithm must recognize the boundary nodes.
 +
 
 +
Considering that the nodes follow a variable <math display="inline">h(x)</math>  distribution, where <math display="inline">h(x)</math> is the minimum distance between two nodes, the following criterion has been used. For each two nodes (three nodes in 3D) a circle (a sphere in 3D) of radius equal to <math display="inline">\alpha h</math> containing the nodes is plotted. ''All nodes laying on a circle (sphere) with a radius  greater than <math>\alpha h</math>, not containing any other node in the interior are considered as boundary nodes''. In practice, <math display="inline">\alpha </math>  is a parameter close to, but greater than one. This criterion coincides with the Alpha Shape concept [39,55].
 +
 
 +
The method  also allows to identify isolated fluid particles (nodes) outside the main fluid domain. These particles are treated as part of the external boundary where the pressure is fixed to the atmospheric value. We note that the “flying particles” are in fact “points” ''with no mass'' which motion is followed by integrating the dynamic equations of continuum mechanics for known values of the mass, the initial velocity and acceleration, gravity body forces and a prescribed zero (atmospheric) pressure.
 +
 
 +
===8.5 Temperature coupled effects===
 +
 
 +
Thermal-mechanical problems can be effectively treated with the PFEM. This requires solving for the temperature field at each time step, accounting for the couplings  induced by the mechanical equations. The effect of temperature in the mechanical problem is introduced via the constitutive equation in the usual manner.
 +
 
 +
The form of the heat transfer equation is identical to that of Eq.(70). Recall  that the  Lagrangian formulation eliminates the convective term in the heat transfer equation and hence the FIC stabilization is here not needed [41].
 +
 
 +
===Remark===
 +
 
 +
The key advantage of the PFEM versus conventional  particle method is that it retains the best feature of the FEM to solve problems in continuum mechanics using a variational approach. A standard finite element mesh is used to discretize in space the problem variables and for solving the governing equations precisely as in the standard FEM. The combination of the lagrangian FEM with the Alpha-Shape approach for identification of the  domain boundary at each time step and the frequent regeneration of the finite element mesh are the distinct features of the PFEM.
 +
 
 +
==9 APPLICATIONS OF THE PFEM TO METAL FORMING PROCESSES==
 +
 
 +
===9.1 Sloshing problem===
 +
 
 +
The sloshing example shown in Figure 19 ilustrates the ability of the PFEM to simulate the flow of liquids within closed cavities with large motions of the free surface. This feature of the PFEM is essential to model the filling of moulds as shown in the next examples. More applications of the  PFEM to sloshing problemas are reported in [39].
 +
 
 +
<div id='img-19'></div>
 +
<div id='img-19'></div>
 +
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
 +
|-
 +
|[[Image:Draft_Samper_805805629-test-Fig5.png|351px|PFEM results for a large amplitude sloshing problem]]
 +
|[[Image:Draft_Samper_805805629-test-Figure19_1.png|390px|]]
 +
|- style="text-align: center; font-size: 75%;"
 +
| (19) PFEM results for a large amplitude sloshing problem
 +
|-
 +
|[[Image:Draft_Samper_805805629-test-Figure19_2.png|390px|]]
 +
|[[Image:Draft_Samper_805805629-test-Figure19_3.png|390px|]]
 +
|-
 +
|[[Image:Draft_Samper_805805629-test-Figure19_4.png|390px|]]
 +
|[[Image:Draft_Samper_805805629-test-Figure19_5.png|390px|]]
 +
|-
 +
| colspan="2"|[[Image:Draft_Samper_805805629-test-Figure19_6.png|390px|a) Filling of a 2D mould with a powder material using the PFEM]]
 +
|- style="text-align: center; font-size: 75%;"
 +
| colspan="2" | '''Figure 19:''' a) Filling of a 2D mould with a powder material using the PFEM
 +
|}
 +
 
 +
A cylinder 100 mm long with a radius of 100 mm is subjected to sidepressing between two plane dies. It is compressed up to 100 mm. The material properties are the following: <math display="inline">E =217</math> GPa, <math display="inline">\nu =0.3</math>, <math display="inline">\rho =7830</math> kg/m<math display="inline">^3</math>, <math display="inline">\sigma _0=170</math> MPa, <math display="inline">H=30</math> MPa, friction coefficient = 0.2. The die velocity is assumed to be 2 m/s. A quarter of a cylinder was discretized using two meshes of 22186 linear tetraedra (4369 nodes) and 1296 hexaedra (1641 nodes) as shown in Figure [[#img-20|20]]. Figure [[#img-21|21]]  shows numerical results obtained with the four node hexaedral mesh and an explicit mixed velocity-pressure formulation. Figure [[#img-22|22]] shows the pressure and effective plastic strain distributions obtained with the mesh of linear tetrahedra and the full explicit algorithm presented. The maximum and minimum time steps for this solution where <math display="inline">0.5 \times 10^{-5}</math> and <math display="inline">0.25 \times 10^{-5}</math>, respectively. Good agreement between the results obtained with the FIC and the mixed formulation is achieved.
 +
 
 +
Very similar results were obtained with the linear tetrahedral mesh using the more expensive semi-explicit scheme.
 +
 
 +
<div id='img-20'></div>
 +
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
 +
|-
 +
|[[Image:Draft_Samper_805805629-test-cyltn3-pi-v5e5-msh.png|330px|]]
 +
|[[Image:Draft_Samper_805805629-test-cylhc-msh-new.png|300px|Sidepressing of a cylinder: (a) initial tetrahedral mesh; (b) initial hexahedral mesh]]
 +
|- style="text-align: center; font-size: 75%;"
 +
| colspan="2" | '''Figure 20:''' Sidepressing of a cylinder: (a) initial tetrahedral mesh; (b) initial hexahedral mesh
 +
|}
 +
 
 +
<div id='img-21'></div>
 +
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
 +
|-
 +
|[[Image:Draft_Samper_805805629-test-cylhc-hp-new.png|330px|]]
 +
|[[Image:Draft_Samper_805805629-test-cylhc-ep-new.png|330px|Sidepressing of a cylinder &#8211; results obtained using mixed formulation: a) pressure distribution, b) effective plastic distribution ]]
 +
|- style="text-align: center; font-size: 75%;"
 +
| colspan="2" | '''Figure 21:''' Sidepressing of a cylinder &#8211; results obtained using mixed formulation: a) pressure distribution, b) effective plastic distribution
 +
|}
 +
 
 +
<div id='img-22'></div>
 +
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
 +
|-
 +
|[[Image:Draft_Samper_805805629-test-cylt-fine07-pc-dtrho.png|330px|]]
 +
|[[Image:Draft_Samper_805805629-test-cylt-fine07-ep-dtrho.png|330px|Sidepressing of a cylinder &#8211; results obtained with the FIC algorithm for mesh of 22186 elements: a) pressure distribution, b) effective plastic distribution ( τ=\frac(∆t)²ρ ) ]]
 +
|- style="text-align: center; font-size: 75%;"
 +
| colspan="2" | '''Figure 22:''' Sidepressing of a cylinder &#8211; results obtained with the FIC algorithm for mesh of 22186 elements: a) pressure distribution, b) effective plastic distribution (<math> \tau =\frac{(\Delta t)^2}{\rho } </math>)
 +
|}
 +
 
 +
<div id="cite-1"></div>
 +
'''[1]'''  A.J. Chorin,  “A numerical solution for solving incompressible viscous flow problems” ''J. Comp. Phys.'', '''2''', 12&#8211;26, 1967.
 +
 
 +
<div id="cite-2"></div>
 +
'''[2]'''  W.R. Briley, S.S. Neerarambam and D.L. Whitfield, “Multigrid algorithm for three-dimensional incompressible high-Reynolds number turbulent flows”, ''AIAA Journal'',  '''33 (1)''', 2073&#8211;2079, 1995.
 +
 
 +
<div id='img-23'></div>
 +
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
 +
|-
 +
|[[Image:Draft_Samper_805805629-test-picture1.png|600px|]]
 +
|[[Image:Draft_Samper_805805629-test-picture50.png|600px|]]
 +
|-
 +
|[[Image:Draft_Samper_805805629-test-picture150.png|600px|]]
 +
|[[Image:Draft_Samper_805805629-test-picture300.png|600px|]]
 +
|-
 +
|[[Image:Draft_Samper_805805629-test-picture600.png|600px|]]
 +
|[[Image:Draft_Samper_805805629-test-picture801.png|600px|]]
 
|-
 
|-
|[[Image:Draft_Samper_805805629-cylhc-ep-new.png|228px|]]
+
| colspan="2"|[[Image:Draft_Samper_805805629-test-medit500.png|300px|Filling of a mould with a casted metal using the PFEM]]
|[[Image:Draft_Samper_805805629-cylhc-hp-new.png|228px|Sidepressing of a cylinder, mixed formulation,  hexahedral mesh: (a) effective plastic strain; (b) pressure distribution ]]
+
 
|- style="text-align: center; font-size: 75%;"
 
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 7:''' Sidepressing of a cylinder, mixed formulation,  hexahedral mesh: (a) effective plastic strain; (b) pressure distribution
+
| colspan="2" | '''Figure 23:''' Filling of a mould with a casted metal using the PFEM
 
|}
 
|}

Revision as of 14:41, 22 October 2018


Abstract

The paper describes some recent developments in finite element and particle methods for analysis of a wide range of bulk forming processes. The developments include new stabilized linear triangles and tetrahedra using finite calculus and a new procedure combining particle methods and finite element methods. Applications of the new numerical methods to casting, forging and other bulk metal forming problems and mixing processes are shown.

keywords Bulk forming processes, stabilized finite element method, particle method, particle finite element method, mixing processes.

1 INTRODUCTION

The development of efficient and robust numerical methods for analysis of bulk forming problems has been a subject of intensive research in recent years [1–7]. Many of these problems require the solution of incompressible fluid flow situations (such as in mould filling problems) whereas in other cases (such as forging, rolling, extrusion, etc.) the numerical method must be able to account for the quasi/fully incompressible behaviour induced by the large plastic deformation. The solution of these problems has motivated the development of the so called stabilized numerical methods overcoming the two main sources of instability in the analysis of incompressible continua, namely those originated by the high values of the convective terms in fluid flow situtions and those induced by the difficulty in satisfying the incompressibility condition.

Different approaches to solve both type of problems in the context of the finite element method (FEM) have been recently developed [8]. Traditionally, the underdiffusive character of the Galerkin FEM for high convection flows has been corrected by adding some kind of artificial viscosity terms to the standard Galerkin equations [8,9].

A popular way to overcome the problems with the incompressibility constraint in the FEM is by introducing a pseudo-compressibility in the continuum and using implicit and explicit algorithms ad hoc such as artificial compressibility schemes [10] and preconditioning techniques [11]. Other FEM schemes with good stabilization properties for the convective and incompressibility terms in fluid flows are based in Petrov-Galerkin (PG) techniques. The background of PG methods are the non-centred (upwind) schemes for computing the first derivatives of the convective operator in FD and FV methods [8,9,12]. A general class of Galerkin FEM has been developed where the standard Galerkin variational form is extended with adequate residual-based terms in order to achieve a stabilized numerical scheme. Among the many FEM of this kind we can name the Streamline Upwind Petrov Galerkin (SUPG) method [8,13–16], the Galerkin Least Square (GLS) method [8,17], the Taylor-Galerkin method [18], the Characteristic Galerkin method [19] and its variant the Characteristic Based Split (CBS) method [20,21], the pressure gradient operator method [22] and the Subgrid Scale (SS) method [23]. A good review of these methods can be found in [24]. Extensions of the CBS and SS methods to treat incompressible problems in solid mechanics are reported in [25,26,58] and [27–29], respectively.

In this paper a different class of stabilized FEM for quasi and fully incompressible fluid and solid materials applicable to a wide range of bulk forming problems is presented. The starting point is the modified governing differential equations of the continuous problem formulated via a finite calculus (FIC) approach [30,31]. The FIC method is based in invoking the classical balance (or equilibrium) laws in a domain of finite size. This introduces naturally additional terms in the differential equations of infinitesimal continuum mechanics which are a function of the balance domain dimensions. The new terms in the modified governing equations provide the necessary stabilization to the discrete equations obtained via the standard Galerkin FEM. One of the main advantages of the FIC formulation versus other alternative approaches (such as mixed FEM, etc.) is that it allows to solve incompressible fluid problems using low order finite elements (such as linear triangles and tetrahedra) with equal order approximations for the velocity and pressure variables [32–35]. The FIC formulation has been successfully used for analysis of fully or quasi incompressible solids [36,37].

The layout of the paper is the following. In the next section the basic FIC equations for incompressible flow problems formulated in an Eulerian frame are presented. The finite element discretization is introduced and the resulting discretized equations are detailed. A fractional step scheme for the transient solution is presented.

The stabilized Eulerian formulation is extended to account for thermal effects and the transport of the free surface which are needed for mould filling processes.

The following sections outline the FIC formulation for analysis of quasi/fully incompressible solids using a Lagrangian description and linear triangles and tetrahedra. An explicit algorithm for integrating in time the equations of motion of elasto-plastic solids in large-strain problems involving frictional contact is described. Examples of application to a casting problem and some bulk forming processes are presented.

In the last part of the paper a Lagrangian formulation for fluid flow analysis is presented as a straightforward extension of the formulation for solid mechanics. The procedure, called the Particle Finite Element Method (PFEM) [38–40], treats the mesh nodes in the fluid and solid domains as dimensionless particles which can freely move, an even separate from the main fluid domain, representing, for instance, the effect of liquid drops. A finite element mesh connects the nodes defining the discretized domain where the governing equations are solved in the standard FEM fashion. The main advantage of the Lagrangian flow formulation is that the convective terms do not enter in the fluid equations. The difficulty is however transferred to the problem of adequately (and efficiently) moving the mesh nodes. The final examples show that the PFEM is a promising method to solve mould filling and casting problems, material mixing processes and many other bulk metal forming problems involving the interaction between solids and fluids which can be treated with the same formulation.

2 FIC EQUATIONS FOR VISCOUS INCOMPRESSIBLE FLOWS

The FIC governing equations for a viscous incompressible fluid can be written in an Eulerian frame of reference as [30–34]

Momentum

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): r_{m_i} - \underline{{1\over 2} h_j{\partial r_{m_i} \over \partial x_j}}{1\over 2} h_j{\partial r_{m_i} \over \partial x_j}=0 \quad \hbox{in }\Omega
(1)

Mass balance

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): r_d - \underline{{1\over 2} h_j {\partial r_d \over \partial x_j}}{1\over 2} h_j {\partial r_d \over \partial x_j}=0 \quad \hbox{in }\Omega
(2)

where

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): r_{m_i} = \rho \left({\partial v_i \over \partial t}+v_j{\partial v_i \over \partial x_j}\right)+ {\partial p \over \partial x_i}- {\partial s_{ij} \over \partial x_j}-b_i (3)
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): r_d = {\partial v_i \over \partial x_i}\qquad i,j = 1, n_d (4)

Above Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \Omega }

is the analysis domain which can evolve with time, Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle n_d}
is the number of space dimensions (Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle n_d=3}
for 3D problems), Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle v_i}
is the velocity along the ith global axis, Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \rho }
is the (constant) density of the fluid, Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle p}
is the absolute pressure (defined positive in compression), Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle b_i}
are the body forces and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle s_{ij}}
are the  deviatoric stresses related to the viscosity Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \mu }
by the standard expression
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): s_{ij}=2\mu \left(\dot \varepsilon _{ij} - \delta _{ij} {1\over 3} {\partial v_k \over \partial x_k}\right)
(5)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \delta _{ij}}

is the Kronecker delta and the strain rates Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \dot \varepsilon _{ij}}
are
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \dot \varepsilon _{ij}={1\over 2} \left({\partial v_i \over \partial x_j}+{\partial v_j \over \partial x_i}\right)
(6)

The boundary conditions are written in the FIC approach as

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): n_j \sigma _{ij} -t_i^p + \underline{{1\over 2} h_j n_j r_{m_i}}{1\over 2} h_j n_j r_{m_i}=0 \quad \hbox{on }\Gamma _t
(7)
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): v_j - v_j^p =0 \quad \hbox{on }\Gamma _u
(8)

and the initial condition is Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle v_j =v_j^0}

for Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle t=t_0}

.

Summation convention for repeated indexes in products and derivatives is used unless otherwise specified.

In Eqs.(7) and (8) Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle t_i^p}

and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle v_j^p}
are surface tractions and prescribed velocities on the boundaries Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \Gamma _t}
and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \Gamma _u}

, respectively, Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle n_j}

are the components of the unit normal vector to the boundary and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \sigma _{ij}}
are the total stresses given by Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \sigma _{ij}=s_{ij}-\delta _{ij}p}

.

Eqs. (1) and (2) are obtained by invoking the classical balance equations in fluid mechanics in a domain of finite size and retaining higher order terms [30]. The Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle h_i's}

in above equations are characteristic lengths of the domain where the balance of momentum and mass is enforced. In Eq.(7) these lengths define the domain where equilibrium of boundary tractions is established. In the discretized problem the Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle h_{i}}
coincide with a typical element dimension, as described in Section 5. Note that by making Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle h_i =0}
in these equations the standard infinitesimal form of the fluid mechanics equations is recovered [8,9,24].

Eqs.(1)–(8) are the starting point for deriving stabilized FEM for solving the incompressible Navier-Stokes equations using equal order interpolation for the velocity and pressure variables [32–35]. Application of the FIC formulation to meshless analysis of fluid flow problems using the finite point method can be found in [42].

Remark

In most metal forming processes the viscosity Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \mu }

is a non linear function of the strain rate and the yield stress of the material (which also depends on the temperature) [8,43,44]. This dependence adds another non linearity to the problem.

Remark

In Eqs.(1)–(7) Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \Omega }

and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \Gamma _t}
respectively denote the actual volume and the boundary of the domain where the governing equations are solved at each instant of the forming processes. This domain is considered here to be fixed in space (Eulerian approach). Moving free surfaces, such as in the case of mold filling processes, can be modelled by using standard level set and volume of fluid (VOF) techniques [47,48]. An alternative is to use an arbitrary Lagragian-Eulerian (ALE) description or even a fully Lagrangian description to follow the motion of the particles during the forming process.  The Lagrangian description is typical in the case of solid mechanics problems (Section 5). A particular Lagrangian formulation for fluid flow problems involving large motion of the free surface, named the Particle Finite Element Method, is presented in Section 8.

2.1 Stabilized integral forms

From the momentum equations it can be obtained [32–34]

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\partial r_d \over \partial x_i}\simeq {h_j\over 2a_i} {\partial r_{m_i} \over \partial x_j}\quad ,\quad \hbox{no sum in }i
(9)

where

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): a_i = {2\mu \over 3} +{\rho v_i h_i\over 2}\quad ,\quad \hbox{no sum in }i
(10)

Substituting Eq.(9) into Eq.(2) and retaining the terms involving the derivatives of Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle r_{m_i}}

with respect to Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle x_i}
only, leads to the following expression for the stabilized mass balance equation
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): r_d - \sum \limits _{i=1}^{n_d} \tau _i {\partial r_{m_i} \over \partial x_i}=0
(11)

with

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \tau _i = \left({8\mu \over 3h_i^2}+{2\rho v_i\over h_i}\right)^{-1}
(12)

The Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \tau _i} 's in Eq.(12), when scaled by the density, are termed in the stabilization literature intrinsic time parameters.

The weighted residual form of the momentum and mass balance equations (Eqs.(1) and (11)) is written as

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \int _\Omega \delta v_i \left[r_{m_i} - {h_j\over 2} {\partial r_{m_i} \over \partial x_j}\right]d\Omega + \int _{\Gamma _t} \delta v_i (\sigma _{ij} n_j - t_i^p + {h_j \over 2} n_j r_{m_i}) d\Gamma =0
(13)
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \int _\Omega q \left[r_d - \sum \limits _{i=1}^{n_d} \tau _i {\partial r_{m_i} \over \partial x_i}\right]d\Omega =0
(14)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \delta v_i}

and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle q}
are arbitrary weighting functions representing virtual velocities and virtual pressure fields. Integrating by parts the Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle r_{m_i}}
terms  leads to
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \int _\Omega \delta v_i r_{m_i}d\Omega + \int _{\Gamma _t} \delta v_i (\sigma _{ij} n_j - t_i)d\Gamma + \int _{\Omega } {h_j\over 2}{\partial \delta v_i \over \partial x_j} r_{m_i} d\Omega =0
(15a)
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \int _\Omega q r_d d\Omega + \int _\Omega \left[\sum \limits _{i=1}^{n_d}\tau _i {\partial q \over \partial x_i}r_{m_i} \right]d\Omega - \int _\Gamma \left[\sum \limits _{i=1}^{n_d} q \tau _i n_i r_{m_i}\right]d\Gamma =0
(15b)

We will neglect hereonwards the third integral in Eq.(15b) by assuming that Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle r_{m_i}}

is negligible on the boundaries. The deviatoric stresses and the pressure terms in the first integral of Eq.(15a) are integrated by parts in the usual manner. The resulting momentum and mass balance equations are
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \begin{array}{r} \displaystyle \int _\Omega \left[\delta v_i\rho \left({\partial v_i \over \partial t}+v_j {\partial {v_i} \over \partial x_j}\right)+ {\partial \delta v_i \over \partial x_j}\left(\mu {\partial v_i \over \partial x_j}-\delta _{ij}p \right)\right]d\Omega - \int _{\Omega } \delta v_i b_i d\Omega - \int _{\Gamma _t} \delta v_i t_i^p d\Gamma +\qquad \\ \displaystyle + \int _{\Omega } {h_j\over 2}{\partial \delta v_i \over \partial x_j} r_{m_i} d\Omega =0\qquad \end{array}
(16a)
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \int _\Omega q {\partial v_i \over \partial x_i} d\Omega + \int _\Omega \left[\sum \limits _{i=1}^{n_d} \tau _i {\partial q \over \partial x_i} r_{m_i}\right]d\Omega =0
(16b)

In the derivation of the viscous term in Eq.(16a) we have used the following identity (prior to the integration by parts)

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\partial s_{ij} \over \partial x_j}=2\mu {\partial \varepsilon _{ij} \over \partial x_j}=\mu {\partial ^2v_i\over \partial x_j \partial x_j}
(17)

Eq.(17) is identically true for the exact incompressible limit Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \left({\partial v_i \over \partial x_i}=0\right)} .

2.2 Convective and pressure gradient projections

The computation of the residual terms can be simplified if we introduce now the convective and pressure gradient projections Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle c_i}

and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \pi _i}

, respectively defined as

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \begin{array}{l}\displaystyle c_i = r_{m_i} -\rho v_j {\partial v_i \over \partial x_j}\\ \displaystyle \pi _i = r_{m_i} - {\partial p \over \partial x_i} \end{array}
(18)

We can express Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle r_{m_i}}

in  Eqs.(16a) and (16b) in terms of Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle c_i}
and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \pi _i}

, respectively which then become additional variables. The system of integral equations is now augmented in the necessary number of equations by imposing that the residual Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle r_{m_i}}

vanishes (in average sense) for both forms given by Eqs.(18). This gives the final system of governing equation as:
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \int _\Omega \left[\delta v_i\rho \left({\partial v_i \over \partial t}+v_j {\partial {v_i} \over \partial x_j}\right)+ {\partial \delta v_i \over \partial x_j}\left(\mu {\partial v_i \over \partial x_j}-\delta _{ij}p \right)\right]d\Omega - \int _{\Omega } \delta v_i b_i d\Omega - \int _{\Gamma _t} \delta v_i t_i^p d\Gamma +
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): + \int _{\Omega } {h_k\over 2}{\partial (\delta v_i) \over \partial x_k} \left(\rho v_j {\partial {v_i} \over \partial x_j} + c_i\right)d\Omega =0\qquad
(19)
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \int _\Omega q {\partial v_i \over \partial x_i} d\Omega + \int _\Omega \sum \limits _{i=1}^{n_d} \tau _i {\partial q \over \partial x_i} \left({\partial p \over \partial x_i}+\pi _i\right)d\Omega =0
(20)
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \int _\Omega \delta c_i \rho \left(\rho v_j {\partial {v_i} \over \partial x_j} + c_i\right)d\Omega =0 \qquad \hbox{no sum in }i
(21)
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \int _\Omega \delta \pi _i \tau _i \left({\partial p \over \partial x_i}+\pi _i\right)d\Omega =0\qquad \hbox{no sum in }i
(22)

with Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle i,j,k=1,n_d} . In Eqs.(21) and (22) Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \delta c_i}

and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \delta \pi _i}
are appropriate weighting functions and the Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \rho }
and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \tau _i}
weights are introduced for convenience.

We note that accounting for the convective and pressure gradient projections enforces the consistency of the formulation as it ensures that the stabilization terms in Eqs.(19–22) have a residual form which vanishes for the “exact” solution. Neglecting these terms lowers the accuracy of the numerical solution and it makes the formulation more sensitive to the value of the stabilization parameters as shown in [34,36,37].

3 FINITE ELEMENT DISCRETIZATION

We choose Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle C^\circ }

continuous linear interpolations of the velocities, the pressure, the convection projections Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle c_i}
and the pressure gradient projections Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \pi _i}
over three node triangles (2D) and four node tetrahedra (3D) [8]. The linear interpolations are written as
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \begin{array}{l}\displaystyle v_i = N^k \bar v_i^k \quad , \quad p = N^k \bar p^k\\ \displaystyle c_i = N^k \bar c_i^k \quad , \quad \pi _i = N^k \bar \pi _i^k \end{array}
(23)

where the sum goes over the number of nodes of each element Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle n}

(Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle n=3/4}
for triangles/tetrahedra), Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \bar {(\cdot )}^k}
denotes nodal variables and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle N^k}
are the linear shape functions [8].

Substituting the approximations (23) into Eqs.(19–22) and choosing the Galerking form with Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \delta v_i =q=\delta c_i=\delta \pi _i =N^i}

leads to following system of discretized equations
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \displaystyle {\boldsymbol M}\dot{\bar {\boldsymbol v}} + {\boldsymbol H} \bar {\boldsymbol v} - {\boldsymbol G}\bar {\boldsymbol p}+{\boldsymbol C}\bar {\boldsymbol c}={\boldsymbol f}
(24a)
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \displaystyle {\boldsymbol G}^T \bar {\boldsymbol v} + \hat{\boldsymbol L}\bar {\boldsymbol p}+{\boldsymbol Q}\bar {\boldsymbol \pi }={\boldsymbol 0}
(24b)
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \displaystyle \hat {\boldsymbol C}\bar{\boldsymbol v}+ {\boldsymbol M}\bar {\boldsymbol c}={\boldsymbol 0}
(24c)
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \displaystyle {\boldsymbol Q}^T \bar {\boldsymbol p} + \hat {\boldsymbol M}\bar {\boldsymbol \pi }={\boldsymbol 0}
(24d)

The form of the different matrices is given in the Appendix.

The solution in time of the system of Eqs.(24) can be written in general form as

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\boldsymbol M} \displaystyle{1\over \Delta t} (\bar{\boldsymbol v}^{n+1}-\bar{\boldsymbol v}^{n}) + {\boldsymbol H}^{n+\theta }\bar{\boldsymbol v} ^{n+\theta } - {\boldsymbol G}\bar{\boldsymbol p}^{n+\theta }+ {\boldsymbol C}^{n+\theta }\bar {\boldsymbol c}^{n+\theta }={\boldsymbol f}^{n+\theta }
(25a)
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\boldsymbol G}^T \bar {\boldsymbol v}^{n+\theta }+\hat {\boldsymbol L}^{n+\theta } \bar {\boldsymbol p}^{n+\theta }+{\boldsymbol Q}\bar {\boldsymbol \pi }^{n+\theta }={\boldsymbol 0}
(25b)
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \hat {\boldsymbol C}^{n+\theta } \bar {\boldsymbol v}^{n+\theta }+ {\boldsymbol M} \bar {\boldsymbol c}^{n+\theta }={\boldsymbol 0}
(25c)
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\boldsymbol G}^T \bar {\boldsymbol p}^{n+\theta }+\hat {\boldsymbol M}^{n+\theta }\bar {\boldsymbol \pi }^{n+\theta }={\boldsymbol 0}
(25d)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {\boldsymbol H}^{n+\theta }={\boldsymbol H} (\bar{\boldsymbol v}^{n+\theta })} , Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \bar{\boldsymbol v}^{n+\theta }}

are the velocities evaluated at time Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle n+\theta }
and the parameter Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \theta \in [0,1]}

. The direct monolitic solution of Eqs.(25) is possible using an adequate iterative scheme. However, we have found more convenient to use a fractional step method as described in the next section.

3.1 Fractional step method

A fractional step scheme is derived by noting that the discretized momentum equation (25a) can be split into the two following equations

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\boldsymbol M} \displaystyle{1\over \Delta t} (\tilde{\boldsymbol v}^{n+1}-\bar{\boldsymbol v}^{n}) + {\boldsymbol H}^{n+\theta }\bar{\boldsymbol v} ^{n+\theta } - \alpha {\boldsymbol G}\bar{\boldsymbol p}^{n} + {\boldsymbol C}^{n+\theta }\bar {\boldsymbol c}^{n+\theta }={\boldsymbol f}^{n+\theta }
(26)
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\boldsymbol M} \displaystyle{1\over \Delta t} (\bar{\boldsymbol v}^{n+1}-\tilde{\boldsymbol v}^{n+1})- {\boldsymbol G}(\bar{\boldsymbol p}^{n+1}-\alpha \bar {\boldsymbol p}^{n})={\boldsymbol 0}
(27)

In above equation Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \tilde{\boldsymbol v}^{n+1}}

is a predicted value of the velocity at time Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle n+1}
and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \alpha }
is a variable whose values of interest are zero and one. For Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \alpha =0}
(first order scheme) the splitting error is of order Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle  0 (\Delta t)}

, whereas for Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \alpha =1}

(second order scheme) the error is of order Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle 0 (\Delta t^2)}
[45].

Eqs.(26) and (27) are completed with the following three equations emanating from Eqs.(25b-d)

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\boldsymbol G}^T\bar{\boldsymbol v}^{n+1}+\hat {\boldsymbol L}^n \bar{\boldsymbol p}^{n+1}+{\boldsymbol Q}\bar {\boldsymbol \pi }^{n}= {\boldsymbol 0}
(28a)
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \hat {\boldsymbol C}^{n+1} \bar{\boldsymbol v}^{n+1} + {\boldsymbol M} \bar {\boldsymbol c}^{n+1}= {\boldsymbol 0}
(28b)
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\boldsymbol Q}^T \bar{\boldsymbol p}^{n+1}+ \hat {\boldsymbol M}^{n+1} \bar {\boldsymbol \pi }^{n+1}= {\boldsymbol 0}
(28c)

The value of Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \bar{\boldsymbol v}^{n+1}}

obtained from Eq.(27) is substituted into Eq.(28a) to give
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\boldsymbol G}^T\tilde{\boldsymbol v}^{n+1}+ \Delta t {\boldsymbol G}^T {\boldsymbol M}^{-1} {\boldsymbol G} (\bar{\boldsymbol p}^{n+1} - \alpha \bar{\boldsymbol p}^n)+ \hat {\boldsymbol L}^n {\boldsymbol p}^{n+1}+{\boldsymbol Q}\bar {\boldsymbol \pi }^n={\boldsymbol 0}
(29)

The product Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {\boldsymbol G}^T {\boldsymbol M}^{-1} {\boldsymbol G}}

can be approximated by a laplacian matrix, i.e.
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\boldsymbol G}^T {\boldsymbol M}^{-1} {\boldsymbol G}={1\over \rho } {\boldsymbol L}\quad \hbox{with }L^{ab}=\int _{ \Omega ^e} {\boldsymbol \nabla }^T N^a {\boldsymbol \nabla } N^b d\Omega
(30)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle L^{ab}}

are the element contributions to Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {\boldsymbol L}}
(see Appendix).

The steps of the fractional step scheme chosen here are:

Step 1

The fractional nodal velocities Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \tilde {\boldsymbol v}^{n+1}}

can be  explicitely computed from Eq.(26) by
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \tilde{\boldsymbol v}^{n+1} = \bar{\boldsymbol v}^{n} - \Delta t {\boldsymbol M}^{-1}_d [\tilde{\boldsymbol H}^{n}\bar{\boldsymbol v}^{n} - \alpha {\boldsymbol G} \bar{\boldsymbol p}^n + {\boldsymbol C}^{n}\bar{\boldsymbol c}^{n} - {\boldsymbol f}^{n}]
(31)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {\boldsymbol M}_d}

is the diagonal form of M obtaining by lumping the row terms into the corresponding diagonal terms.


Step 2 Compute Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \bar{\boldsymbol p}^{n+1}}

from Eq.(29) as
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \bar{\boldsymbol p}^{n+1}= -[\hat{\boldsymbol L}^n + {\Delta t \over \rho } {\boldsymbol L}]^{-1} [{\boldsymbol G}^T\tilde{\boldsymbol v}^{n+1} - \alpha {\Delta t \over \rho }{\boldsymbol L}\bar{\boldsymbol p}^n +{\boldsymbol Q} \bar{\boldsymbol \pi }^{n}]
(32)

Step 3 Compute Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \bar{\boldsymbol v}^{n+1}}

explicitely from Eq.(28a) as
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \bar {\boldsymbol v}^{n+1}=\tilde{\boldsymbol v}^{n+1}+ \Delta t {\boldsymbol M}_d^{-1} {\boldsymbol G} (\bar {\boldsymbol p}^{n+1}- \alpha \bar {\boldsymbol p}^n)
(33)


Step 4 Compute Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \bar{\boldsymbol c}^{n+1}}

explicitely from Eq.(28b) as
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \bar{\boldsymbol c}^{n+1}=- {\boldsymbol M}_d^{-1}\hat {\boldsymbol C}^{n+1}\bar{\boldsymbol v}^{n+1}
(34)

Step 5 Compute Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \bar{\boldsymbol \pi }^{n+1}}

explicitely from Eq.(28c) as
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \bar{\boldsymbol \pi }^{n+1}=- \hat {\boldsymbol M}_d^{-1} {\boldsymbol Q}^T \bar {\boldsymbol p}^{n+1}
(35)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \hat {\boldsymbol M}_d}

is the lumped form of Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \hat {\boldsymbol M}}

. A standard diagonal lumping procedure based in summing up the terms of each row has been used.

Note that all steps can be solved explicitely except for the computation of the pressure in Eq.(32) which requires to invert the sum of two laplacian matrices. This can be effectively performed using an iterative solution scheme such as the conjugate gradient method.

Above algorithm has improved stabilization properties versus the standard segregation methods due to the introduction of the laplacian matrix Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \hat {\boldsymbol L}}

in Eq.(32). This matrix emanates from the FIC formulation.

The boundary conditions are applied as follows. No condition is applied in the computation of the fractional velocities Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \tilde{\boldsymbol v}^{n+1}}

in Eq.(31). The prescribed velocities at the boundary are applied when solving for Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \bar{\boldsymbol v}^{n+1}}
in  step 3. The prescribed pressures at the boundary are imposed by making Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \bar{\boldsymbol p}^n}
equal to the prescribed pressure values when solving Eq.(32).

3.2 Stokes flow

Many metal forming processes can be simulated with the assumption that the convective terms are negligible (Stokes flow) [8,43,44]. The Stokes FIC formulation can be readily obtained simply by neglecting the convective terms in the Navier-Stokes formulation presented earlier. This also implies neglecting the convective stabilization terms in the momentum equations and, consequently, the convective projection variables are not larger necessary. Also the intrinsic time parameters Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \tau _i}

take now the simpler form (see Eq.(12)):
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \tau _i={3h_i^2\over 8\mu }
(36)

We note again that in metal forming problems the viscosity Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \mu }

will typically be a function of the strain rate and the yield stress of the material [8,43,44].

The resulting discretized system of equations can be written as (see Eqs.(28))

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \begin{array}{l}\displaystyle {\boldsymbol M}\dot{\bar {\boldsymbol v}} + {\boldsymbol K}\bar{\boldsymbol v} - {\boldsymbol G}\bar {\boldsymbol p}={\boldsymbol f}\\ \displaystyle {\boldsymbol G}^T \bar {\boldsymbol v} + \hat{\boldsymbol L}\bar {\boldsymbol p}+{\boldsymbol Q}\bar {\boldsymbol \pi }={\boldsymbol 0}\\ \displaystyle {\boldsymbol Q}^T \bar {\boldsymbol p} + \hat {\boldsymbol M}\bar {\boldsymbol \pi }={\boldsymbol 0} \end{array}
(37)

The algorithm of previous section can now be implemented.

The steady-state form of Eqs.(37) can be expressed in matrix form as

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \left[\begin{matrix}{\boldsymbol K}&-{\boldsymbol G}&{\boldsymbol 0}\\ -{\boldsymbol G}^T & - \hat{\boldsymbol L} &-{\boldsymbol Q}\\ {\boldsymbol 0}& -{\boldsymbol Q}^T&-\hat {\boldsymbol M}\\\end{matrix}\right]\left\{\begin{matrix}\bar {\boldsymbol v}\\ \bar {\boldsymbol p}\\ \bar {\boldsymbol \pi }\\\end{matrix}\right\}= \left\{\begin{matrix}{\boldsymbol f}\\ {\boldsymbol 0}\\{\boldsymbol 0}\\ \end{matrix}\right\}
(38)

The system is symmetric and always positive definite and therefore leads to a non singular solution. This property holds for any interpolation function chosen for Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \bar {\boldsymbol v},\bar {\boldsymbol p}}

and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \bar {\boldsymbol \pi }}

, therefore overcoming the Babuŝka-Brezzi (BB) restrictions [8].

A reduced velocity-pressure formulation can be obtained by elliminating the Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \bar \pi }

variables from the last row of Eq.(38)  [37].

The FIC formulation for Stokes flows is applicable for analysis of quasi/fully incompressible solids. An analogous formulation based on solid mechanics concepts is derived in Section 6.

4 THERMAL-COUPLED FLOWS. TREATMENT OF THE FREE BOUNDARY

4.1 Thermal coupled flow

The effect of temperature can be easily introduced by solving the equation for heat transport coupled to the fluid flow equations.

The equation of balance of heat is written in the FIC formulation as [30]

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): r_\phi - {h_j\over 2} {\partial r_\phi \over \partial x_j}=0\qquad ,\quad j=1,n_d
(39a)

with

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): r_\phi :=-\rho c \left({\partial \phi \over \partial t}+v_i {\partial \phi \over \partial x_i}\right)+ {\partial \over \partial x_k} \left(k{\partial \phi \over \partial x_k}\right)+Q
(39b)

In above Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \phi }

is the temperature, Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle c}
and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle k}
are the specific heat and the thermal conductivity of the material, respectively, Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle Q}
is the heat source and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle h_j}
are the characteristic length distances which are typical of the FIC formulation [30].

Eq.(39a) is completed with the Dirichlet and Neuman boundary conditions for the heat problem. For details see [30,46].

The convective velocities Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle v_i}

in Eq.(39b) are provided by the solution of the fluid flow problem. As usual in metal forming processes, the heat source Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle Q}
is a function of the mechanical work generated in the flow of the material during the forming process. The temperature field affects in turn the flow viscosity via its dependence with the yield stress which is very sensitive to the temperature changes. The solution of the heat transfer equation is therefore fully coupled with that of the fluid flow problem [44].

4.2 Treatment of the free boundary motion

For the treatment of the free boundary we use a standard volume of fluid (VOF) technique, also known as pseudo-concentration method or level set technique [22,47,48]. In the VOF method the motion of the free boudary is followed by solving the following transport equation (written using the FIC formulation)

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): r_\beta - \underline{{\bar h_j\over 2} {\partial r_\beta \over \partial x_j}}{\bar h_j\over 2} {\partial r_\beta \over \partial x_j}=0 \quad \hbox{in }\Omega
(40a)

where

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): r_\beta :={\partial \beta \over \partial t}+v_k {\partial \beta \over \partial x_k}
(40b)

In above Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \beta }

is an auxiliary variable which takes a value equal to one on the free surface and zero elsewhere. The solution of Eq.(40a) in time in the discretized fluid-air domain allows to track the free surface which is characterized by the contours of Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \beta }
which have a unit value.

The underlined term in Eq.(40a) emerges from the FIC formulation [33]. Again, the Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \bar h_j} 's in Eq.(40a) are characteristic length distances of the order of the element size. This term introduces the necessary stabilization in the transient solution of Eq.(40a).

4.3 Fully coupled algorithm for mould filling problems

The analysis of mould filling problems typically requires the solution of the coupled-thermal flow problem accounting for the transport of the free surface.

The simplest explicit algorithm is as follows:

  • For each time step:
  • Compute the velocities and the pressure in the fluid Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle (\bar{v}^{n+1}_i, \bar{p}^{n+1})}
using the fractional step scheme of Section 3.1.
  • Compute the nodal temperatures Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle (\phi ^{n+1})}
solving Eq.(39a).
  • Compute the new position of the free surface Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle (\beta ^{n+1})}
solving Eq.(40a).

A number of implicit versions of the algorithm are possible and they all involve an iteration loop within each time step until convergence for the flow variables, the temperature and the free surface position is found. For details see [47,48].

5 COMPUTATION OF THE CHARACTERISTIC LENGTHS

The evaluation of the stabilization parameters is one of the crucial issues in stabilized methods. Excellent results have been obtained in all problems solved using linear tetrahedra with the characteristic length vector defined by

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\boldsymbol h}=h_s {{\boldsymbol v}\over {v}}+h_{c} {{\boldsymbol \nabla } v\over \vert{\boldsymbol \nabla }v\vert }
(41a)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle v=\vert {\boldsymbol v}\vert }

and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle h_s}
and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle h_{c}}
are the “streamline” and “cross wind” contributions given by
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): h_s=\max ({\boldsymbol l}^T_j {\boldsymbol v})/{v}
(41b)
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): h_{c}=\max ({\boldsymbol l}^T_j {\boldsymbol \nabla }v)/ \vert {\boldsymbol \nabla }v\vert \quad , \quad j=1,n_s
(41c)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {\boldsymbol l}_j}

are the vectors defining the element sides (Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle n_s=6}
for tetrahedra).

As for the free surface equation the following value of the characteristic length vector Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \bar{\boldsymbol h}}

in Eq.(40a) has been taken
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \bar{\boldsymbol h} =\bar h_s {{\boldsymbol v}\over {v}}+\bar h_c {{\boldsymbol \nabla }\beta \over \vert {\boldsymbol \nabla }\beta \vert }
(42a)

The streamline parameter Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \bar h_s}

has been obtained by Eq.(41b) whereas the cross wind parameter Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \bar h_c}
has been computed by
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \bar h_c = \max [{\boldsymbol l}_j^T {\boldsymbol \nabla }\beta ] {1\over \vert {\boldsymbol \nabla }\beta \vert } \quad ,\quad j=1,2,3
(42b)

The cross-wind terms in eqs.(41a) and (42a) account for the effect of the gradient of the solution in the stabilization parameters. This is a standard assumption in most “shock-capturing” stabilization procedures [8,24].

6 FIC FORMULATION FOR QUASI/FULLY INCOMPRESSIBLE SOLIDS

6.1 Equilibrium equations

As usual the governing equations of solid mechanics are written in a Lagrangian reference frame. Following the arguments of Section 2 the equilibrium equations for a solid are written using the FIC technique as [37]

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {r}_{i}- \frac{{h}_k}{2}{\partial {r}_{i} \over \partial {x}_k}=0 \quad \mbox{in} \;\,\Omega \qquad k=1,n_d
(43)

where for the dynamic case

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {r}_{i}: =-\rho {\partial ^2 u_i\over \partial t^2} + {\partial {\sigma }_{ij} \over \partial {x}_j}+{b}_i\qquad , \quad j=1,n_d
(44)

In Eqs.(43) and (44) Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle u_i}

are the displacements, Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\sigma }_{ij}

and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {b}_i are the stresses and the body forces, respectively and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {h}_k are characteristic length distances of an arbitrary prismatic domain where equilibrium of forces is considered.

Equations (43) and (44) are completed with the boundary conditions on the displacements Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle u_i}


Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {u}_i-{u}_i^p =0\quad \mbox{on} \;\,\Gamma _u
(45)

and the equilibrium of surface tractions

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\sigma }_{ij}{n}_j- t_i^p - \underline{\frac{1}{2}{h}_k{n}_k{r}_{i}}\frac{1}{2}{h}_k{n}_k{r}_{i} =0\quad \mbox{on} \;\,\Gamma _t
(46)

In the above Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {u}_i^p}

and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle t_i}
are prescribed displacements and tractions over the boundaries Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \Gamma _u}

and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \Gamma _t} , respectively and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle n_i}

are the components of the unit normal vector.

Note that for consistency with the fluid flow equations of previous section and differently from the tradition in solid mechanics, the compression pressure has been taken as positive.

6.2 Stabilized pressure constitutive equations

For simplicity the treatment of the constitutive equations will be explained for the linear elastic model. The approach extends naturally to the non linear elasto-plastic/viscoplastic constitutive equations typical of metal forming problems.

As usual in quasi-incompressible problems the stresses are split into deviatoric and volumetric (pressure) parts

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\sigma }_{ij}= s_{ij} -p{\delta }_{ij}
(47)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\delta }_{ij} is the Kronnecker delta function. The linear elastic constitutive equations for the deviatoric stresses Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle s_{ij}}

are written as

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): s_{ij}\!\!=\!\! 2G \left({\varepsilon }_{ij}-\frac{1}{3}{\varepsilon }_{v}{\delta }_{ij}\right)
(48)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle G} is the shear modulus,

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\varepsilon }_{ij}= \frac{1}{2} \left( {\partial {u}_i \over \partial {x}_j}+ {\partial {u}_j \over \partial x_i} \right) \quad \mbox{and} \quad {\varepsilon }_{v}={\varepsilon }_{ii}\,.
(49)

The constitutive equation for the pressure Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle p}

can be written  for an arbitrary domain of finite size  of volume Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle V}
as

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {1\over K} p^{av}= - {\Delta V\over V}
(50)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle K}

is the bulk modulus of the material and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle p^{av}}
is the average value of the pressure over  domain Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle V}

.

The value of Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle p^{av}}

can be approximated as [37]

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): p^{av}: ={1\over V} \int _V p\,dV =p - {h_k\over 2}{\partial p \over \partial x_k} + O (h_k)^2\quad ,\quad k=1,n_d
(51)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle p}

is the pressure at an arbitrary point within the domain Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle V}
and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle h_k}
are characteristic lengths of such a domain.

The ratio Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {\Delta V\over V}}

can be expressed as

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\Delta V\over V} = {\varepsilon }_{v}- {h_k\over 2}{\partial {\varepsilon }_{v} \over \partial x_k} + O (h_k)^2\quad ,\quad k=1,n_d
(52)

Substituting Eqs.(51) and (52) into Eq.(50) and neglecting second order terms in Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle h_k}

gives the FIC constitutive equation for the pressure as

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \left({p\over K} + {\varepsilon }_{v}\right)- \frac{{h}_k}{2}{\partial \over \partial {x}_k} \left({p\over K} + {\varepsilon }_{v}\right)=0\, ,\quad k=1,n_d
(53)

Note that for Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle h_k \to 0}

the standard relationship between the pressure and the volumetric strain of the infinitesimal theory Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle (p=-K\varepsilon _v )}
is found.

For an incompressible material Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle K\to \infty }

and Eq.(53) yields
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \varepsilon _v - \frac{{h}_k}{2}{\partial {\varepsilon }_{v} \over \partial {x}_k}=0
(54)

Eq.(54) expresses the limit incompressible behaviour of the solid. This equation is identical to Eq.(2) for incompressible flow problems and there arises from the mass continuity condition [30,32].

By combining Eqs. (43), (44), (47), (48) and (53) a mixed displacement–pressure formulation can be written as

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): -\rho {\partial ^2 u_i\over \partial t^2} + {\partial \sigma _{ij} \over \partial x_i}+{b}_i-\frac{{h}_k}{2} {\partial r_i \over \partial {x}_k}=0
(55)

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \left(\frac{p}{K}+{\varepsilon }_{v}\right)-\frac{{h}_k}{2}{\partial \over \partial {x}_k}\left(\frac{p}{K}+{\varepsilon }_{v}\right)=0
(56)

From the observation of Eq.(55) we can obtain after some algebra [37]

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\partial \over \partial x_k}\left(\frac{p}{K}+{\varepsilon }_{v}\right)\simeq {3h_j\over 8G} {\partial r_k \over \partial x_j}
(57)

Substituting Eq.(57) into (53) the mass balance equation can be written as

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \left(\frac{p}{K}+{\varepsilon }_{v}\right)- \sum \limits _{i=1}^{n_d} \tau _i {\partial r_i \over \partial x_i}=0 \quad \hbox{with }\quad \tau _i = \displaystyle{3h_i^2\over 8G}
(58)

In the derivation of Eq.(58) we have neglected the terms involving products Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle h_ih_j}

for Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle h_i \not = h_j}

.

The pressure gradient projections Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \pi _i}

are introduced now as
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \pi _i \equiv r_i - {\partial p \over \partial x_i}
(59)

The final system of governing equations is

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \rho {\partial ^2 u_i\over \partial t^2} - {\partial \sigma _{ij} \over \partial x_i} -{b}_i=0
(60)
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\Delta p\over k}+ {\partial (\Delta u_i) \over \partial x_i} + \sum \limits _{i=1}^{n_d} \tau _i \left({\partial p \over \partial x_i}+\pi _i\right)=0
(61)
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): r_i \equiv {\partial p \over \partial x_i}+\pi _i =0
(62)

In Eqs.(60)–(62) we note the following:

  • The stabilization terms have been neglected in Eq.(60). These terms were useful to derive the stabilized form of Eq.(61). However they are not longer necessary as the convective terms do not appear in the equilibrium equations.
  • The constitutive equation for the pressure (Eq.(61)) has been written in an incremental form. This is more convenient for non linear material behaviour typical of metal forming situations. Here, appropriate elasto-plastic or elasto-viscoplastic constitutive laws relating stresses and strains in incremental or rate forms must be used [1–5,43,44].

6.3 Finite element equations

The weighted residual form of the governing equations can be written as (after integration by parts of the relevant terms)

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \int _{\Omega }\delta u_i \rho \frac{\partial ^2{u_i}}{\partial{t}^2}\,d\Omega + \int _{\Omega }\delta {\varepsilon }_{ij}{\sigma }_{ij}\,d\Omega{-} \int _{\Omega }\delta u_i b_i \,d\Omega - \int _{\Gamma _t}\delta {u}_it_i^p \,d\Gamma _t=0
(63a)
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \int _{\Omega }q \left({\Delta p\over K} + {\partial (\Delta u_i) \over \partial x_i} \right)\,d\Omega + \int _\Omega \left[\sum \limits _{i=1}^{n_d}{\partial q \over \partial x_i} \tau _{i}\left({\partial p \over \partial x_i}+\pi _i \right)\right]d\Omega =0
(63b)
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \int _{\Omega }\delta \pi _i \tau _{i} \left({\partial p \over \partial x_i}+\pi _i \right)d\Omega =0\quad \hbox{no sum in }i
(63c)

where again the Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \tau _i}

are introduced in Eq.(63c) for symmetry reasons.

The finite element discretization of the displacements, the pressure and the pressure gradient projections is written by expressions identical to Eqs.(23) with the nodal variables now being a function of the time Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle t} . Substituting the approximations into eqs.(63) and using the Galerkin form gives the following system of discretized equations

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\boldsymbol M} \ddot{\bar {\boldsymbol u}} + {\boldsymbol g} - {\boldsymbol f}={\boldsymbol 0}
(64a)
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\boldsymbol G}^T \Delta \bar {\boldsymbol u} +{\boldsymbol C} \Delta \bar {\boldsymbol p}+ {\boldsymbol L} \bar {\boldsymbol p} + {\boldsymbol Q}\bar {\boldsymbol \pi }= {\boldsymbol 0}
(64b)
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\boldsymbol Q}^T \bar {\boldsymbol p}+ \bar {\boldsymbol C} \bar {\boldsymbol \pi }= {\boldsymbol 0}
(64c)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \ddot{\bar {\boldsymbol u}}}

is the nodal acceleration vector,
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): M_{ij}=\int _{\Omega ^e} \rho N_i N_j\,d\Omega
(65)

is the mass matrix,

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\boldsymbol g}= \int _\Omega {\boldsymbol B}^T {\boldsymbol \sigma } d\Omega
(66)

is the internal nodal force vector and the rest of matrices and vectors are defined in the Appendix. Note that the expression of Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {\boldsymbol g}}

of eq.(66) is adequate for non linear  analysis.

A four steps semi-explicit time integration algorithm can be derived as follows

step

Step step.step

  • Compute the nodal velocities Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \dot{\bar {\boldsymbol u}}^{n+1/2}}
from Eq.(64a)
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \dot{\bar {\boldsymbol u}}^{n+1/2} = \dot{\bar {\boldsymbol u}}^{n-1/2}+\Delta t {\boldsymbol M}^{-1}_d ({\boldsymbol f}^n - {\boldsymbol g}^n)
(67a)
  • Compute the nodal displacements Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \bar {\boldsymbol u}^{n+1}}


Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \bar {\boldsymbol u}^{n+1} = \bar{\boldsymbol u}^n + \Delta t \dot{\bar {\boldsymbol u}}^{n+1/2}
(67b)
  • Compute the nodal pressures Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \bar {\boldsymbol p}^{n+1}}
from Eq.(64b)
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \bar {\boldsymbol p}^{n+1}=-[ {\boldsymbol C} + {\boldsymbol L}]^{-1} [\Delta t {\boldsymbol G}^T \dot{\bar {\boldsymbol u}}^{n+1/2}- {\boldsymbol C} \bar {\boldsymbol p}^n + {\boldsymbol Q} \bar {\boldsymbol \pi }^n]
(67c)
  • Compute the nodal projected pressure gradients Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \bar {\boldsymbol \pi }^{n+1}}
from Eq.(64c)
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \bar {\boldsymbol \pi }^{n+1}=- \bar {\boldsymbol C}_d^{-1} {\boldsymbol Q}^T \bar {\boldsymbol p}^{n+1}
(67d)

In above, all matrices are evaluated at Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle t^{n+1}} , Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle (\cdot )_d }

 denotes a lumped diagonal matrix and
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\boldsymbol g}^n =\int _{\Omega ^e} [{\boldsymbol B}^T {\boldsymbol \sigma }]^n\,d\Omega
(68)

where the stresses Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {\boldsymbol \sigma }^n}

are obtained by consistent integration of the adequate (non linear) constitutive law.

Note that steps 1, 2 and 4 are fully explicit as a diagonal form of matrices Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {\boldsymbol M}}

and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \bar{\boldsymbol C}}
has been chosen. The solution of step 3 requires invariably the inverse of a Laplacian matrix. This can be an inexpensive process using an iterative equation solution method (e.g. a preconditioned conjugate gradient method).

For the full incompressible case Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle K=\infty }

and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {\boldsymbol C}=0 }
in all above equations.

The critical time step Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \Delta t}

is taken as that of the standard explicit dynamic scheme (see Section 6.5 and [9,37]).

Fully explicit algorithm

A fully explicit four steps algorithm can be obtained by computing Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \bar {\boldsymbol p}^{n+1}}

from step 3 in eq.(67c) as follows
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \bar{\boldsymbol p}^{n+1}= - {\boldsymbol C}_d^{-1}[\Delta t {\boldsymbol G}^T \dot{\bar {\boldsymbol u}}^{n+1/2} - ({\boldsymbol C}_d - {\boldsymbol L}) \bar{\boldsymbol p}^n + {\boldsymbol Q} \bar{\boldsymbol \pi }^n]
(69)

Note that the explicit algorithm is not applicable in the full incompressible limit as the solution of Eq.(69) breaks down for Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle K=\infty }

and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {\boldsymbol C}=0}

. The explicit form can however be used with success in problems where quasi-incompressible regions exist adjacent to standard “compressible” zones. Examples of this kind are shown in the next section. In both cases the semi-explicit and fully explicit schemes gave identical results with important savings in both computer time and memory storage requirements obtained when using the explicit form.

6.4 Thermal coupled effects

The effect of temperature can be easily accounted for by solving the heat transfer equation formulated in a Lagrangian frame as

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \rho c {\partial \phi \over \partial t}-{\partial \over \partial x_j}\left(k{\partial \phi \over \partial x_j}\right)- Q=0
(70)

As usual, the source Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle Q}

is dependent of the mechanical work generated during the forming process.

Note that the convective terms do not enter into Eq.(70) This also eliminates the need to stabilize the numerical solution.

The coupled thermal-mechanical problem requires the computation of the temperature Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \phi }

at each time step using a transient solution scheme [47,48,50].

Remark

Above formulation is similar to that developed by Chiumenti et al. [27,29] and Cervera et al. [28] for analysis of incompressible problems in solid mechanics using a sub-grid scale approach.

6.5 Computation of the intrinsic time parameter for quasi/fully incompressible solids

In solid mechanics applications it is usual to accept that all Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \tau _i}

are identical and constant within each element and given by
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \tau ^{(e)} = {3(h^{(e)})^2\over 8G} \quad \hbox{with}\quad h^{(e)} =[V^{(e)}]^{1/n_d}\, T
(71)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle V^{(e)}}

is the element volume (or the element area for 2D problems). This expression for Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \tau ^{(e)}}
does not take into account the element distorsions along a particular direction during the deformation process.

The correct value of the shear modulus in the expression of Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \tau ^{(e)}}

is another sensitive issue as, obviously, for non linear problems the value of Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle G}
will differ from the elastic modulus. This fact has been identified by Cervera et al. [28] for non linear analysis of incompressible problems using linear triangles.

A useful alternative to compute Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \tau ^{(e)}}

for explicit non linear transient situations is to make use of the value of the speed of sound in an elastic solid, defined by
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): c=\sqrt{{E\over \rho }}
(72)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle E}

is the Young modulus. The stability condition for explicit dynamic computations is given by the Courant condition defined as [8]
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \Delta t^{(e)}\le \Delta t^{(e)}_c ={h^{(e)}\over c}
(73)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \Delta t^{(e)}_c}

is the critical time step for the element.

Accepting that Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle G\simeq {E\over 3}}

for the incompressible case and using Eqs.(72) and (73) (assuming the identity in Eq.(73)) an alternative expression for the element intrinsic time parameter in terms of the critical time step can be found as
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \tau ^{(e)}={[\Delta t^{(e)}_c]^2\over \rho }
(74)

Eq.(74) shows clearly that the intrinsic time parameter varies across the mesh as a function of the critical time step for each element.

7 APPLICATIONS TO CASTING AND BULK METAL FORMING PROCESSES

7.1 Aluminium casting simulation

A numerical simulation of an aluminium casting process is presented as a demonstration of the accuracy of the stabilized formulation. The computations are performed with the finite element code VULCAN where the stabilized FEM presented has been implemented [56].

The analysis simulates the casting process of an aluminium (AlSi7Mg) specimen in a steel (X40CrMoV5) mould. Material behavior of aluminium casting has been modeled by a fully coupled thermo-viscoplastic model, while the steel mould has been modeled by a simpler thermo-elastic model [51]. Geometrical and material data were provided by the foundry RUFFINI. Figure 1 shows the finite element mesh used for the part and the cooling system. The full mesh, including the mould has 380.000 four node tetrahedra. The pouring temperature is 650Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle ^\circ } C. Initial temperature for the mould is obtained through a thermal die-cycling simulation. Figure 2 shows the evolution of the mould temperature after 6 cycles. The cooling system has been kept at 20Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle ^\circ } C. Filling evolution has been simulated as in a pressure die-casting process using the stabilized VOF technique described in [47,48]. Figure 3 shows different time steps of the simulation.

Finite element discretization of the aluminium casting.
Figure 1: Finite element discretization of the aluminium casting.
Temperature die-cycling.
Figure 2: Temperature die-cycling.

The final temperature field obtained after the filling simulation is taken as the initial condition for the solidification and cooling analysis. Temperature and liquid-fraction distributions during solidification are shown in Figures 4 and 5, respectively. The heat transfer coefficient takes into account the air-gap resistance due to the casting shrinkage during the solidification process. Figure 6 shows volumetric and von Mises deviatoric stress distributions in a x-y section. The figures also show the air-gap between the part and the mould, responsible of a non-uniform heat flux at the contact interface.

Other examples of application of the stabilized formulation to casting problems can be found in [47–51].

Draft Samper 805805629-test-MATER.png Draft Samper 805805629-test-MATER-01.png
Draft Samper 805805629-test-MATER-02.png Draft Samper 805805629-test-MATER-03.png
Filling evolution: pressure die-casting simulation.
Figure 3: Filling evolution: pressure die-casting simulation.
Draft Samper 805805629-test-tempe 1.png Draft Samper 805805629-test-tempe 2.png Temperature evolution during cooling phase.
Figure 4: Temperature evolution during cooling phase.
Draft Samper 805805629-test-solid frac.png Draft Samper 805805629-test-solid frac-01.png
Draft Samper 805805629-test-solid frac-02.png Liquid-fraction evolution during phase-change.
Figure 5: Liquid-fraction evolution during phase-change.
Aluminium casting. Stress-trace and von Mises deviatoric stress indicator during  phase-change (plane xy).
Figure 6: Aluminium casting. Stress-trace and von Mises deviatoric stress indicator during phase-change (plane xy).

7.2 Side pressing of a cylinder

A cylinder 100 mm long with a radius of 100 mm is subjected to sidepressing between two plane dies. It is compressed to 100 mm. The material properties are the following: Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle E =217}

GPa, Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \nu =0.3}

, Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \rho =7830}

kg/mFailed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle ^3}

, Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \sigma _0=170}

MPa, Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle H=30}
MPa,  friction coefficient = 0.2. The die velocity is assumed to be 2 m/s. Initial set-up is shown in Figure 7. A quarter of a cylinder was discretized with tri-linear hexahedra and linear tetrahedra meshes.

Figure 8 shows the results obtained using the hexahedral mesh and a standard mixed formulation. The results show the distribution of the effective plastic strain and pressure on the deformed shape. The sensitivity of the FIC results with the expression of the intrinsic time parameter of Eq.(71) was studied by defining Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \tau ^{(e)} =a {[h^{(e)}]^2\over G}}

and solving the problem for different values of Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle a}

. The results obtained with the FIC method are shown in Figs. 9 and 10 for Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle a =0.1}

and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle 0.03}

, respectively. The alternative expression for Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \tau ^{(e)}}

of Eq.(74) has also been studied. The results for this case are shown in Fig. 11. Quite a good agreement can be seen between the FIC solutions and the reference solution with the best results for the effective plastic strain obtained for Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \tau ^{(e)}}
calculated according to Eq. (71) with Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle a=0.03}

. The results for the two alternative formulae for Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \tau ^{(e)}}

are similar, but those obtained using Eq. (71) seem to be slightly better. In any case, the results on the pressure distribution are quite insentitive to the value of Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \tau ^{(e)}}

. This is also confirmed in Fig. 12b, which displays the distribution of the pressure along the line Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle ABCDEA}

defined in Fig. 12a. A small perturbation can be seen at the sharp edges of the deformed body.

All the calculations here have been carried out using fully explicit version of the algorithm which is more efficient than the semi-implicit one with giving practically the same results.

Draft Samper 805805629-test-cyltn3-pi-v5e5-msh.png Sidepressing of a cylinder: (a) initial tetrahedral mesh; (b) initial hexahedral mesh
Figure 7: Sidepressing of a cylinder: (a) initial tetrahedral mesh; (b) initial hexahedral mesh
Draft Samper 805805629-test-cylhc-ep-new.png Sidepressing of a cylinder, mixed formulation,  hexahedral mesh: (a) effective plastic strain; (b) pressure distribution
Figure 8: Sidepressing of a cylinder, mixed formulation, hexahedral mesh: (a) effective plastic strain; (b) pressure distribution
Draft Samper 805805629-test-cylt-fine-fast-fic01-eps.png Sidepressing of a cylinder, FIC algorithm  (α=0.1),  tetrahedra, mesh  of 22186 elements: (a) effective plastic strain; (b) pressure distribution
Figure 9: Sidepressing of a cylinder, FIC algorithm (Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \alpha=0.1

), tetrahedra, mesh of 22186 elements: (a) effective plastic strain; (b) pressure distribution

Draft Samper 805805629-test-cylt-fine-fast-fic-eps.png Sidepressing of a cylinder, FIC algorithm  (a=0.03),  tetrahedra, mesh  of 22186 elements: (a) effective plastic strain; (b) pressure distribution
Figure 10: Sidepressing of a cylinder, FIC algorithm (Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): a=0.03

), tetrahedra, mesh of 22186 elements: (a) effective plastic strain; (b) pressure distribution

Draft Samper 805805629-test-cylt-fine07-fast-v1e12-ab2.png Sidepressing of a cylinder, FIC algorithm  (τ(e) calculated according to Eq. (74)),  tetrahedra, mesh  of 22186 elements: (a) effective plastic strain; (b) pressure distribution
Figure 11: Sidepressing of a cylinder, FIC algorithm (Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \tau ^{(e)}
calculated according to Eq. (74)),  tetrahedra, mesh  of 22186 elements: (a) effective plastic strain; (b) pressure distribution
Draft Samper 805805629-test-cylt-linia.png a) Definition of the  line for comparison of pressure distribution b) Pressure distribution along the line ABCDEA
Figure 12: a) Definition of the line for comparison of pressure distribution b) Pressure distribution along the line Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): ABCDEA

7.3 Backward extrusion

Backward extrusion of a cylinder made of steel 16MNCr5 has been analysed using an axi-symmetric formulation. This is a benchmark example of the finite element program for forming simulation MARC/Autoforge [marc]. The tooling and billet geometry are given in Fig. 13a. Initial material dimensions are the following: length 30 mm and diameter 30 mm. The punch of diameter 20 mm has a prescribed stroke of 28 mm. Material properties are as follows: Young's modulus Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle E=3.24\cdot 10^5}

MPa, Poissson's coefficient Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \nu=0.3}

, material density Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \rho=8120}

kg/mFailed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle ^3}

, yield stress Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \sigma _{Y0}=300}

MPa and hardening modulus Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle H=50}
MPa. Friction between the material and tools is defined by the Coulomb friction coefficient Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \mu=0.1}

.

Draft Samper 805805629-test-bext-geo.png Draft Samper 805805629-test-bext4a-quad-nurbs-sprn9-eps.png
Backward extrusion a) geometry definition. Final deformed shape with effective plastic strain distribution; b)  solution with quadrilaterals and mixed formulation; c)  solution with triangles and the FIC algorithm
Figure 13: Backward extrusion a) geometry definition. Final deformed shape with effective plastic strain distribution; b) solution with quadrilaterals and mixed formulation; c) solution with triangles and the FIC algorithm

The simulation of the backward extrusion process was carried out with a particularization of the FIC formulation of Section 6 for axisymmetric solids. Regeneration of the mesh was performed when element distorsion was excessive. Figures 13b and 13c show the results in the form of the final deformed shape with the distribution of the effective plastic strain obtained using quadrilaterals with a mixed formulation, and using triangles and the FIC algorithm, respectively. The results are in a good agreement with the solution given in [marc]. This example demonstrates the efficiency of the FIC algorithm for simulation of bulk forming processes. Different stages of the forming process are shown in Fig. 14.

Draft Samper 805805629-test-bext4a 002-eps.png Draft Samper 805805629-test-bext4a 014-eps.png
Backward extrusion – deformed shapes with effective plastic strain distribution at different stages of forming. Solution with triangles and the FIC algorithm
Figure 14: Backward extrusion – deformed shapes with effective plastic strain distribution at different stages of forming. Solution with triangles and the FIC algorithm

7.4 Forming of a hose-clamp band

The stabilized formulation was applied to model the manufacturing of a hose-clamp band of steel AISI 409L (Fig. 15). The initial set-up of the tooling and band is shown in Fig. 16. A series of grooves are forged in the band by the roll passing over the band placed on the toothed punch. The band thickness is 0.7 mm and its width 8 mm. Plane strain conditions have been assumed. Material properties are as follows: Young's modulus Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle E=2.1\cdot 10^5}

MPa, Poissson's coefficient Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \nu=0.33}

, material density Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \rho=7800}

kg/mFailed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle ^3}

, the true stress–true strain relationship is given by the power Ludwik–Nadai equation Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \sigma _{Y}=623(0.36822\cdot 10^{-2}+\varepsilon _p)^{0.1362}}

MPa, friction between the material and tools is defined by the Coulomb friction coefficient Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \mu=0.1}

.

A hose clamp
Figure 15: A hose clamp
A hose clamp – initial set-up
Figure 16: A hose clamp – initial set-up
Draft Samper 805805629-test-p1222d 12-eps.png A hose clamp – deformed shapes at different stages of forming with distribution of effective plastic strain
Figure 17: A hose clamp – deformed shapes at different stages of forming with distribution of effective plastic strain

Simulation was carried out using linear triangular elements and the FIC formulation. As in the previous example the meshes where regenerated when element distorsion was excessive. The purpose of the simulation was to check if the expected groove depth and tooth height in the band were obtained. Figures 17a and 17b show the results in the form of the deformed shape with the distribution of the effective plastic strain at different stages of forming. The results are in good agreement with the real process. Figure 18 shows a detail of the deformed shape with finite element discretization and distribution of effective plastic strain. In this figure the obtained dimensions are compared with the required ones shown in brackets. Effects of elastic springback are also clearly seen.

Detail of the deformed shape with finite element discretization and  distribution of effective plastic strain
Figure 18: Detail of the deformed shape with finite element discretization and distribution of effective plastic strain

Remark

The examples presented show that the FIC/FEM formulation is an effective procedure for solving bulk metal forming problems involving full or quasi-incompressible situations. The key advantage of the FIC approach versus more standard mixed FEM formulations is that it provides a natural theoretical framework for equal order finite element interpolations for the velocity and pressure variables, both in the context of implicit and explicit solution schemes. We note the simplicity and effectiveness of the full explicit algorithm as demonstrated in the examples presented. The FIC formulation reproduces also the best feature of the so called stabilized FEM method for incompressible problems, such as the CBS scheme [20,21,25,26,58], the pressure gradient operator method [22] and the subgrid scale method [23] among others.

8 LAGRANGIAN FLOWS. THE PARTICLE FINITE ELEMENT METHOD

8.1 The Particle Finite Element Method (PFEM)

The Lagrangian formulation is an excellent procedure for treating bulk forming processes involving the interaction of fluids and solids using a unified formulation. An important advantage of the Lagrangian formulation is that both the motions of the solid and the fluid are defined in the same frame of reference and modelled with the some governing equations.

The Lagrangian fluid flow equations can be simply obtained by noting that the velocity of the mesh nodes and that of the fluid particles are the same. Hence the the convective terms vanish in the momentum equations, while the rest of the fluid flow equations remain unchanged. The resulting governing equations have an identical form as those of the Stokes flow problem, with the motion of the flow particles being referred now to a Lagrangian coordinate frame.

The FEM algorithms for solving the Lagrangian flow equations are very similar to those presented for incompressible solids in a previous section. Here we present a particular class of Lagrangian formulation called the particle finite element method (PFEM) [38–41]. The PFEM treats the mesh nodes in the fluid and solid domains as dimensionless particles which can freely move an even separate from the main fluid domain representing, for instance, the effect of fluid drops. A finite element mesh connects the nodes defining the discretized domain where the governing equations are solved in the standard FEM fashion.

The quality of the numerical solution depends on the discretization chosen as in the standard FEM. Adaptive mesh refinement techniques can be used to improve the solution in zones where large motions of the fluid or the structure occur.

A typical solution with the PFEM involves the following steps.

  1. Identify the external boundaries for both the fluid and solid domains. This is an essential step as some boundaries (such as the free surface in the fluid) may be severely distorted during the solution process including separation and re-entring of nodes. The Alpha Shape method is used for the boundary definition (see Section 8.4).
  2. Discretize the fluid and solid domains with a finite element mesh. For the mesh generation process we use and extended Delaunay technique [52].
  3. Solve iteratively the coupled Lagrangian equations of motion for the fluid and the solid domains. Compute the relevant state variables in both domains at each time step: velocities, pressure and viscous stresses in the fluid and displacements, stresses and strains in the solid.
  4. The iterative solution scheme chosen is a particularization of the fractional step algorithm of Section 3.1 for Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \alpha =1} . In summary the solution steps are the following.

    Step 3.1. Compute the predicted velocities (viz, Eq.(31))

    Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \tilde{\boldsymbol v}^{n+1,i+1} = \bar{\boldsymbol v}^n - \Delta t {\boldsymbol M}^{-1}_d [{\boldsymbol K} \bar{\boldsymbol v}^{n} - {\boldsymbol G}\bar {\boldsymbol p}^{n} -{\boldsymbol f}^{n+1}]
    (75)

    Step 3.2. Compute Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \bar {\boldsymbol p}^{n+1,i+1}}

    from Eq.(32) for Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \alpha =1}
    

    .

    Step 3.3. Compute explicitely Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \bar{\boldsymbol v}^{n+1,i+1}}

    from Eq.(33) with Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \alpha =1}
    

    .

    Step 3.4. Compute Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \bar{\boldsymbol \pi }^{n+1,i+1}}

    explicitely from Eq.(35).
    

    Step 3.5. Solve for the motion of the solid. This can be performed by integrating the dynamic equations of motion in the solid. Here the algorithm of Section 6 can be used, as it applies to both “compressible” and incompressible materials.

    Step 3.6. Estimate a new position of the mesh nodes in terms of the time increment size as

    Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\boldsymbol x}_j^{n+1,i+1} = {\boldsymbol x}_j^{n}+\bar {\boldsymbol u}_j^{n+1,i+1} \Delta t
    (76)

    It is important to note that all matrices in Steps 3.1–3.5 are evaluated at the last predicted configuration Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \Omega ^{n+1,i}} .

    In steps 3.1–3.6 superindex Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle i}

    denotes the iteration within each time step.
    

    Step 3.7. Check the convergence of the velocity and pressure fields in the fluid and the displacements, strains and stresses in the solid. If convergence is achieved froze the final position of the mesh nodes and move to the next time increment (step 4), otherwise return to step 3.1 for the next iteration.

  5. Go back to step 1 and repeat the solution process for the next time step.

Above algorithm can be found to be analogous to the standard updated lagrangian scheme typically used in non linear solid mechanics problems [8,57].

Despite the motion of the nodes within the iterative process, in general there is no need to regenerate the mesh at each iteration. In the examples presented in the paper the mesh in the fluid domain has been regenerated at each time step. A cheaper alternative is to generate a new mesh only after a prescribed number of converged time steps, or when the nodal displacements induce significant geometrical distorsions in some elements.

The boundary conditions are applied as described in Section 3.1.

In the examples presented in the paper the time increment size has been chosen as

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \Delta t =\min (\Delta t_i ) \quad \hbox{with}\quad \Delta t_i ={\vert {\boldsymbol v}\vert \over h_i^{\min }}
(77)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle h_i^{\min }}

is the distance between node Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle i}
and the closest node in the mesh.

8.2 Treatment of contact between fluid and solid interfaces

The condition of prescribed velocities or pressures at the solid boundaries in the PFEM are applied in strong form to the boundary nodes. These nodes might belong to fixed external boundaries or to moving boundaries linked to the interacting solids. In some problems it is useful to define a layer of nodes adjacent to the external boundary in the fluid where the condition of prescribed velocity is imposed. These nodes typically remain fixed during the solution process. Contact between liquid particles and the solid boundaries is accounted for by the incompressibility condition which naturally prevents the liquid nodes to penetrate into the solid boundaries [38–41].

8.3 Generation of a new mesh

One of the key points for the success of the PFEM is the fast regeneration of a finite element mesh at every time step on the basis of the position of the nodes in the space domain. In our work the mesh is typically generated using the so called extended Delaunay tesselation (EDT) [38,52]. The EDT allows one to generate non standard meshes combining elements of arbitrary polyhedrical shapes (triangles, quadrilaterals and other polygons in 2D and tetrahedra, hexahedra and arbitrary polyhedra in 3D) in a computing time of order Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle n} , where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle n}

is the total number of nodes in the mesh. The Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle C^\circ }
continuous shape functions of the elements are obtained using the so called meshless finite element interpolation (MFEM) [53].

8.4 Identification of boundary surfaces

The PFEM requires the correct definition of the boundary domain. Sometimes, boundary nodes are clearly distinguished from internal nodes. In other cases, the total set of nodes is the only information available and the algorithm must recognize the boundary nodes.

Considering that the nodes follow a variable Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle h(x)}

 distribution, where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle h(x)}
is the minimum distance between two nodes, the following criterion has been used. For each two nodes (three nodes in 3D) a circle (a sphere in 3D) of radius equal to Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \alpha h}
containing the nodes is plotted. All nodes laying on a circle (sphere) with a radius  greater than Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \alpha h

, not containing any other node in the interior are considered as boundary nodes. In practice, Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \alpha }

 is a parameter close to, but greater than one. This criterion coincides with the Alpha Shape concept [39,55].

The method also allows to identify isolated fluid particles (nodes) outside the main fluid domain. These particles are treated as part of the external boundary where the pressure is fixed to the atmospheric value. We note that the “flying particles” are in fact “points” with no mass which motion is followed by integrating the dynamic equations of continuum mechanics for known values of the mass, the initial velocity and acceleration, gravity body forces and a prescribed zero (atmospheric) pressure.

8.5 Temperature coupled effects

Thermal-mechanical problems can be effectively treated with the PFEM. This requires solving for the temperature field at each time step, accounting for the couplings induced by the mechanical equations. The effect of temperature in the mechanical problem is introduced via the constitutive equation in the usual manner.

The form of the heat transfer equation is identical to that of Eq.(70). Recall that the Lagrangian formulation eliminates the convective term in the heat transfer equation and hence the FIC stabilization is here not needed [41].

Remark

The key advantage of the PFEM versus conventional particle method is that it retains the best feature of the FEM to solve problems in continuum mechanics using a variational approach. A standard finite element mesh is used to discretize in space the problem variables and for solving the governing equations precisely as in the standard FEM. The combination of the lagrangian FEM with the Alpha-Shape approach for identification of the domain boundary at each time step and the frequent regeneration of the finite element mesh are the distinct features of the PFEM.

9 APPLICATIONS OF THE PFEM TO METAL FORMING PROCESSES

9.1 Sloshing problem

The sloshing example shown in Figure 19 ilustrates the ability of the PFEM to simulate the flow of liquids within closed cavities with large motions of the free surface. This feature of the PFEM is essential to model the filling of moulds as shown in the next examples. More applications of the PFEM to sloshing problemas are reported in [39].

PFEM results for a large amplitude sloshing problem Draft Samper 805805629-test-Figure19 1.png
(19) PFEM results for a large amplitude sloshing problem
Draft Samper 805805629-test-Figure19 2.png Draft Samper 805805629-test-Figure19 3.png
Draft Samper 805805629-test-Figure19 4.png Draft Samper 805805629-test-Figure19 5.png
a) Filling of a 2D mould with a powder material using the PFEM
Figure 19: a) Filling of a 2D mould with a powder material using the PFEM

A cylinder 100 mm long with a radius of 100 mm is subjected to sidepressing between two plane dies. It is compressed up to 100 mm. The material properties are the following: Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle E =217}

GPa, Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \nu =0.3}

, Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \rho =7830}

kg/mFailed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle ^3}

, Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \sigma _0=170}

MPa, Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle H=30}
MPa, friction coefficient = 0.2. The die velocity is assumed to be 2 m/s. A quarter of a cylinder was discretized using two meshes of 22186 linear tetraedra (4369 nodes) and 1296 hexaedra (1641 nodes) as shown in Figure 20. Figure 21  shows numerical results obtained with the four node hexaedral mesh and an explicit mixed velocity-pressure formulation. Figure 22 shows the pressure and effective plastic strain distributions obtained with the mesh of linear tetrahedra and the full explicit algorithm presented. The maximum and minimum time steps for this solution where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle 0.5 \times 10^{-5}}
and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle 0.25 \times 10^{-5}}

, respectively. Good agreement between the results obtained with the FIC and the mixed formulation is achieved.

Very similar results were obtained with the linear tetrahedral mesh using the more expensive semi-explicit scheme.

Draft Samper 805805629-test-cyltn3-pi-v5e5-msh.png Sidepressing of a cylinder: (a) initial tetrahedral mesh; (b) initial hexahedral mesh
Figure 20: Sidepressing of a cylinder: (a) initial tetrahedral mesh; (b) initial hexahedral mesh
Draft Samper 805805629-test-cylhc-hp-new.png Sidepressing of a cylinder – results obtained using mixed formulation: a) pressure distribution, b) effective plastic distribution
Figure 21: Sidepressing of a cylinder – results obtained using mixed formulation: a) pressure distribution, b) effective plastic distribution
Draft Samper 805805629-test-cylt-fine07-pc-dtrho.png Sidepressing of a cylinder – results obtained with the FIC algorithm for mesh of 22186 elements: a) pressure distribution, b) effective plastic distribution ( τ=\frac(∆t)²ρ )
Figure 22: Sidepressing of a cylinder – results obtained with the FIC algorithm for mesh of 22186 elements: a) pressure distribution, b) effective plastic distribution (Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \tau =\frac{(\Delta t)^2}{\rho }

)

[1] A.J. Chorin, “A numerical solution for solving incompressible viscous flow problems” J. Comp. Phys., 2, 12–26, 1967.

[2] W.R. Briley, S.S. Neerarambam and D.L. Whitfield, “Multigrid algorithm for three-dimensional incompressible high-Reynolds number turbulent flows”, AIAA Journal, 33 (1), 2073–2079, 1995.

Draft Samper 805805629-test-picture1.png Draft Samper 805805629-test-picture50.png
Draft Samper 805805629-test-picture150.png Draft Samper 805805629-test-picture300.png
Draft Samper 805805629-test-picture600.png Draft Samper 805805629-test-picture801.png
Filling of a mould with a casted metal using the PFEM
Figure 23: Filling of a mould with a casted metal using the PFEM
Back to Top

Document information

Published on 01/01/2006

DOI: 10.1016/j.cma.2004.10.018
Licence: CC BY-NC-SA license

Document Score

0

Times cited: 6
Views 31
Recommendations 0

Share this document

claim authorship

Are you one of the authors of this document?