Non-newtonian laminar 2D swirl flow design by the topology optimization method

Abstract

The performance of fluid devices, such as channels, valves, nozzles, and pumps, may be improved by designing them through the topology optimization method. There are various fluid flow problems that can be elaborated in order to design fluid devices and among them there is a specific type which comprises axisymmetric flow with a rotation (swirl flow) around an axis. This specific type of problem allows the simplification of the computationally more expensive 3D fluid flow model to a computationally less expensive 2D swirl flow model. The topology optimization method applied to a Newtonian fluid in 2D swirl flow has already been analyzed before, however not all fluids feature Newtonian (linear) properties, and can exhibit non-Newtonian (nonlinear) effects, such as shear-thinning, which means that the fluid should feature a higher viscosity when under lower shear stresses. Some fluids that exhibit such behavior are, for example, blood, activated sludge, and ketchup. In this work, the effect of a non-Newtonian fluid flow is considered for the design of 2D swirl flow devices by using the topology optimization method. The non-Newtonian fluid is modeled by the Carreau-Yasuda model, which is known to be able to accurately predict velocity distributions for blood flow. The design comprises the minimization of the relative energy dissipation considering the viscous, porous, and inertial effects, and is solved by using the finite element method. The traditional pseudo-density material model for topology optimization is adopted with a nodal design variable. A penalization scheme is introduced for 2D swirl flow in order to enforce the low shear stress behavior of the non-Newtonian viscosity inside the modeled solid material. The optimization is performed with IPOPT (Interior Point Optimization algorithm). Numerical examples are presented for some 2D swirl flow problems, comparing the non-Newtonian with the Newtonian fluid designs.

Appendices

Appendix 1. Continuous adjoint model for the non-Newtonian 2D swirl flow

The continuous adjoint model for the 2D swirl flow problem is derived as follows. The adjoint (dual) equations for Navier-Stokes flow have already been deduced in Brandenburg et al. (2009). However, in this Appendix, they are particularized for the 2D swirl flow model in a rotating reference frame, and an approach for dealing with the non-Newtonian viscosity is suggested. In the following development, 2D coordinates are considered in the equations, the domain is given in cylindrical coordinates (in which the differential volume and area are given by, respectively, 2πrdΩ and 2πrdΓ), axisymmetry is considered ($\frac {\partial (\ )}{\partial \theta } = 0$), and the differential operators correspond to their cylindrical coordinate system versions (Lai et al. 2009).

The adjoint equation is first presented in Section 4.5 and is based on the Lagrangian function of the optimization problem, which is given by

$$ L((\textbf{\textit{v}}, p), \alpha, (\boldsymbol{\lambda}_{v}, \lambda_{p})) = J((\textbf{\textit{v}}, p), \alpha) - F((\textbf{\textit{v}}, p), \alpha, (\boldsymbol{\lambda}_{v}, \lambda_{p})) $$

(21)

where (v, p) are the state (primal) variables (velocity and pressure), α is the design variable, (λ_v, λ_p) are the adjoint (dual) variables (adjoint velocity and adjoint pressure) (that is, the adjoint variable presented in Section 4.5 separated in its components: λ_J = (λ_v, λ_p)), J((v, p), α) = Φ_rel((v, p), α) is the objective function (relative energy dissipation), and F((v, p), α, (λ_v, λ_p)) is given in (12) (i.e., (10) and (11) without the division by 2π, and with the test functions w_v and w_p replaced by the adjoint variables λ_v and λ_p, respectively).

Then, in order to obtain the weak form of the adjoint equation (F_λ = 0), the equations that compose (21) need to be derived in function of the state variables (v, p), as shown in Section 4.5. This is given by the directional derivative of the Lagrangian function (21), with respect to the state variables (v and p) and in the directions given by $\textit {\textbf {w}}_{\lambda ,v} =\left [\begin {array}{ll}{\textit {w}}_{\lambda ,v,r}\\{\textit {w}}_{\lambda ,v,\theta }\\ {\textit {w}}_{\lambda ,v,z} \end {array}\right ]$ and w_{λ, p} (test functions for the adjoint equations), respectively for each state variable:

