Simplex space-time meshes in thermally coupled two-phase flow simulations of mold filling

TL;DR

This paper introduces a novel 4D simplex space-time discretization method for simulating thermally coupled two-phase flow during mold filling, improving accuracy and efficiency in predicting injection molding defects.

Contribution

It presents a new discretization approach using simplex space-time meshes for adaptive local refinement in complex mold filling simulations.

Findings

Enhanced numerical accuracy near the melt front

Efficient local temporal refinement in 4D space-time grids

Successful application to 2D and 3D benchmark cases

Abstract

The quality of plastic parts produced through injection molding depends on many factors. Especially during the filling stage, defects such as weld lines, burrs, or insufficient filling can occur. Numerical methods need to be employed to improve product quality by means of predicting and simulating the injection molding process. In the current work, a highly viscous incompressible non-isothermal two-phase flow is simulated, which takes place during the cavity filling. The injected melt exhibits a shear-thinning behavior, which is described by the Carreau-WLF model. Besides that, a novel discretization method is used in the context of 4D simplex space-time grids [2]. This method allows for local temporal refinement in the vicinity of, e.g., the evolving front of the melt [10]. Utilizing such an adaptive refinement can lead to locally improved numerical accuracy while maintaining the highest possible computational efficiency in the remaining of the domain. For demonstration purposes, a set of 2D and 3D benchmark cases, that involve the filling of various cavities with a distributor, are presented.

Figures

Tables

Table 1. Table 1: Non-isothermal Carreau-Yasuda-WLF parameters

Carreau-Yasuda-WLF parameter	Value
low shear-rate viscosity, ${\mu }_0$	$\mathrm{50\hspace{0.17em}000}\mathrm{missing}\left[\frac{g}{cm\cdot s}\mathrm{missing}\right]$
high shear-rate viscosity, ${\mu }_{\infty }$	$350\mathrm{missing}\left[\frac{g}{cm\cdot s}\mathrm{missing}\right]$
Power-law coefficient, n	$0.32[-]$
Time constant, $\lambda$	$2.67\left[s\right]$
Scale factor, $a$	$1.915[-]$
WLF parameter, $c_1$	$50[-]$
WLF parameter, $c_2$	$1600[-]$
WLF reference temperature, $T_{ref}$	$390\left[K\right]$

Table 2. Table 2: Material properties

Material property	Value
Density of paste, ${\rho }_{\mathrm{𝑃𝑍𝑇}}$	$4.5\mathrm{missing}\left[\frac{g}{c{m}^3}\mathrm{missing}\right]$
Viscosity of paste, ${\left({\mu }_{\mathrm{𝑒𝑓𝑓}}\right)}_{\mathrm{𝑃𝑍𝑇}}$	See Section 2.2
Heat capacity of paste, ${\left(C_p\right)}_{\mathrm{𝑃𝑍𝑇}}$	$\mathrm{6\hspace{0.17em}200\hspace{0.17em}000}\mathrm{missing}\left[\frac{c{m}^2}{s^2\cdot K}\mathrm{missing}\right]$
Heat conductivity of paste, ${\left(k_0\right)}_{\mathrm{𝑃𝑍𝑇}}$	$\mathrm{61\hspace{0.17em}000}\mathrm{missing}\left[\frac{g\cdot cm}{K}\mathrm{missing}\right]$
Density of air, ${\rho }_{air}$	$0.001\mathrm{missing}\left[\frac{g}{c{m}^3}\mathrm{missing}\right]$
Viscosity of air, ${\left({\mu }_{\mathrm{𝑒𝑓𝑓}}\right)}_{air}$	$0.1\mathrm{missing}\left[\frac{g}{cm\cdot s}\mathrm{missing}\right]$
Heat capacity of air, ${\left(C_p\right)}_{air}$	$\mathrm{10\hspace{0.17em}060\hspace{0.17em}000}\mathrm{missing}\left[\frac{c{m}^2}{s^2\cdot K}\mathrm{missing}\right]$
Heat conductivity of air, $k_{air}$	$2623\mathrm{missing}\left[\frac{g\cdot cm}{K}\mathrm{missing}\right]$

Table 3. Table 3: Typical computational performance of 2D mold filling computations.

	Time	Nodes	Elements	Total time for	Total time for
	Steps	per step	per step	system formation ($s$)	system solution ($s$)
FST (fine)	$1000$	$\mathrm{28\hspace{0.17em}024}$	$\mathrm{27\hspace{0.17em}334}$	$555.09$	$\mathrm{16\hspace{0.17em}186.13}$
FST (coarse)	$200$	$\mathrm{28\hspace{0.17em}024}$	$\mathrm{27\hspace{0.17em}334}$	$76.61$	$4100.53$
SST	$200$	$\sim \mathrm{31\hspace{0.17em}872}$	$\sim \mathrm{104\hspace{0.17em}546}$	$86.42$	$4351.73$

