On the virtual element method for topology optimization of non-Newtonian fluid-flow problems

Abstract

This paper presents some virtual element method (VEM) applications for topology optimization of non-Newtonian fluid-flow problems in arbitrary two-dimensional domains. The objective is to design an optimal layout for the incompressible non-Newtonian fluid flow, governed by the Navier–Stokes–Brinkman equations, to minimize the viscous drag. The porosity approach is used in the topology optimization formulation. The VEM is used to solve the governing boundary value problem. The key feature distinguishing the VEM from the classical finite element method is that the local basis functions in the VEM are only implicitly known. Instead, the VEM uses local projection operators to describe each element’s rigid body motion and constant strain components. Therefore, the VEM can handle meshes with arbitrarily shaped elements. Several numerical examples are provided to demonstrate the efficacy and efficiency of the VEM for the topology optimization of fluid-flow problems. A MATLAB code for reproducing the results provided in this paper is freely available at https://github.com/mampueros/VEM_TopOpt_FluidFlow.

Appendix A: virtual element projections

1.1 Virtual element projection ${\Pi }_E^0 \nabla v$

Based on the works of Gain et al. [24], Beirão da Veiga et al. [50] and Brenner et al. [59], we consider the lower order element $(i.e., k=1,n_{p_{k}}=3$ and $n_{p_{k-1}}=2)$ in which the set of basis functions for $P_k (E)$, (linear polynomial space): $m_\alpha ^{(k)},\alpha =1,\ldots , n_{p_{k}}$ is defined as

$$\begin{aligned} \begin{aligned} m_1^{(1)}=1, \ m_2^{(1)}= \displaystyle \frac{ x-x_C}{ h_E}, \ m_3^{(1)}= \displaystyle \frac{y-y_C}{h_E}, \end{aligned} \end{aligned}$$

(51)

the gradients of $P_k (E)$ are

$$\begin{aligned} \begin{aligned} \nabla m_1^{(1)}&=\begin{bmatrix} 0\\ 0\\ \end{bmatrix}, \ \nabla m_2^{(1)}&=\begin{bmatrix} 1\\ 0\\ \end{bmatrix}, \ \nabla m_3^{(1)}&=\begin{bmatrix} 0\\ 1\\ \end{bmatrix} , \end{aligned} \end{aligned}$$

(52)

and the basis functions $\varvec{m}_\alpha ^{(k-1)}$, for the two-dimensional vector polynomial space $[P_{k-1}(E)]^2$ with $\alpha =1,\dots ,n_{p_{k-1}}$ are defined as

$$\begin{aligned} \begin{aligned} \varvec{m}_1^{(0)}&=\begin{bmatrix} 1\\ 0\\ \end{bmatrix} \text {and} \ \varvec{m}_2^{(0)}&=\begin{bmatrix} 0\\ 1\\ \end{bmatrix} , \end{aligned} \end{aligned}$$

(53)

where $x_c$ and $y_c$ are the coordinates of the centroid of the element E, |E| is the element area, and $h_E=|E|^{1/2}$ is the average element size (Fig. 24).

Fig. 24

Local VEM spaces and degrees of freedom of a given element E

The first projection operator ${\Pi }_E^0 \nabla v$, which projects the gradient of v onto $[P_{k-1}(E)]^2$, satisfies the orthogonality condition as

$$\begin{aligned} \begin{aligned} <(\nabla v-{\Pi }_E^{0} \nabla v),\varvec{p}> = 0, \ \forall \varvec{p} \in [P_{k-1}(E)]^2 \end{aligned} \end{aligned}$$

from the inner product and applying the divergence theorem, we have

$$\begin{aligned} \begin{aligned}&\int _{E}{\Pi }_E^0 \nabla {v}.\varvec{p}{\rm d}\varvec{x}=\int _E \nabla v .\varvec{p}{\rm d}\varvec{x} = \oint _{\partial E} v.\varvec{p}.\varvec{n}{\rm d}s \\&\quad -\int _E v\nabla . \varvec{p} {\rm d}\varvec{x} , \ \forall \varvec{p} \in [P_{k-1}(E)]^2. \end{aligned} \end{aligned}$$

(54)

By introducing a set of shape functions for the local VEM space $V_h (E)$, $\phi _i (x),i=1,\ldots ,n_v$, we can express ${\Pi }_E^0 \nabla v$ as