$$ F_{\lambda} = L((\textbf{\textit{v}}, p); \underline{(\textbf{\textit{w}}_{\lambda,v}, {\textit{w}}_{\lambda,p})}) = L(\textbf{\textit{v}}; \underline{\textbf{\textit{w}}_{\lambda,v}}) + L(p; \underline{{\textit{w}}_{\lambda,p}}) $$

(22)

where $L(\textit {\textbf {v}}; \underline {{\textit {\textbf {w}}}_{\lambda ,v}})$ is the directional derivative of L with respect to v in the direction of w_{λ, v}, and $L(p; \underline {\textit {w}_{\lambda ,p}})$ is the directional derivative of L with respect to p in the direction of w_{λ, p}.

The dependency of the non-Newtonian viscosity (μ) with respect to the state variables (v, p) can be separated by applying the chain rule. In this case,

$$ \begin{array}{@{}rcl@{}} F_{\lambda} &=& L((\textbf{\textit{v}}, p), \mu(\textbf{\textit{v}}, p); \underline{(\textbf{\textit{w}}_{\lambda,v}, {\textit{w}}_{\lambda,p})}) \\ &=&L((\textbf{\textit{v}}, p); \underline{(\textbf{\textit{w}}_{\lambda,v}, {\textit{w}}_{\lambda,p})})\\ &&+ \frac{\partial L(\mu)}{\partial \mu} \mu((\textbf{\textit{v}}, p); \underline{(\textbf{\textit{w}}_{\lambda,v}, {\textit{w}}_{\lambda,p})}) \end{array} $$

(23)

where $\mu ((\textit {\textbf {v}}, p); \underline {(\textit {\textbf {w}}_{\lambda ,v}, {\textit {w}}_{\lambda ,p})})$ is the directional derivative of μ with respect to (v, p) in the direction of (w_{λ, v}, w_{λ, p}). In the case of the non-Newtonian viscosity given by (7), which does not depend on p, $\mu ((\textit {\textbf {v}}, p); \underline {(\textit {\textbf {w}}_{\lambda ,v}, \textit {w}_{\lambda ,p})}) = \mu (\textit {\textbf {v}}; \underline {\textit {\textbf {w}}_{\lambda ,v}})$. Note that, if $\mu = \mu _{\infty }$ (Newtonian viscosity), $\mu ((\textit {\textbf {v}}, p); \underline {(\textit {\textbf {w}}_{\lambda ,v}, \textit {w}_{\lambda ,p})}) = 0$.

The first term of the weak form of the adjoint equation (23) ($L((\textit {\textbf {v}}, p); \underline {(\textit {\textbf {w}}_{\lambda ,v}, {\textit {w}}_{\lambda ,p}}))$) becomes, after dividing the forward equations by 2π:

(24)

where T(w_{λ, v}, w_{λ, p}) is the stress tensor (3) calculated by substituting v and p by w_{λ, v} and w_{λ, p}, respectively; and the symbols “(+)” and “(–)” serve to indicate the signal that is already considered in the equation and that is multiplying each adjoint form.

The second term of the weak form of the adjoint equation (23) ($\frac {\partial L(\mu )}{\partial \mu } \mu (\textit {\textbf {v}}; \underline {\textit {\textbf {w}}_{\lambda ,v}})$) can be calculated from (11) as, after dividing the forward equations by 2π:

(25)

