Skip to main content
Log in

Fluid flow topology optimization in PolyTop: stability and computational implementation

  • EDUCATIONAL ARTICLE
  • Published:
Structural and Multidisciplinary Optimization Aims and scope Submit manuscript

Abstract

We present a Matlab implementation of topology optimization for fluid flow problems in the educational computer code PolyTop (Talischi et al. 2012b). The underlying formulation is the well-established porosity approach of Borrvall and Petersson (2003), wherein a dissipative term is introduced to impede the flow in the solid (non-fluid) regions. Polygonal finite elements are used to obtain a stable low-order discretization of the governing Stokes equations for incompressible viscous flow. As a result, the same mesh represents the design field as well as the velocity and pressure fields that characterize its response. Owing to the modular structure of PolyTop, incorporating new physics, in this case modeling fluid flow, involves changes that are limited mainly to the analysis routine. We provide several numerical examples to illustrate the capabilities and use of the code. To illustrate the modularity of the present approach, we extend the implementation to accommodate alternative formulations and cost functions. These include topology optimization formulations where both viscosity and inverse permeability are functions of the design; and flow control where the velocity at a certain location in the domain is maximized in a prescribed direction.

This is a preview of subscription content, log in via an institution to check access.

Access this article

Price excludes VAT (USA)
Tax calculation will be finalised during checkout.

Instant access to the full article PDF.

Fig. 1
Fig. 2
Fig. 3
Fig. 4
Fig. 5
Fig. 6
Fig. 7
Fig. 8
Fig. 9
Fig. 10
Fig. 11

Notes

  1. 1 Note, however, that not all Voronoi tessellations yield stable discretizations. For example, a structured quadrilateral mesh, corresponding to a Voronoi tessellation of a uniform grid of seeds, does not lead to a stable method.

  2. 2 The area/volume of the set A is denoted by \(\left |A\right |\).

  3. 3 It can also be seen from the weak form that terms \({\int }_{\Omega }p\,\text {div}\,\mathbf {v}\mathrm {d}\mathbf {x}\) and \({\int }_{\Omega }q\,\text {div}\,\mathbf {u}\mathrm {d}\mathbf {x}\) are not affected by the addition or subtraction of a constant to p or q. Indeed, for qc, a constant, and \(\mathbf {v}\in \mathcal {U}\) or \(\mathbf {v}\in \mathcal {V}\), we have \({\int }_{\Omega }q~\text {div}~\mathbf {v}\mathrm {d}\mathbf {x}=c{\int }_{\Omega }\mathbf {v}\cdot \mathbf {n}\mathrm {d}s=0\).

  4. 4 In Talischi et al. (2012b), P is used to denote the filtering matrix. In this paper, we will not use filtering so no confusion should arise.

  5. 5 In Matlab, if FreeDofs and PresDofs are arrays containing the indices of free and prescribed degrees of freedom, then M f U and M p U are simply obtained by U(FreeDofs,:) and U(PresDofs,:).

  6. 6 The strength of the singularity of the exact solution governs the convergence rates here.

  7. 7 We must note that, unlike optimization problem (12), this problem is ill-posed unless additional regularity is imposed on ρ, for instance, via filtering (Wiker et al. 2007). In the numerical result presented here, we do not consider filtering as we are mainly interested in comparing the results to the solutions of the optimization problem (12). However, we note that filtering can be easily enabled in PolyTop through the filtering matrix opt.P.

References

Download references

Acknowledgments

Glaucio H. Paulino acknowledges support from the US National Science Foundation through grants #1559594 (formerly #1335160) and #1624232 (formerly #1437535). Ivan F. M. Menezes and Anderson Pereira acknowledge the financial support provided by Tecgraf/PUC-Rio (Group of Technology in Computer Graphics), Rio de Janeiro, Brazil. The information presented in this publication is the sole opinion of the authors and does not necessarily reflect the views of the sponsors or sponsoring agencies.

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Glaucio H. Paulino.

Additional information

Electronic supplementary material The online version of this article (doi:10.1007/s00158-014-1182-z) contains supplementary material, which is available to authorized users.

Appendices

Appendix A:: Additions to PolyMesher

The polygonal finite element meshes used in this work are generated by PolyMesher (Talischi et al. 2012a) but with one main modification that allows to specify the exact location of the inlets and outlets. To this effect, a new function is added to the kernel in order to generate meshes that included user-defined fixed points as vertices in the mesh. The function PolyMshr_FixedPoints uses the property of the Voronoi diagram that an empty circle through three or more points defines a vertex. As shown in Fig. 12, consider seeds x and y and their reflections R Ω(x) and R Ω(y) are the two closest points in the mesh to a given fixed point z . If y has a smaller distance, then the function moves the seed x, R Ω(x) and R Ω(y) on a circle of radius \(r=\left |\mathbf {y}-\mathbf {z}^{*}\right |\) and center z . The new points \(\bar {\mathbf {x}}\), \(\bar {R}_{\Omega }(\mathbf {x})\) and \(\bar {R}_{\Omega }(\mathbf {y})\) are now positioned such that the desired point z is a vertex of the resulting Voronoi diagram.

Fig. 12
figure 12

Fixed vertices on the boundary are specified by adjusting the location of the nearby seeds and their reflections. Two scenarios are shown in the figure

The user must define a list of fixed points within the Domain function, more specifically, in a function FixedPoints, which returns a two-column array containing a list with the coordinates of the fixed points. The PolyMesher kernel function is also modified on line 10 as follows:

$$\begin{array}{@{}rcl@{}} \texttt{BdBox = Domain('BdBox');}\\ \texttt{PFix = Domain('PFix');} \end{array} $$

During each iteration of the Lloyd’s algorithm, a call is made to the functionPolyMshr_FixedPoints. In particular, the line P=Pc (line 15) is replaced by the following line of code:

$$\texttt{P = PolyMshr\_FixedPoints(Pc,PFix);} $$

The code PolyMshr_FixedPoints function is short and is provided below:

figure a

In order to enforce velocity boundary conditions, the domain function must identify the nodes that lie on the boundary Ω and assign prescribed velocity values to the associated degrees of freedom. This information is provided by a call to the auxiliary function PolyBoundary. This function uses the Matlab built-in function freeBoundary that queries the edges of a triangulation and outputs those incident on exactly one triangle. In PolyBoundary, each polygonal element of the mesh is first triangulated with respect to its centroid. The resulting triangulation is used as the input to freeBoundary. The list of boundary nodes is the list of nodes corresponding to the boundary edges output by freeBoundary.

Appendix B:: PolyScript

figure b

Appendix C:: PolyTop

figure c
figure d
figure e
figure f

Rights and permissions

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Pereira, A., Talischi, C., Paulino, G.H. et al. Fluid flow topology optimization in PolyTop: stability and computational implementation. Struct Multidisc Optim 54, 1345–1364 (2016). https://doi.org/10.1007/s00158-014-1182-z

Download citation

  • Received:

  • Revised:

  • Accepted:

  • Published:

  • Issue Date:

  • DOI: https://doi.org/10.1007/s00158-014-1182-z

Keywords

Navigation