$$\begin{aligned} \begin{aligned} {\Pi }_E^0 \nabla {v} = \sum _{i=1}^{n_v}{\Pi }_E^0 \nabla {\phi _i} (\varvec{x})V_i. \end{aligned} \end{aligned}$$

(55)

Therefore, Eq. (54) can be rewritten as

$$\begin{aligned} \begin{aligned}&\int _{E}{\Pi }_E^0 \nabla {\phi _i}.\varvec{m}_\alpha ^{(k-1)} {\rm d}\varvec{x}=\int _E \nabla \phi _i .\varvec{m}_\alpha ^{(k-1)} {\rm d}\varvec{x} \\&\quad = \oint _{\partial E} \phi _i.\varvec{m}_\alpha ^{(k-1)}.\varvec{n}ds -\int _E \phi _i\nabla . \varvec{m}_\alpha ^{(k-1)} {\rm d}\varvec{x}. \end{aligned} \end{aligned}$$

(56)

We can also express ${\Pi }_E^0 \nabla \phi _i$ using the set of basis $\varvec{m}_\alpha ^{(k-1)}$ for $[P_{k-1}(E)]^2$ as

$$\begin{aligned} \begin{aligned} {\Pi }_E^0 \nabla {\phi _i}(\varvec{x}) = \sum _{\beta =1}^{n_{p_{k-1}}}S_{i\beta }\varvec{m}_\beta ^{(k-1)}(\varvec{x}). \end{aligned} \end{aligned}$$

(57)

Finally, Eq. (56) can be rewritten as

$$\begin{aligned} \begin{aligned}&\sum _{\beta =1}^{n_{p_{k-1}}}S_{i\beta }\underbrace{\int _{E}{\varvec{m}_\beta }^{(k-1)}.{\varvec{m}_\alpha }^{(k-1)}{\rm d}\varvec{x}}_{\varvec{M}}\\&\quad = \underbrace{\oint _{\partial E} \phi _i.\varvec{m}_\alpha ^{(k-1)}.\varvec{n}{\rm d}s -\int _E \phi _i.\nabla . \varvec{m}_\alpha ^{(k-1)} {\rm d}\varvec{x}. }_{\varvec{R}} \end{aligned} \end{aligned}$$

(58)

From Eq. (58), we can form the matrices $\varvec{M}$ and $\varvec{R}$ and compute the matrix $\varvec{S}$ as

$$\begin{aligned} \begin{aligned} \varvec{S}=\varvec{R}\varvec{M}^{-1}. \end{aligned} \end{aligned}$$

(59)

1.2 Virtual element projection ${\Pi }_E^{\nabla } v$

The second projection operator ${\Pi }_E^{\nabla } v$, projects v onto lineal polynomial space $P_k (E)$, satisfies the orthogonality condition as

$$\begin{aligned} \begin{aligned} <\nabla (v-{\Pi }_E^{\nabla } v),\nabla p> = 0, \ \forall p \in P_{k}(E) \end{aligned} \end{aligned}$$

from the inner product and applying the divergence theorem, we have

$$\begin{aligned} \begin{aligned}&\int _{E}\nabla {\Pi }_E^{\nabla } v.\nabla p {\rm d}\varvec{x}=\int _E {\nabla } v.\nabla p {\rm d}\varvec{x} = \oint _{\partial E} v \nabla {p}.\varvec{n}{\rm d}s \\&\quad -\int _E v\Delta p {\rm d}\varvec{x} , \ \forall p \in P_{k}(E). \end{aligned} \end{aligned}$$

(60)

We express ${\Pi }_E^{\nabla } v$ using the shape functions $\phi _i$ as

$$\begin{aligned} \begin{aligned} {\Pi }_E^{\nabla } v= \sum _{i=1}^{n_v}{\Pi }_E^{\nabla } {\phi _i} (\varvec{x})V_i. \end{aligned} \end{aligned}$$

(61)

Therefore, Eq. (60) can be rewritten as