Table 4. Table 4: Typical performance of 3D mold filling computations.

	Time	Nodes	Elements	Total time for	Total time for
	Steps	per step	per step	system formation ($s$)	system solution ($s$)
FST (fine)	$1000$	$\mathrm{70\hspace{0.17em}542}$	$\mathrm{176\hspace{0.17em}398}$	$9298.96$	$7534.43$
FST (coarse)	$200$	$\mathrm{70\hspace{0.17em}542}$	$\mathrm{176\hspace{0.17em}398}$	$1923.78$	$1666.57$
SST	$200$	$\sim \mathrm{100\hspace{0.17em}790}$	$\sim \mathrm{1\hspace{0.17em}325\hspace{0.17em}318}$	$2624.87$	$3307.29$

Equations

$${\rho }_i\left(\frac{\partial \mathbf{u}}{\partial t}+\mathbf{u}\cdot \nabla \mathbf{u}-\mathbf{f}\right)-\nabla \cdot {\mathbf{\sigma }}_i=0,$$

$$\nabla \cdot \mathbf{u}=0,$$

$${\mathbf{\sigma }}_i(\mathbf{u},p)=-p\mathbf{I}+2({\mu }_{eff}{)}_i\mathbf{\epsilon }(\mathbf{u}),$$

$$\mathbf{\epsilon }(\mathbf{u})=\frac{1}{2}(\nabla \mathbf{u}+(\nabla \mathbf{u}{)}^T),$$

$${\rho }_i{\left(C_p\right)}_i(\frac{\partial T}{\partial t}+\mathbf{u}\cdot \nabla T)=k_i{\nabla }^{2}T+{\Phi }_i.$$

$${\Phi }_i=({\mu }_{eff}{)}_i(\nabla \mathbf{u}+(\nabla \mathbf{u}{)}^T):(\nabla \mathbf{u}+(\nabla \mathbf{u}{)}^T).$$

$$\mathbf{n}\cdot [\mathbf{\sigma }{]}_{{\Gamma }_{int}(t)}$$

$$\mathbf{t}\cdot [\mathbf{\sigma }{]}_{{\Gamma }_{int}(t)}$$

$$[\mathbf{u}{]}_{{\Gamma }_{int}(t)}$$

$$\frac{\partial \varphi }{\partial t}+\mathbf{u}\cdot \nabla \varphi =0\mathrm{in}{\Omega }_i(t),\forall t\in (0,t_{total}],$$

$${\mathbf{n}}_s\cdot \mathbf{u}$$

$${\mathbf{t}}_s\cdot \mathbf{\sigma }(\mathbf{u},p)\cdot \mathbf{u}$$

$${\mathbf{b}}_s\cdot \mathbf{\sigma }(\mathbf{u},p)\cdot \mathbf{u}$$