where $\frac {\partial \textit {\textbf {T}}(\textit {\textbf {v}}, p)}{\partial \mu } = \nabla \textit {\textbf {v}} + \nabla \textit {\textbf {v}}^{T}$. Since it may be difficult and laborious to calculate $\mu (\textit {\textbf {v}}; \underline {\textit {\textbf {w}}_{\lambda ,v}})$ (sensitivity of the non-Newtonian viscosity) from (7) or (14 (or any other non-Newtonian fluid with a more complex constitutive equation), this term may be calculated by automatic differentiation, such as the algorithm used in the FEniCS platform.

From (23), it would be possible to derive the strong form of the adjoint equation (by applying Gauss’s Theorem of Divergence). However, it would require the sensitivity of the non-Newtonian viscosity to be analytically evaluated. Since the finite element method only requires the weak form of the adjoint equation, this step does not need to be performed.

The Dirichlet boundary conditions from (9) assume constant velocity values, which means that their corresponding adjoint boundary conditions are homogeneous (i.e., equal to zero). By also including the Neumann boundary condition, the adjoint boundary conditions become, from (9):

$$ \begin{array}{@{}rcl@{}} \boldsymbol{\lambda}_{v} = \textbf{{0}} ~~~~~~~~~~~~~~~&& \text{ on } \mathrm{\Gamma}_{\text{in}}\\ \boldsymbol{\lambda}_{v} = \textbf{{0}} ~~~~~~~~~~~~~~~&& \text{ on } \mathrm{\Gamma}_{{{\text{wall}}}}\\ {\lambda_{v,r}} = {0} \text{ and } \frac{\partial (\ )}{\partial r } = {0} && \text{ on } \mathrm{\Gamma}_{\text{sym}}\\ \textbf{\textit{T}}({\textbf{\textit{w}}}_{\lambda, v}, {{\textit{w}}}_{\lambda, p}) \bullet \textbf{\textit{n}} = \textbf{{0}} && \text{ on } \mathrm{\Gamma}_{\text{out}} \end{array} $$

(26)

where λ_{v, r} is the radial component of the adjoint velocity (λ_v = (λ_{v, r}, λ_{v, 𝜃}, λ_{v, z})). It can also be mentioned that, in (25), the term which relies on $\left (\frac {\partial \textit {\textbf {T}}(\textit {\textbf {v}}, p)}{\partial \mu }{\bullet }\boldsymbol {\lambda }_{v}\right ){\bullet }\textit {\textbf {n}}$ is zero on Γ_out, because T(v, p)∙n = 0 on Γ_out (9) ($\left (\frac {\partial \textit {\textbf {T}}(\textit {\textbf {v}}, p)}{\partial \mu }{\bullet }\boldsymbol {\lambda }_{v}\right ){\bullet }\textit {\textbf {n}}$$= \left (\frac {\partial \textit {\textbf {T}}(\textit {\textbf {v}}, p)}{\partial \mu }{\bullet }\textit {\textbf {n}}\right ){\bullet }\boldsymbol {\lambda }_{v} = \left (\frac {\partial [\textit {\textbf {T}}(\textit {\textbf {v}}, p){\bullet }\textit {\textbf {n}}]}{\partial \mu }\right ){\bullet }\boldsymbol {\lambda }_{v} = \textbf {{0}}$, since T is symmetric and n does not depend on μ).

Appendix 2. Simulation of the effect of the non-Newtonian penalization

In order to check the effect of the non-Newtonian penalization, a test example of a channel with an obstacle in the middle is simulated (see Fig. 31). The obstacle is modeled by using a rotating material model (ω_mat = ω₀) located in the middle of the channel, the flow rate is 0.5 L/min, the rotation is 20 rpm, and the dimensions are given by: H = 15.0 mm, R = 10.0 mm, r_o = 2.5 mm, and h_o = 2.5 mm. The mesh is the same as the one shown in Fig. 25. The outlet boundary condition is a weak imposition of zero pressure on the outlet, imposing the radial and axial components of the velocity (v_r, v_z) to be perpendicular to the outlet section (i.e., with zero tangential component) (Dirichlet boundary condition), and imposing zero normal stress on the interface (n∙Tn = 0) (Neumann boundary condition). The demonstration of this boundary condition is shown in Alonso et al. (2018) for 2D swirl flow, and in Barth and Carey (2007) for 2D flow.

Fig. 31

Computational domain for the obstacle simulation

When the material model is blocking the fluid flow (i.e., with a sufficiently high value in $\kappa _{\max \limits }$, chosen in this case as $\kappa _{\max \limits } = 2.5 \times 10^{8} \mu _{\infty }$), there is little difference between using only the inverse permeability (“κ(α)”), or using both inverse permeability and non-Newtonian penalization (“κ(α), μ(α)”). This is shown in Fig. 32a. Figure 32b shows the difference fractions for the modeled obstacles, which correspond to the absolute difference of the variable in the cut and modeled obstacles divided by the range of the variable: $f_{x} = \frac {|x_{\text {material}} - x_{\text {cut}}|}{\max \limits (x_{\text {cut}}) - \min \limits (x_{\text {cut}})}$, where x is the variable being considered: relative tangential velocity (v_𝜃), pressure (p) or non-Newtonian viscosity (μ). For the magnitude of the radial-axial velocity ($v_{rz,\text {mag}} = \|\textit {\textbf {v}}_{rz}\| = \|(v_{r}, v_{z})\| = \sqrt {{v_{r}^{2}} + {v_{z}^{2}}}$), the definition is changed to $f_{{v}_{rz,\text {mag}}} = \frac {\|\textit {\textbf {v}}_{rz,\text {material}} - \textit {\textbf {v}}_{rz,\text {cut}}\|}{\max \limits (\|\textit {\textbf {v}}_{rz,\text {cut}}\|) - \min \limits (\|\textit {\textbf {v}}_{rz,\text {cut}}\|)}$. From Fig. 32b: $f_{{v}_{rz,\text {mag}}}$ is mostly small, but features a 0.33 peak near the edge for “κ(α), μ(α)”; $f_{{v}_{\theta }}$ is even smaller, with a maximum value of 0.091; f_p is higher near the edge of the obstacle, reaching a relatively high peak of 0.5 in “κ(α), μ(α)”; and f_μ is higher around the obstacle for “κ(α), μ(α)”. As can be noticed, the highest values of the difference fractions are concentrated near the obstacle / edge of the obstacle. These differences, and, more specifically, the higher values of f_μ on the obstacle, are mainly due to the nodal interpolation used for the design variable (α), which does not exactly match the effect of the “cut obstacle,” such that it “slightly softens” the effect of the edge because of the linear interpolation, and “forces,” when considering the non-Newtonian penalization (“κ(α), μ(α)”), the nodal values of the non-Newtonian viscosity (μ) to the maximum dynamic viscosity value (μ₀).

Fig. 32

Plots of cut and modeled obstacles ($\kappa _{\max \limits } = 2.5 \times 10^{8} \mu _{\infty }$, q = 1000, $\kappa _{\min \limits } = 0$)

The difference between simulation results is more apparent for lower values of $\kappa _{\max \limits }$, which may occur during topology optimization due to the interpolation of the material model (“gray values”). Therefore, $\kappa _{\max \limits }$ is reduced to $\kappa _{\max \limits } = 2.5 \times 10^{6} \mu _{\infty }$ in Fig. 33. From Fig. 33, when including the non-Newtonian penalization (“κ(α), μ(α)”), the radial-axial velocity ((v_r, v_z)) is more reduced inside the modeled obstacle than in the case without the non-Newtonian penalization (“κ(α)”). This is due to the higher and uniform non-Newtonian viscosity inside the modeled obstacle. The relative tangential velocity (v_𝜃) does not show significant differences between both modeled obstacles.

Fig. 33

Plots of cut and modeled obstacles ($\kappa _{\max \limits } = 2.5 \times 10^{6} \mu _{\infty }$, q = 1000, $\kappa _{\min \limits } = 0$)

Therefore, by imposing a uniform non-Newtonian viscosity inside the solid material of the obstacle, the non-Newtonian penalization seems to show a more noticeable effect in the radial-axial velocity ((v_r, v_z)) during topology optimization rather than in its end.