$$\begin{aligned} \begin{aligned}&\int _{E}\nabla {\Pi }_E^{\nabla } \phi _i.\nabla m_\alpha ^{(k)} {\rm d}\varvec{x}=\int _E {\nabla } \phi _i.\nabla m_\alpha ^{(k)} {\rm d}\varvec{x} \\&\quad = \oint _{\partial E} \phi _i\nabla m_\alpha ^{(k)}.\varvec{n}ds -\int _E \phi _i\Delta m_\alpha ^{(k)} {\rm d}\varvec{x}. \end{aligned} \end{aligned}$$

(62)

We can also express ${\Pi }_E^{\nabla } \phi _i$ in terms of the set of basis $m_\alpha ^{(k)}$ for $P_{k}(E)$ as

$$\begin{aligned} \begin{aligned} {\Pi }_E^{\nabla } \phi _i= \sum _{\alpha =1}^{n_{p_{k}}}S_{i\alpha }^{\nabla } m_\alpha ^{(k)}. \end{aligned} \end{aligned}$$

(63)

By Combining Eqs. (62) and (63), we obtain

$$\begin{aligned} \begin{aligned} \sum _{\beta =2}^{n_{p_{k}}} S_{i\beta }^{\nabla }\underbrace{\int _{E}\nabla m_\beta ^{(k)}.\nabla m_\alpha ^{(k)} {\rm d}\varvec{x}}_{\varvec{M}^\nabla } = \underbrace{ \oint _{\partial E} \phi _i\nabla m_\alpha ^{(k)}.\varvec{n}{\rm d}s -\int _E \phi _i\Delta m_\alpha ^{(k)} {\rm d}\varvec{x}.}_{\varvec{R}^\nabla }. \end{aligned} \end{aligned}$$

(64)

From Eq. (64), we can form matrices $\varvec{M}^\nabla$ and $\varvec{R}^\nabla$ and compute the matrix $\varvec{S}^\nabla$ as

$$\begin{aligned} \begin{aligned} \varvec{S}^\nabla =\varvec{R}^\nabla (\varvec{M}^{\nabla })^{-1}. \end{aligned} \end{aligned}$$

(65)

Additionally, we can express ${\Pi }_E^{\nabla } \phi _i(\varvec{x})$ in terms of the shape functions $\phi _i$ as

$$\begin{aligned} \begin{aligned} {\Pi }_E^{\nabla } \phi _i (\varvec{x})= \sum _{i=1}^{n_v}P_{ij}^{\nabla } \phi _j (\varvec{x}), \end{aligned} \end{aligned}$$

(66)

and we can express the set of basis functions $m_\alpha ^{(k)}$ for $P_{k}(E)$ as

$$\begin{aligned} \begin{aligned} m_{\alpha }^{(k)}(\varvec{x})= \sum _{j=1}^{n_v} G_{\alpha j}^{\nabla } \phi _j (\varvec{x}). \end{aligned} \end{aligned}$$

(67)

By substituting Eq. (67) into Eq. (63), we obtain

$$\begin{aligned} \begin{aligned} {\Pi }_E^{\nabla } \phi _i = \sum _{\alpha =1}^{n_{pk}} S_{i\alpha }^{\nabla } \sum _{j=1}^{n_{v}} G_{\alpha j}^{\nabla } \phi _j (\varvec{x}) = \sum _{\alpha =1}^{n_{pk}} \sum _{j=1}^{n_{v}} \underbrace{S_{i\alpha }^{\nabla }G_{\alpha j}^{\nabla }}_{P_{ij}^{\nabla }} \phi _j (\varvec{x}). \end{aligned} \end{aligned}$$

(68)

Then, $\varvec{P}^{\nabla } = \varvec{S}^{\nabla }\varvec{G}^{\nabla }$, where

$$\begin{aligned} \begin{aligned} \varvec{G}^{\nabla } = \begin{bmatrix} 1 &{} \ldots &{} 1 \\ m_{2}^{(1)}(\varvec{x}_1) &{} \ldots &{} m_{2}^{(1)}(\varvec{x}_{n_v}) \\ m_{3}^{(1)}(\varvec{x}_1) &{} \ldots &{} m_{3}^{(1)}(\varvec{x}_{n_v}) \\ \end{bmatrix} . \end{aligned} \end{aligned}$$

It is important to mention that $m_\alpha ^{(1)} (\varvec{x}_i)$ indicates the $\alpha$th basis function for $P_1 (E)$ evaluated at position $\varvec{x}_i$ of the vertex i.