$$\beta ={\beta }_{\infty }\cdot \delta \left(\varphi \right)+{\beta }_0,\delta \left(\varphi \right)=\{\begin{array}{l}\begin{array}{rlr}0,& |\varphi |\le a& \\ \frac{|\varphi |}{a}-1,& a<|\varphi |<2a& \\ 1,& |\varphi |\ge 2a& \end{array}\end{array}$$

$$({\mu }_{eff}{)}_{PZT}={\mu }_{\infty }+\left({\mu }_0-{\mu }_{\infty }\right){\left(1+{\left(\lambda G_f\right)}^a\right)}^{\frac{\left(n-1\right)}{a}},$$

$$({\mu }_{eff}{)}_{PZT}=a_T\left[{\mu }_{\infty }+\left({\mu }_0-{\mu }_{\infty }\right){\left(1+{\left(a_T\lambda G_f\right)}^a\right)}^{\frac{\left(n-1\right)}{a}}\right].$$

$$a_T=exp\left[\frac{c_1\left(T_{ref}-T\right)}{c_2+T-T_{ref}}\right],$$

$$k_{PZT}=\frac{C_p{\mu }_{\infty }}{Pr}.$$

$$Pr=\frac{C_p{\mu }_{\infty }}{(k_0{)}_{PZT}}.$$

$$({\mathcal{S}}_{\mathbf{u}}^h{)}_n$$

$$({\mathcal{V}}_{\mathbf{u}}^h{)}_n$$

$$({\mathcal{S}}_p^h{)}_n$$

$$({\mathcal{S}}_T^h{)}_n$$

$$({\mathcal{V}}_T^h{)}_n$$

$$({\mathcal{S}}_{\varphi }^h{)}_n$$

$$({\mathcal{V}}_{\varphi }^h{)}_n$$

$${\int }_{Q_n}{\mathbf{w}}^h$$

$$+{\int }_{Q_n}{q}^h\nabla \cdot {\mathbf{u}}^{h}dQ+{\int }_{{\Omega }_n}({\mathbf{w}}^h{)}_n^{+}\cdot {\rho }_i(({\mathbf{u}}^h{)}_n^{+}-({\mathbf{u}}^h{)}_n^{-})d\Omega$$

$$+\sum _{e=1}^{(n_{el}{)}_n}{\int }_{Q_n^e}{\tau }_{MOM}\frac{1}{{\rho }_i}\left[{\rho }_i\left(\frac{\partial {\mathbf{w}}^h}{\partial t}+{\mathbf{u}}^h\cdot \nabla {\mathbf{w}}^h\right)-\nabla \cdot {\mathbf{\sigma }}_i({\mathbf{w}}^h,q^h)\right]$$

$$\cdot \left[{\rho }_i\left(\frac{\partial {\mathbf{u}}^h}{\partial t}+{\mathbf{u}}^h\cdot \nabla {\mathbf{u}}^h-\mathbf{f}\right)-\nabla \cdot {\mathbf{\sigma }}_i({\mathbf{u}}^h,p^h)\right]dQ$$

$$+\sum _{e=1}^{(n_{el}{)}_n}{\int }_{Q_n^e}{\tau }_{CONT}\nabla \cdot {\mathbf{w}}^h{\rho }_i\nabla \cdot {\mathbf{u}}^{h}dQ$$

$$={\int }_{(P_{int}{)}_n}{\mathbf{w}}^h\cdot \gamma \kappa \mathbf{n}dP+{\int }_{(P_n{)}_h}{\mathbf{w}}^h\cdot {\mathbf{h}}^{h}dP.$$

$${\int }_{Q_n}{v}^h$$

$$+{\int }_{{\Omega }_n}(v^h{)}_n^{+}((T^h{)}_n^{+}-(T^h{)}_n^{-})d\Omega$$

$$+\sum _{e=1}^{n_{el}}{\int }_{Q_n^e}{\tau }_{TEMP}\left[\frac{\partial v^h}{\partial t}+u^h\cdot \nabla v^h-\frac{k_i}{{\rho }_i{\left(c_v\right)}_i}{\nabla }^2{v}^h\right]$$

$$\left[\frac{\partial T^h}{\partial t}+u^h\cdot \nabla T^h-\frac{k_i}{{\rho }_i{\left(c_v\right)}_i}{\nabla }^2{T}^h-{\Phi }_i(u^h)\right]dQ=0.$$

$${\int }_{Q_n}{\psi }^h$$

$$+\sum _{e=1}^{n_{el}}{\int }_{Q_n^e}{\tau }_{LEV}\left[\frac{\partial {\psi }^h}{\partial t}+u^h\cdot \nabla {\psi }^h\right]\left[\frac{\partial {\varphi }^h}{\partial t}+u^h\cdot \nabla {\varphi }^h\right]dQ=0.$$

$$({\bullet }^h{)}_n^{\pm }$$

$${\int }_{Q_n}\dots dQ$$

$${\int }_{P_n}\dots dP$$

$${\int }_{(P_{int}{)}_n}{\mathbf{w}}^h\cdot \gamma \kappa \mathbf{n}dP$$

$$=-{\int }_{(P_{int}{)}_n}\gamma \underline{\nabla }{\mathbf{id}}_{(P_{int}{)}_n}:\underline{\nabla }{\mathbf{w}}^{h}dP,$$

Full text

Keywords simplex space-time finite elements  $\cdot$ level-set method  $\cdot$ non-isothermal flow  $\cdot$ injection molding  $\cdot$ shear-thinning  $\cdot$ Carreau-WLF

1 Introductions

The thorough understanding of the filling stage in injection molding has been the subject of many experimental and numerical studies since it takes place in a wide variety of manufacturing procedures [11]. During this process, a highly viscous molten material is guided into a cavity and replaces a quiescent gas present inside the mold, meaning that several physical aspects need to be taken into consideration (e.g., the interplay of the two phases due to the discontinuities in their properties, the surface tension, the heat transfer, as well as the wetting, and shear-thinning effects) [6]. A similar range of phenomena governs other production processes, such as die casting [22]. Consequently, numerical methods and material models that are capable of modeling mold filing processes accurately and efficiently are inherently complex.

In this work, a new method is presented, taking advantage of a 4D fully unstructured finite element (FE) formulation [2], founded on simplex space-time discretization. This method is suitable for interface capturing based on the level-set method and designed for two-phase incompressible flows in the context of mold filling. It allows for temporal refinement near areas of interest within the computational domain [10]. With this method, the local computational accuracy is increased while maintaining, as much as possible, the computational efficiency. The current state of the art shows an increasing interest in 4D simplex space-time discretizations, where the time is considered as the fourth dimension. Noteworthy reasons for exploring unstructured space-time discretizations are the space-time mesh adaptivity, utilization of parallelism in 4D, and handling of topological changes. The following references are exploiting some of the advantages as mentioned above. In [7], the discontinuous Petrov Galerkin (DPG) method is used for the solution of a transient acoustic wave system of equations in arbitrary dimensions. In addition, 4D FE discretizations are studied in [25] for the heat equation, the scalar conservation law and the wave equation, which are reformulated as constrained first-order system least-squares (CFOSLS). Regarding the linear algebra aspects, algebraic multigrid (AMG) preconditioned GMRES methods are investigated in [23] for solving linear systems which arise from a space-time FE discretization of the heat equation in 3D and 4D space-time domains.

In the following sections of this publication, the governing equations, the material model and the corresponding solution techniques are presented, followed by a discussion of simplex space-time meshes when used for adaptive temporal refinement. Additionally, a set of numerical examples highlight the performance of the presented methods. Finally, some concluding remarks are given.

2 Theoretical background of mold filling

2.1 Governing equations

During the filling stage of injection molding, a highly viscous incompressible non-isothermal two-phase flow takes place. Such flows are governed by the transient incompressible Navier-Stokes equations, coupled with the heat equation.

The computational domain $\Omega$, on which the governing equations need to be solved, is a subset of $\mathbb{R}^{n_{sd}}$ and contains two immiscible fluids represented by the subdomains $\Omega_{1}(t)$ and $\Omega_{2}(t)$. Here, $\Omega_{1}(t)\cup\Omega_{2}(t)=\Omega$. The boundary of $\Omega$ is given by $\Gamma=\partial{\Omega}$, while the interface between the two immiscible phases, i.e., melt and air, is defined as $\Gamma_{\mathit{int}}(t)=\partial{\Omega_{1}(t)}\cap\partial{\Omega_{2}(t)}$.

At each instant $t\in(0,t_{total}]$, the velocity, $\mathbf{u}(\mathbf{x},t)$, the pressure, $p(\mathbf{x},t)$, and the temperature, $T(\mathbf{x},t)$, in each phase $i$ in subdomain $\Omega_{i}(t)$ are governed by the following equations:

[TABLE]

where $\rho_{i}$ is the density, $(\mu_{\mathit{eff}})_{i}$ is the dynamic effective viscosity, $\boldsymbol{\sigma}_{i}$ denotes the stress tensor, $\boldsymbol{\varepsilon}(\mathbf{u})$ stands for the rate-of-strain tensor and ${\left(C_{p}\right)}_{i}$ implies the isobaric heat capacity, which is assumed to be constant. Moreover, $k_{i}$ is the thermal conductivity. Note that $(\mu_{\mathit{eff}})_{i}$ is a function of the shear rate, defined as $G_{f}=\sqrt{2\ \boldsymbol{\varepsilon}(\mathbf{u}):\boldsymbol{\varepsilon}(\mathbf{u})}$ for viscous incompressible fluids [17]. The indexing $i=1,2$ represents the current phase. Finally, the viscous dissipation $\Phi_{i}$ is defined, as follows:

[TABLE]

The traction acting on the phase interface $\Gamma_{\mathit{int}}(t)$ depends on the surface tension $\gamma(\mathbf{x})$, and the interface curvature $\kappa$. This leads to the following interface conditions $\forall\ t\in(0,t_{total}]$:

[TABLE]

where $[\mathbf{\bullet}]$ is the notation of the usual jump operator across $\Gamma_{\mathit{int}}(t)$. Here, $\mathbf{n}$ and $\mathbf{t}$ denote the outward unit normal and tangent vector on $\Gamma_{\mathit{int}}(t)$, respectively. The connection operator on $\Gamma_{\mathit{int}}(t)$ is represented by $\mathbf{\nabla}^{\Gamma_{\mathit{int}}(t)}$ in Equation (8). Although the surface tension $\gamma(\mathbf{x})$ depends on a local surfactant concentration and temperature, as stated in [21], in the current work $\gamma(\mathbf{x})$ is assumed to be constant. Finally, compatibility between the two phases is enforced by requiring the velocities to be continuous across the interface. For the energy equation, the effects of surface tension, as described in [3], are neglected.

For capturing the evolution of the interface between the two phases, as presented in [14], the level-set transport equation:

[TABLE]

is solved. Here, $\phi$ is a signed-distance function and $\mathbf{u}$ is the velocity field obtained from the Navier-Stokes equations (1) – (2).

After defining the level-set function, the boundary conditions regarding the velocity field on the walls $\Gamma$ of $\Omega$ needs to be described. Instead of using no-slip boundary conditions, which do not allow for the desired wetting during the filling stage of injection molding, the Navier-slip boundary condition is applied to the solid walls. It is a Robin-type boundary condition, as described in [1]:

[TABLE]

with $\mathbf{n}_{s}$ the normal, $\mathbf{t}_{s}$ the tangent, $\mathbf{b}_{s}$ the bi-tangent vector at the boundary and $\beta$ being the Navier-slip coefficient. Following the example of [18], this coefficient is varying as follows:

[TABLE]

with $\beta_{0}$ being the wetting coefficient at the contact line, $\beta_{\infty}$ the far-field wetting coefficient, and $a$ the characteristic length scale around the interface.

2.2 Non-Newtonian fluids

In contrast to Newtonian fluids, the stresses in a non-Newtonian fluid are not linearly proportional to the local rate-of-strain. In the numerical results presented in this work, the molten material is assumed to be non-Newtonian. More specifically, a PZT ceramic paste (paste of lead zirconate titanate) is used and shows shear-thinning behavior. Therefore, the viscosity must be related to the shear rate through an appropriate model.

Following the example of [18], the Carreau-Yasuda model has been used for the melts considered here:

[TABLE]

where $\mu_{0}$ is the viscosity at zero shear rate, $\mu_{\infty}$ the viscosity at infinite shear rate, $\lambda$ the relaxation time, $n$ the power-law index and $a$ the temperature sensitivity scaling factor. As stated in [18], the Carreau-Yasuda material model is well-behaved at the zero shear-rate limit, in contrast to the power-law models. However, values for the low and high shear-rate viscosities must be defined.

Since the model depends on the fluid temperature, the equation needs to be corrected by a temperature shift factor $a_{T}$, as stated in [15]. Thus, the corrected Carreau-Yasuda model, as specified by [18], has the form:

[TABLE]

Although several different formulations for $a_{T}$ exist, the Williams-Landel-Ferry (WLF) formulation is used here:

[TABLE]

with $c_{1}$, $c_{2}$ and $T_{ref}$ being the WLF parameters and the reference temperature, respectively. They are empirically obtained via regression to experimental data. The WLF shift factor allows for the estimation of the viscosity curves for different temperatures than those at which the molten plastic was tested.

The parameter values for the Carreau-Yasuda-WLF model are presented in Table 1 and are those used later in the numerical examples. Correcting the Carreau-Yasuda model with the WLF formulation leads to the shear rate and temperature dependency of the viscosity as presented in Figure 1.

Last but not least, concerning our material model, we need to draw the attention to the thermal conductivity coefficient of the PZT ceramic paste $k_{\mathit{PZT}}$, which is coupled to the dynamic viscosity through the Prandtl number as follows:

[TABLE]

This leads to a better convergence of the heat equation. The Prandtl number is computed based on the parameters, described in Tables 1 and 2:

[TABLE]

2.3 Solution technique

In order to discretize equations (1), (2), (5) and (10), $P1P1$ finite elements are used. The Galerkin/least-squares (GLS) stabilization method is applied for this purpose. In the GLS method, the stabilization term is a least-squares form of the original differential equation, weighted element by element [4]. For the creation of the finite element function spaces, which are used for the space-time discretization, the time interval is divided into subintervals $I_{n}=(t_{n},t_{n+1})$, with $t_{n}$ and $t_{n+1}$ being two consecutive time levels which belong to $0=t_{0}<t_{1}<\cdots<t_{N}=T$. If $\Omega_{n}=\Omega({t_{n}})$, the space-time slab $Q_{n}$ is defined as the region confined by the surfaces $\Omega_{n},\Omega_{n+1}$ and the space-time boundary, described by the spatial boundary $\Gamma(t)$ evolved over $I_{n}$, referred to as $P_{n}$, and contains the advancing space-time front, $(P_{\mathit{int}})_{n}$. The finite element interpolation and weighting function spaces which are set for the variables of the equation system (velocity, pressure, level-set function, and temperature) are:

[TABLE]

for every space-time slab. Here, $H^{1h}\subset H^{1}$ is a finite dimensional Sobolev space consisting of functions which are square-integrable in $Q_{n}$ and have square-integrable first derivatives. The trial function spaces for the velocity, level-set function and temperature, denoted by $(\mathbf{\mathscr{S}}_{\bullet}^{h})_{n}$, must additionally satisfy the Dirichlet boundary conditions. Similar test function spaces $(\mathscr{V}_{\bullet}^{h})_{n}$ are chosen, but it is required that the test functions vanish on the Dirichlet boundary $P_{n}$. The requirements for the pressure trial function space are less restrictive, as no derivatives of the pressure appear and no explicit pressure boundary conditions exist. Therefore, the pressure trial and test functions are chosen from $L^{2h}\subset L^{2}$, being the finite dimensional Hilbert space of square-integrable functions. The interpolation functions in the elements constitute first-order polynomials, which are continuous in space, but discontinuous in time for every region of the domain.

The Navier-Stokes equations (1) and (2), the heat equation (5) and the level-set transport equation (10) have the following space-time discretized forms, respectively, after stabilization:

Given $(\mathbf{u}^{h})_{n}^{-}$, find $\mathbf{u}^{h}\in(\mathscr{S}_{\mathbf{u}}^{h})_{n}$ and $p^{h}\in(\mathscr{S}_{p}^{h})_{n}$ such that $\forall\mathbf{w}^{h}\in(\mathscr{V}_{\mathbf{u}}^{h})_{n},\forall q^{h}\in(\mathscr{V}_{p}^{h})_{n}$:

[TABLE]

Given $(T^{h})_{n}^{-}$, find $T^{h}\in(\mathscr{S}_{T}^{h})_{n}$ such that $\forall v^{h}\in(\mathscr{V}_{T}^{h})_{n}$:

[TABLE]

Given $(\phi^{h})_{n}^{-}$, find $\phi^{h}\in(\mathscr{S}_{\phi}^{h})_{n}$ such that $\forall\psi^{h}\in(\mathscr{V}_{\phi}^{h})_{n}$:

[TABLE]

In the above equations, the notation below is applied:

[TABLE]

The problem is solved sequentially for each space-time slab, starting with $(\mathbf{\bullet}^{h})^{+}_{n}=\mathbf{\bullet}_{0}$ at $t_{0}$. We use a metric stabilization approach. Further information regarding this stabilization method and its parameters $\tau_{MOM}$ and $\tau_{CONT}$ is given in [24] and [16], while the parameters $\tau_{LEV}$ and $\tau_{TEMP}$ are similarly computed. Moreover, the surface tension term in (2.3) can be transformed through the Laplace-Beltrami technique, as suggested in [8, 20, 5], into:

[TABLE]

with $\underline{\Delta}$ being the Laplace-Beltrami operator, $\underline{\nabla}$ the tangential gradient and $\mathbf{id}$ the identity mapping on the space-time evolving interface $(P_{\mathit{int}})_{n}$.

The last term $\int_{(P_{n})_{h}}\mathbf{w}^{h}\cdot\mathbf{h}^{h}\ dP$ in (2.3) is of primary importance for applying the Navier-slip boundary condition to the cavity walls (as stated in [9]), as it is a formulation to apply traction $\mathbf{h}^{h}$ to a boundary of the domain. The term $\mathbf{h}$ corresponds to the imposed traction, whereas the superscript $h$ stands for discretized quantities.

2.4 Coupling between Navier-Stokes, level-set and heat equations

In this section, it is worth pointing out the mutual interaction between the fluid velocity, level-set, and temperature fields. On the one hand, the quantities that depend on the level-set field are the surface tension coefficient, the density, the viscosity, the isobaric heat capacity, and the thermal conductivity. On top of that, the viscosity of the melt is also a function of the temperature field, since the melt is considered to behave as a shear-thinning fluid (cf. Section 2.2). On the other hand, the level-set and temperature fields are advected by the fluid velocity. For this paper, a partitioned approach for the coupling has been chosen, meaning that the flow field is computed using a fixed melt-air interface. The temperature distribution is then calculated utilizing the already determined flow field and assuming a fixed melt-air front. Finally, the velocity field advects the melt-air interface. If this procedure is performed once per time step, then we refer to the coupling approach as weak. If the previously described process is repeated until all fields are in equilibrium and the desired convergence is reached, as described schematically in Figure 2, then a strong coupling is obtained. The coupling is repeated at every time step, $n$, until the total number of time steps, $n_{ts}$, is reached. If a strongly coupled partitioned approach is used, the level-set, $\phi^{\theta}$, and temperature field, $T^{\theta}$, is a relaxation between the level-set and temperature field from the previous coupling iteration, $\phi^{j-1}$ and $T^{j-1}$, and those from the current iteration, $\phi^{j}$ and $T^{j}$, respectively.

On the one hand, a strongly coupled partitioned approach is more precise, providing more stable results in comparison to a weakly coupled partitioned approach. On the other hand, it requires more computation time, as stated in [20, 19].

3 Simplex space-time finite elements

An algorithm for the construction of simplex space-time meshes is proposed by [2]. If $n_{sd}$ is the number of spatial dimensions, an $\left(n_{sd}+1\right)$-simplex is a space-time element, created after extruding the spatial simplex element in time and subdividing it into simplices, and contains $\left(n_{sd}+2\right)$ vertices. Furthermore, we need to point out that an $n_{sd}$-dimensional spatial simplex extruded in time can be subdivided into $n_{sd}+1$ simplicial elements.

The simplex space-time elements allow for local temporal refinement, by adding nodes between the top and the bottom of the space-time slab. In the case of two-phase flow problems, an algorithm that applies temporal refinement in the vicinity of the evolving interface is described in [10].

For setting the boundary conditions (11) – (13), the computation of the local normal–tangent–bi-tangent coordinate system at the solid boundaries of the mold must be performed, as described in Section 2.1. This computation is based on a generalized method for calculating the cross-product in higher dimensions, as stated in [12, 13].

4 Numerical Results

4.1 2D cavity with a distributor

As a first test case, we simulate the filling stage of a 2D geometry consisting of a pipe and a cavity, connected through a distributor. This benchmark case was also simulated in [18] with the software GOMA, developed by Sandia National Laboratories.

Figure 3 illustrates the computational domain. The spatial discretization of the domain consists of $27\,334\text{\,}$ triangular elements. The time-slab size varies between $\Delta t=$0.001\text{,}\mathrm{s}$-$0.005\text{,}\mathrm{s}. Navier-slip boundary condition is assumed for the velocity degrees of freedom on the walls, except for the inflow (leftmost vertical boundary) and outflow (rightmost vertical boundary). The molten PZT ceramic paste enters the cavity with a parabolic inflow velocity profile ($u=168\frac{y}{0.2}\left(1-\frac{y}{0.2}\right)$\text{\,}\frac{cm}{s}) and displaces the air, which is initially quiescent. Traction-free boundary condition is used at the outflow boundary. With respect to the temperature degree of freedom, homogeneous Neumann boundary conditions are imposed on the mold walls, except for the inflow one, where the injection temperature is equal to $405\text{\,}\mathrm{K}$. Although we consider non-isothermal conditions, the phase-change effects are disregarded. The material properties correspond to those of [18] and are shown in Table 2. The gravitational acceleration is neglected, but the surface tension coefficient remains constant and equal to $\gamma(\mathbf{x})=$42.4\text{,}\frac{g}{s^{2}}$$.

The PZT ceramic paste fills the cavity for $0.1\text{\,}\mathrm{s}$. We use three different temporal discretizations here. The first one is a flat space-time (FST) discretization, consisting of $200\text{\,}$ time slabs ($\Delta t=$0.0005\text{,}\mathrm{s}$$), in contrast to an FST discretization of $1000\text{\,}$ time slabs. Furthermore, we apply adaptive temporal refinement, as described in [10], which leads to a hybrid simplex space-time (SST) discretization, with each of the $200\text{\,}$ time slabs being $0.0005\text{\,}\mathrm{s}$ thick and discretized differently with $1\text{\,}$ to $5\text{\,}$ elements in the time direction (see Figure 4). Note that the temporal accuracy is increased only close to the propagating interface, as shown in Figure 4.

Figure 5 illustrates the front position of the molten material at various time instances. The results obtained with the fine and coarse FST discretization are then compared with those computed with the SST discretization while applying temporal refinement. As we can realize from Figure 6, the wetting, and the viscosity profile are better resolved by the SST and the finer FST discretization compared to the coarse FST discretization because they provide higher temporal accuracy.

Furthermore, as Figure 7 indicates, the simulation with the FST discretization leads to a viscosity distribution qualitatively similar to that presented in [18]. Inside the pipe, the viscosity remains high in the center and low near the walls, whereas inside the cavity, the viscosity increases progressively. This especially holds close to the interface, in both simulations. It is worth mentioning that a linearized Blake model in combination with a Navier-slip boundary condition regarding the velocity degrees of freedom is used in [18]. Furthermore, in contrast to the results presented here, a heat flux boundary condition was applied to the walls of the steel mold. Another difference is that, in this work, a finer spatial discretization for the computations is used. All the factors mentioned above could explain the quantitative differences, which exist between the two results.

In Table 3, we present the performance of the simulations, using three different temporal discretizations, as described above. All simulations were performed in parallel (MPI paralleliztion), utilizing $64\text{\,}$ cores on the RWTH Aachen University IT Center cluster.

We perform Newton-Raphson iterations at every time step, because of the nonlinearity of the Navier-Stokes equations and the heat equation, due to the viscous dissipation term. Strong coupling iterations are also executed due to the mutual dependence between the Navier-Stokes, the level-set, and the heat equation. A GMRES solver is employed for solving the resulting linear systems of equations, in combination with an ILUT factorization. To summarize, the use of SST discretization combined with adaptive temporal refinement leads to a similar resolution of the space-time evolving interface and the viscosity distribution as the one achieved by the fine FST discretization, while keeping the total time for forming and solving the system around $70\%$ lower than in the case of the fine FST simulation (cf. Fig. 6).

4.2 3D cavity with a distributor

As a second test case, we compute the filling process of a 3D coat hanger distributor and a cavity. The inflow is at the small end of the pipe, whereas the outflow is at the back of the cavity on the opposite side. Thanks to the geometry’s symmetry relative to the center plane, only half of the geometry is simulated, saving time and resources. The same test case was simulated in [18] with the software ARIA.

In this benchmark, we observe the temperature and viscosity profile, as well as the shape of the evolving interface. Although a Newtonian melt without any temperature gradients is assumed in [18], shear-thinning effects are taken into consideration in our first computation and represented with the Carreau-WLF model. The domain is illustrated in Figure 8. For our first simulation, we use a very fine spatial discretization, consisting of $7\,119\,241\text{\,}$ tetrahedral elements and the time-slab thickness is $\Delta t=$0.0005\text{,}\mathrm{s}. Navier-slip boundary condition is assumed on the walls for the velocity degrees of freedom, except for the inflow, outflow and symmetry plane. The Blake wetting condition is not taken into consideration here. A constant parabolic velocity profile, which leads to an average velocity of $190\text{\,}\frac{cm}{s}$, is imposed at the inflow boundary, whereas traction-free boundary conditions are used at the outflow boundary. A slip boundary condition is used for the symmetry plane. Concerning the temperature degree of freedom, we impose a Dirichlet boundary condition on the mold walls $T=$298\text{\,}\mathrm{K}, the injection temperature is equal to $405\text{\,}\mathrm{K}$ at the inflow, and we assume homogeneous Neumann boundary conditions at the symmetry plane. Although we consider a non-isothermal condition, the phase-change effects are disregarded. The material properties are given in Table 2. The gravitational acceleration is equal to $f_{y}=-g=$-981\text{,}\frac{cm}{s^{2}}$$ and the surface tension coefficient $\gamma$ is set to $42.4\text{\,}\frac{g}{s^{2}}$.

In Figure 9, we show the velocity, pressure and temperature fields on a slice plane inside the cavity close to the inflow from the distributor at $t=$0.125\text{,}\mathrm{s}$$. As expected, the highest velocity occurs at the inflow of the cavity, whereas at the walls the velocity remains close to zero due to the Navier-slip boundary condition. The temperature distribution exhibits similar behavior since we assume adiabatic walls. We observe a gradient of pressure at the position, where the distributor opens to the mold, with the highest pressure happening at the center of the channel.

Figure 10 shows the shape of the front between the molten PZT and air as well as the viscosity field at $t=$0.09\text{,}\mathrm{s}. As seen, the viscosity values close to the walls of the mold are over $4000\text{\,}\frac{g}{cm\cdot s}$. This observation is of great significance because the viscosity values on the cavity boundaries influence the wetting behavior of the melt the most, as already noted in [[18](#bib.bib18)]. As a consequence, we decided to perfom the same filling simulation, considering this time isothermal conditions and the PZT ceramic paste as Newtonian fluid with constant viscosity of value $4000\text{\,}\frac{g}{cm\cdot s}$. For maintaining a good convergence rate, we used the following properties for air: $\rho_{air}=$0.0045\text{\,}\frac{g}{cm^{3}} and $(\mu_{\mathit{eff}})_{air}=$4.0\text{,}\frac{g}{cm\cdot s}$$.

For the later simulations, the spatial discretization is much coarser and consists of $176\,398\text{\,}$ tetrahedral elements. We made use of three different temporal discretizations, which are the same as in the 2D benchmark case. In Figure 11, we compare the molten material interface at two different time instances, as obtained with our in-house solver utilizing three different temporal discretizations, and the one shown in reference data [18]. In Table 4, we present the timings of our 3D simulations, employing three different temporal discretizations. All three simulations were computed again in parallel (MPI parallelization) utilizing $64\text{\,}$ cores on the RWTH Aachen University IT Center cluster. Similar to the 2D filling case, the use of SST discretization combined with adaptive temporal refinement reduces the total time for forming and solving the system by roughly $55\%$, compared to the case of the fine FST simulation.

5 Conclusions

In this paper, the non-isothermal two-phase flow of a highly viscous fluid has been computed and shows shear-thinning effects, during the filling process of injection molding. The Carreau-WLF material model described the behavior of the melt. For achieving a better wetting on the mold walls, the Navier-slip boundary condition has been used.

Besides, a novel discretization approach has been presented, which allows arbitrary temporal refinement of the space-time slabs in the vicinity of the evolving front of the highly viscous molten material during injection molding. Future work includes the combination of arbitrary temporal refinement with arbitrary spatial refinement close to the propagating interface. The refinement criterion should be based on an appropriate a-posteriori error estimate of local or global quantities of interest, such as the front curvature, pressure jumps close to the interface, material discontinuities and the gradient in the viscosity field.

Acknowledgment

The authors gratefully acknowledge the support of the German Research Foundation (DFG) under program SFB 1120 “Precision Melt Engineering". The computations were conducted on computing clusters provided by the RWTH Aachen University IT Center and by the Jülich Aachen Research Alliance (JARA). Furthermore, we would like to thank the student Efstratios Moskofidis for his contribution by simulating different benchmark cases and our colleague, Max von Danwitz for the fruitful discussions.