This is an outdated version published on September 21, 2023. Read the most recent version.

Optimisation of pipes with constant diameter using the heuristic optimality criterion

Article

Authors & roles Authors

Peer review record Peer review

Metrics

Cite

Section

Research Article

Published

September 21, 2023

DOI

https://doi.org/10.12688/openreseurope.15943.1

Versions

April 9, 2026 (Version of Record 1.0)
September 21, 2023 (Published Manuscript Under Review 1.0)

European Commission

Grant ID(s): 831851

Views

180

Downloads

Citations

See full metrics

Lattice Boltzmann method, shape optimization, heuristic optimality criterion, laminar internal fluid flow, pressure drop minimization, Cellular Automaton

Abstract

Background: Minimising internal pressure drop in pipes is crucial for energy efficiency of fluid flow applications. Numerous computational optimisation tools that are capable of modifying flow geometries to reduce the pressure drop have been developed. Among these is a comparably simple heuristic optimisation algorithm which mimics erosion and sedimentation processes based on the shear stress in the vicinity of the domain boundaries. Although this method succeeds in modifying flow geometries for reduced pressure drop, it allows the fluid domain to widen during the reshaping process. Therefore, a reported reduction of pressure drop is not only caused by an improvement of the flow path, but also by an increase in the domain width. However, pipes with a constant circular diameter are favoured in many applications because they can be easily manufactured.

Methods: Here we combine the heuristic optimisation approach with a new geometrical constraint that kept the average diameter constant. This way, a reduction of pressure drop caused solely by the modification of the flow path can be assessed. We determined the applicability of the modified algorithm for 2D channel and 3D pipe geometries, conducting numerical simulations using the Lattice Boltzmann method with Reynolds numbers ranging from 40 to 500.

Results: For all investigated cases, potentially optimal shapes could be derived at all tested Reynolds numbers. However, the shape originally derived at a simulated flow with Re = 40 yielded a smaller pressure drop even for flows at Re ≥ 100 than shapes derived specifically at these higher Reynolds numbers.

Conclusions: This observed trend should be kept in mind when employing this approach as a simple way to improve channel and pipe layouts.

Corresponding author

David Blacher

0000-0003-3617-5374

david.blacher@ait.ac.at

Grant information

The work presented here was conducted in the context of the Clean Sky 2 project COMBO3D. This project has received funding from the Clean Sky 2 Joint Undertaking under the European Union’s Horizon 2020 research and innovation programme under grant agreement No 831851. This publication reflects only the authors’ views and the European Union is not liable for any use that may be made of the information contained therein.

The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

This work is licensed under a Creative Commons Attribution 4.0 International License.

First version published

September 21, 2023 (https://doi.org/10.12688/openreseurope.15943.1)

Latest version published

April 9, 2026 (https://doi.org/10.12688/openreseurope.15943.2)

Plain language summary

A high pressure drop in pipes limits the energy efficiency of machines, and there are computer tools that can change the pipe design to reduce the pressure drop. One method mimics erosion and sedimentation and is called heuristic optimality criterion. However, this method widens the pipe during the design process. This is problematic, because pipes with varying widths are harder to be produced. In this article, we describe how to restrict the reshaping method to leave the pipe width unchanged. We tested this method using computer simulations for 2D channels and 3D pipes with flows of different behaviours described by their Reynolds number. The Reynolds number describes whether the flow is more likely to be laminar or turbulent. In most cases, the method improved the shape of the pipes while keeping the diameter fixed. However, we found some shortcomings of the method for laminar flows with Reynolds numbers 100 and higher that should be considered for its application.

1 Introduction

1.1 State of the art

Fluid flow applications such as heat exchangers, fluid distributors and gas ducts must be carefully designed to minimise pressure drop in order to maximise energy efficiency. Many approaches to optimise geometries exist, such as the adjoint^1,2 or conjugate gradient³ method. These methods are capable of finding solutions for mathematical optimisation problems and allow for the combination of multiple objective- and cost-functions, making these approaches appropriate for a wide range of applications, including for turbulent flows⁴ and the design of flow manifolds⁵. However, a geometry optimised in such a way might not always be feasible or economically viable because of practical, geometrical or material-specific restrictions. Adding such constraints in a convenient mathematical form can be challenging and requires an experienced user with holistic knowledge about fluid dynamics, mathematical optimisation and the unique demands of the application. Wang et al.⁶ proposed a different reshaping algorithm based on a relatively simple "heuristic optimality criterion". The shape of an initial flow geometry is modified in an iterative process: first, the fluid flow through the geometry is simulated using the Lattice Boltzmann Method (LBM). Then, the shape is slightly altered, mimicking erosion and sedimentation depending on the current flow field. These steps are repeated until a minimum pressure drop hence as optimum is found. The reshaping is governed by the following rules, implemented as a cellular automaton (CA):

The fluid volume is kept constant throughout the reshaping process.
The geometry only changes at the fluid-solid interface.
Erosion occurs where viscous stress near wall nodes is highest. The neighbouring solid nodes become fluid.
Sedimentation is applied to fluid nodes with the smallest value of dynamic pressure.

Wang et al. tested their method on 2D elbows and bifurcations and found a considerable reduction of pressure drop compared to the initial shapes across a range of Reynolds numbers from Re = 0.0267 to Re = 80. The applicability of the method was then proven by Wang et al. for the design of a fluid distributor⁷. Later, Tarlet et al. made use of the same method to derive a suitable geometry for a fluid mixer⁸. A similar approach was taken by Park et al.⁹, but they chose the strain rate as criterion for both erosion and sedimentation. Park et al. tested their method on 2D and 3D rectangular ducts at Reynolds numbers up to Re = 100. The resulting 2D elbow shapes are in close accordance with those Wang et al. published. Tao et al.¹⁰ also implemented the heuristic optimality criterion but they combined it with additional optimisation approaches. They start with an empty box as initial fluid domain, within which an optimal flow can form unhindered by wall boundaries. Then, the fluid domain is shrunken by substituting fluid cells with solid zones where dynamic pressure is lowest until a certain volume is reached. During a second phase, the shape is further enhanced by erosion and sedimentation, with the volume being held constant similar to the heuristic optimality criterion as proposed by Wang et al. Unfortunately, Tao et al. did not specify how this volume was chosen in their publication, which is why a direct comparison of the results is non-trivial. Interestingly, Tao et al. adapted the method for ANSYS Fluent, a finite volume based solver. While the use of a finite volume based solver has some advantages over the Lattice Boltzmann Method (LBM see section subsection 2.1) such as fast convergence for steady-state simulations, it requires meshing of the fluid domain after each reshaping iteration. This time-consuming operation is not necessary when coupling a reshaping algorithm with the LBM.

1.2 Motivation - constant diameter constraint

The methods developed by Wang et al.⁶, Park et al.⁹ and Tao et al.¹⁰ offer simple ways to improve flow geometries. However, there is a major drawback limiting the applicability of the resulting geometries: in the presented cases, applying the reshaping algorithm results in a widening of the fluid domain diameter. This is because the path between inlet and outlet becomes shorter as the reshaping is applied. Due to the constant fluid volume constraint, the accumulating surplus of fluid nodes causes a widened cross section. This also implies that the reduction of pressure drop reported in the course of the reshaping was not solely caused by an improvement of the flow path, but also by the increased domain width. Further increasing the width would likely cause the pressure drop to reduce even more. However, these geometries would not be favoured in real life applications where the diameter of pipes is often limited by a number of factors: available space, requirements regarding the stability of the pipes and material costs limit the possible size of the geometry. Moreover, the manufacturing process of pipes with varying diameter would require advanced techniques. Therefore, a new geometrical constraint is presented to counteract the diameter dilation. In the first section, the Lattice Boltzmann Method (LBM) is briefly introduced. Then, the numerical setup is presented followed by a detailed description of the reshaping algorithm and the new reshaping constraint. The optimisation method was implemented as cellular automaton (CA) into the existing open-source LBM library called Palabos¹¹.

2 Mathematical model

Mathematical symbols used throughout the text are summarised in Table 1.

Table 1. Mathematical symbols frequently used.

Bold letters represent vectors (e.g. r). Capital letters with two subscript indices represent second-order tensors (e.g. S_ij), lower-case letters their elements (e.g. s_ij). Symbols with a tilde (e.g. t̃) denote dimensionless quantities.

Symbol	Unit	Meaning
Δp	1	relative pressure drop
δt	1	numerical time step
ρ	kg m^-3	volumetric mass density
ξ_i	m s^-1	particle velocity set of lattice model
d	m	diameter
d₀	m	initial diameter
f_e	1	fraction of eroding nodes
f_s	1	fraction of sedimenting nodes
t_a	t	adaptation time
f_i (r , ξ , t)	s³ m^-6	particle distribution function
M	1	surface area in voxel units
p	N m^-2	pressure
r	m	position
Re	1	Reynolds number
S_ij	s^-1	strain rate tensor
∥S_{i, j}∥	s^-2	Frobenius norm of the strain rate tensor
u₀	m s^-1	reference velocity

2.1 Lattice Boltzmann Method

An elaborate description of the lattice Boltzmann method (LBM) can be found for example in Sauro Succi's book¹² on which this short introduction is based on.

The LBM is a numerical model used for fluid flow simulations. It models the incompressible Navier-Stokes equations by solving a discretised Boltzmann transport equation, describing the evolution of the particle distribution function f_i due to particle interactions. The particle population f_i at each node r describes the number of particles moving with one velocity ξ_i out of a set of possible velocities. This set is designed in such a way that particles move from one node to a neighbouring node exactly within one numerical time step. For reasons of numerical stability and accuracy, a 3-dimensional model with 19 velocities (D3Q19) was used for all our test cases. The corresponding set of velocities ξ_i and their weights w_i of the D3Q19 lattice model are described by Figure 1 and Table 2. Aside from particles at rest (i = 0), all other velocities (i ∈ [1, 18]) point towards one of 18 neighbouring nodes. During one time step δt, the particles at each node f_i(r) collide and then stream to neighbouring nodes:

f_{i} (r + ξ_{i} δ t, t + δ t) = f_{i} (r, t) + Ω_{i} (r, t) (1)

Figure 1. Discrete velocities of the D3Q19 lattice.

To help visualise the phase space, each velocity is referred to by letters indicating its direction: "N,E,S,W,T,B,C" stand for "North, East, South, West, Top, Bottom, Center".

Table 2. Velocities ξ_i and weights w_i for the D3Q19 model.

Velocities are scaled by the lattice velocity u_lb. The notation (x, y, z)_c indicates all possible permutations of the given tuples.

ξ_i	w_i	i ∈
(0, 0, 0)	1/3	{ C }
(±u_lb, 0, 0)_c	1/18	{ E, N, W, S, T, B }
(±u_lb, ±u_lb, 0)_c	1/36	{ NE, NW, SE, SW, TE, TW, BE, BW, TN, TS, BN, BS }

For practical reasons, Equation 1 is separated in two computational steps. First the post-collision populations $f_{i}^{c}$ are calculated at the current location r:

f_{i}^{c} (r, t + δ t) = f_{i} (r, t) - \frac{1}{τ} (f_{i} - f_{i}^{e q}) (2)

Then, the populations at each node are updated based on the incoming particles from neighbouring nodes:

f_{i} (r + ξ_{i} δ t, t + δ t) = f_{i}^{c} (r, t + δ t) (3)

In Equation 2, the post-collision particle interaction $f_{i}^{c}$ is calculated by approximating the collision operator Ω as local relaxation towards equilibrium, the lattice Maxwell-Boltzmann distribution $f_{i}^{eq}$ . This is known as the BGK model.

f_{i}^{e q} (r) = ρ w_{i} [1 + \frac{1}{c_{s}^{2}} (ξ_{i} \cdot u) + \frac{1}{2 c_{s}^{4}} {(ξ_{i} \cdot u)}^{2} - \frac{1}{2 c_{s}^{2}} (u \cdot u)] (4)

The relaxation time τ depends on the chosen lattice viscosity ν_{l b}.

c_{s} = \frac{1}{\sqrt{3}} τ = \frac{v_{l b}}{c_{s}^{2}} + 1 / 2 (5)

From the distribution function f_i, all relevant macroscopic flow properties can be calculated:

ρ = \sum_{i = 0} f_{i} p = c_{s}^{2} \cdot ρ (6)

u = \frac{1}{ρ} \sum_{i = 0}^{q} f_{i} ξ_{i} S_{i j} = - \frac{1}{2 τ c_{s}^{2} ρ} \sum_{i = 0} (f_{i} - f_{i}^{(e q)}) ξ_{i} ξ_{i} (7)

The incompressible fluid is simulated under isothermal conditions, with isotropic, homogeneous and constant physical properties. The simulations conducted for this paper utilised regularised BGK dynamics, a multiple relaxation time (MRT) model. Latt et al.¹³ showed that this enhanced model can reproduce analytically known solutions more accurately than the basic BGK model by restoring symmetry properties of the populations before every collision step. Table 3 lists the numerical parameters of the simulations for each chosen Reynolds number (Re). To ensure comparability with results by Wang et al. and Park et al., the resolution was set to the same value, given by the number of nodes across the characteristic length (the inlet diameter) N_cells,0 = 1/δ_x = 40.

Table 3. Chosen numerical parameters given in lattice units.

Re	Velocity u_lb	Viscosity ν_lb	N_cells_,0 = 1/δ_x
40	0.04	0.04	40
100	0.05	0.02	40
200	0.05	0.01	40
400	0.05	0.005	40
500	0.05	0.004	40

2.2 Simulation setup

The initial pipe elbow geometry is created by two cylinders of the same diameter aligned perpendicularly to each other. One lies along the x-, the other along the y-axis. They extend from the inlet and outlet to the diagonal plane where they intersect. The simulation domain is a cuboid comprised of N_total = N_x · N_y · N_z nodes. See Table 4 for the size of the domain for each geometry. The XY-plane at z = N_z/2 of the initial 3D elbow corresponds to the initial 2D elbow shape (see Figure 2). In the 2D simulations, a channel with infinite height is simulated by applying periodic boundary conditions in z-direction. At the outlet, a constant pressure condition p̃_out = 1 is imposed. At the inlet, a parabolic velocity profile is imposed. Before the reshaping algorithm is applied, the flow field is initialised: starting from zero, the inlet velocity increases step-wise until it reaches the reference velocity ũ₀ = 1 in its centre. Then the flow is further simulated until the flow field converged to a steady state. This point in time is determined based on the global change of kinetic energy. Between each reshaping iteration, the fluid flow is also simulated for a time t̃_a. This is necessary for the flow to adapt to the new geometry. Ideally, each reshaping step is executed once a steady state has been reached. However, this approach would require an impractical amount of computing time. Instead, a suitable period was determined by looking at the convergence of the strain rate along the walls after a single reshaping step. This lead to the choice of approximately 1.5 times the time needed for one flush through the geometry. The chosen adaptation times for each geometry are listed in Table 4.

Table 4. Optimisation parameters: adaptation time between reshaping iterations t̃_a is given as dimensionless value.

The fraction of wall-neighbouring nodes subject to erosion in each iteration is denoted as f_e. The domain size consists of the total number of computational cells Nx,Ny,Nz.

	Optimisation		Domain size
Type	t̃_a	f_e	Nx	Ny	Nz
2D channel	15	0.02	200	200	4
3D pipe	25	0.01	400	320	120

Figure 2. Initial 2D elbow geometry.

Boundary conditions: prescribed velocity u at inlet; constant pressure at outlet. Red lines indicate where pressure is calculated for evaluation of pressure drop.

Wang et al.⁶ and Park et al.⁹ chose a different approach. After each reshaping step, they initialised the whole flow field again from zero velocity. This lengthy procedure was considered disproportionate regarding the computational costs because the shape only changes little in each reshaping step.

3 Method

The individual steps of the optimisation procedure will be addressed in the following. As a reference, refer to Figure 5 which describes the implemented algorithm in a flowchart.

3.1 Erosion and sedimentation model

The reshaping algorithm allows modifications of the fluid domain only at the fluid-solid interface. Fluid nodes along the interface will be referred to as wall-neighbouring nodes in the remaining text. The reshaping is guided by a CA based on the strain rate. To determine which of the wall-neighbouring nodes should be affected by the reshaping, they are ranked by the local value of ||S_i,j||, the Frobenius norm of the strain rate tensor. Sedimentation will be applied to a fraction f_s of wall-neighbouring nodes where ||S_i,j|| is smallest. Erosion will be applied to a fraction f_e of solid nodes directly neighbouring the fluid nodes which suffer from the greatest stress ||S_i,j||. An eroding solid node joins the fluid domain by changing its dynamics from bounce-back to BGK. For simplicity, the resulting new fluid node is initialised with zero velocity. If a fluid node sediments, it becomes a bounce-back node. In order to change the speed of the entire optimisation process, the fractions of affected nodes can be scaled: allowing more nodes to erode and to sediment in each step reduces the number of overall reshaping steps required to reach an optimised shape. However, Wang et al. recommend to choose a small fraction of affected nodes in order to ensure the shape actually converges to its optimum. Wang et al. chose the fraction of nodes affected by either erosion or sedimentation (f_e + f_s) as 0.04 of all wall-neighbouring interface nodes n_f. Park et al. varied the fraction between 0.02 and 0.05 in their simulations. The fractions of eroding nodes chosen for the present study is listed in Table 4. The ratio of eroding nodes to sedimenting nodes f_e/f_s determines if the fluid geometry may grow or shrink in each iteration. This balance is governed by the chosen reshaping constraint: either keeping the overall fluid volume or the geometry diameter constant throughout the reshaping. In order to keep the fluid volume constant, the total number of fluid nodes must not change, requiring f_e/f_s = 1 throughout the entire optimisation process. A constant geometry diameter can be achieved by adjusting the ratio of eroding nodes to sedimenting nodes during the reshaping process.

3.2 Calculation of current diameter

At each iteration, the diameter is evaluated by considering the relationship between volume and surface area of the fluid domain. This is based on a generalisation of the theorem of Pappus¹⁴, which allows the calculation of surface area and volume for a solid generated by moving a plane figure F (enclosed by curve C) perpendicularly along an arbitrary curve L. Figure 3 depicts how both the initial and an exemplary optimised pipe geometry can be described as such a solid. Its volume V equals the product of the area A of figure F and distance l_F travelled by the centroid of F:

V = A \cdot l_{F} (8)

Figure 3. Geometrical analysis: by tracing figure F (enclosed by curve C) along the path L (which lies in the XY-plane), the 3D pipe geometry can be generated.

(a) Initial geometry. Curve C is depicted 3 times: at the inlet, the intersection of cylinders and at the outlet. (b) Optimised geometry for Re = 40. Curve C is depicted 2 times: at the inlet and outlet.

Provided the traced curve L lies in a plane, the generated body’s surface area M equals the product of the arc length s of C and the distance l_C travelled by the geometric centroid of C:

M = s \cdot l_{C} (9)

Consequently, we can write the ratio of a body’s volume to its surface area as:

V / M = \frac{A \cdot l_{F}}{s \cdot l_{C}} (10)

For the simple cases of figure F being either a circle or rectangle, the centroids of figure F and curve C coincide, and l_F = l_C = l, resulting in:

\Rightarrow V / M = \frac{A}{s} (11)

For a circle with radius r = d/2 we get:

V / M_{◯} = \frac{r^{2} π}{2 π r} = \frac{r}{2} = \frac{d}{4} (12)

This corresponds to a pipe with circular cross section. A channel with infinite height (quasi 2D case) can be approximated as a slender (h >> d) rectangle with width d and height h:

V / M_{□} = \frac{d h}{2 (d + h)} = \frac{d}{2 (\frac{d}{h} + 1)} \approx \frac{d}{2} (13)

In both cases, the ratio V/M is proportional to the diameter d, but independent of the curve. Therefore, the average diameter of these geometries can be gauged by V/M even if the path l changes. This calculation method is not suited for complex geometries, such as bifurcations or 3D geometries with non-circular cross sections. For the conducted simulations, the change of the diameter relative to its initial value d/d₀ was of interest. This ratio can be calculated with knowledge of the initial values of volume V₀ and surface area M₀ and their current values V and M:

\frac{V / M}{V_{0} / M_{0}} = \frac{d}{d_{0}} = \frac{\tilde{d}}{{\tilde{d}}_{0}} (14)

In dimensionless values the initial diameter d̃₀ is 1, and the current diameter is calculated as

\to \tilde{d} = \frac{V / M}{V_{0} / M_{0}} (15)

The volume V is approximated by the sum of all fluid nodes, with a volume of V_n = 1 each: The surface area M is calculated from a generated smoothed STL representation of the fluid-solid interface¹⁵. The STL file format is commonly used to store a 3D representation of a surface comprised of unstructured triangulated sub-surfaces, each defined by its unit normal and vertices.

3.3 Adjusting the erosion/sedimentation ratio for a constant diameter

The measurement of the diameter allows to adjust the erosion/sedimentation ratio to counteract shrinking or widening during the reshaping process. Each reshaping step consists of the erosion process, a re-evaluation of the diameter and, if necessary, one or more layers of sedimentation. Sedimentation is applied incrementally until the diameter reaches its initial value. The number of sedimenting nodes for each increment is proportional to the current deviation of the diameter from its initial value:

f_{s} = {(\frac{d}{d_{0}})}^{3} - 1 for 3 D pipe (16)

f_{s} = {(\frac{d}{d_{0}})}^{2} - 1 for 2 D channel (17)

This ensures convergence of the diameter to its initial value. In order to allow regions with little fluid flow to sediment first, a new layer of solid nodes can be added on top of the last with each increment. Because the fluid-solid interface changes each time, a new ranking by local strain rate ||S_i,j|| is required to determine the nodes for the next increment. However, as no fluid flow is simulated between these incremental sedimentation steps, the ranking of the new wall-neighbouring nodes is based on strain rate information of the previous reshaping step. Figure 4 shows how the initial hard edge of the elbow geometry becomes rounded during the first reshaping step. Such a sedimentation of multiple layers occurs mainly in the first few reshaping steps.

Figure 4. Reshaping of the 2D elbow channel.

During the first two reshaping steps, multiple layers sediment. In the subsequent steps, only few multi-layer sedimentation steps take place.

Figure 5. Flowchart of the reshaping algorithm.

3.4 Pressure drop evaluation

The inlet pressure p̃_in was calculated as the average at the inlet domain boundary. The outlet pressure is defined as the reference pressure p̃_out = 1. The dimensionless absolute pressure drop Δp̃_a of a given fluid flow is the difference between inlet and outlet:

Δ {\tilde{p}}_{a} = {\tilde{p}}_{i n} - {\tilde{p}}_{o u t} = {\tilde{p}}_{i n} - 1 (18)

To allow the comparison of the pressure drop reduction in simulations with different Reynolds numbers, the relative pressure drop was calculated. The relative pressure drop for a given shape is the ratio of its current absolute pressure drop and the absolute pressure drop of the initial elbow shape.

Δ p = \frac{Δ {\tilde{p}}_{a}}{Δ {\tilde{p}}_{a, i n i t}} (19)

By this definition, the optimisation of any shape started at a relative pressure drop of Δp = 1. If the pressure drop continuously increased after hitting a minimum or the shape converged to a final shape, the simulation was terminated and the shape yielding the best result exported.

4 Results and discussion

Table 5 summarises the pressure drop achieved during the reshaping process. For all tested setups, shapes derived under the constant volume constraint yielded a smaller pressure drop at a wider diameter.

Table 5. Results of conducted simulations: Setup parameters are listed on the left hand side.

On the right hand side, information on the optimised shape are given. Optimised pressure drop ∆p_opt and diameter d are relative to values at initial elbow shape. The number of the reshaping steps after which a minimum pressure drop was reached is listed as s_opt. A ‘+’ symbol indicates that the shape was fully converged at this step.

Geometry	Constraint	Re	∆p_opt	d	s_opt
2D channel	V	40	0.322	1.287	960	+
	V	200	0.130	1.240	605	+
	V	500	0.036	1.230	619	+
	d	40	0.739	1.001	379
	d	200	0.472	1.002	771	+
	d	500	0.255	1.002	652	+
3D pipe	d	40	0.902	1.000	658	+
	d	100	0.972	1.000	125
	d	200	0.943	1.000	38
	d	400	0.832	1.000	147

4.1 2D channel elbow

A comparison of the optimised 2D channel at Re = 40 with the shapes reported by Wang et al. and Park et al. is depicted in Figure 6. All three shapes were derived with the same numerical parameters and under a constant volume constraint. An associated relative pressure drop Δp_opt = 67.8% is in good accordance with the values of Wang et al. (68.8%) and Park et al. (66%). This agreement suggests that the optimally criterion implemented in this thesis is valid. Figure 8 visualises the change in pressure drop over the course of the reshaping. In simulations with higher Reynolds numbers, an even smaller of Δp_opt could be achieved. Figure 7 shows that the curvature of the resulting shapes increases with rising Reynolds number. This behaviour was also reported by Wang et al. across a range of lower Reynolds numbers. They found a straight connection between inlet and outlet for Re = 0.0267, and arc like shapes at Re = 20, 40, 80. In all but one of the simulations, the pressure drop converges to a minimum. Only during the optimisation at Re = 40 under the constant diameter constraint, Δp_opt increased again by 2.6% after reaching a minimum at s = 379 reshaping steps. The evolution of shape during this simulation is depicted in Figure 9a. In all other 2D simulations, the geometry had reached a final shape and Δp_opt its final value after approximately s = 500 reshaping iterations. The evolution of shape during one such simulation is depicted in Figure 9b. At s = 1000 reshaping iterations, the simulations were terminated.

Figure 6. 2D elbow: optimised shape at Re = 40 for constant volume.

Comparison with existing publications^6,9. Our result is depicted here and in Figure 7a as solid green line.

Figure 7. 2D elbow: optimised shapes.

(a) Shapes derived under constant volume constraint. (b) Shapes derived under constant diameter constraint.

Figure 8. 2D elbow: evolution of relative pressure during the reshaping (at varying Reynolds numbers at either of the two reshaping constraints).

For each simulation, the pressure drop values are relative to their initial values. Simulations at the same Reynolds number start with the same absolute pressure drop (the difference between inlet and outlet). (a) Constant volume constraint. (b) Constant diameter constraint.

Figure 9. 2D elbow: snapshots taken at selected reshaping steps s during reshaping.

(a) Evolution of shape at Re = 40 and constant diameter. (b) Evolution of shape at Re = 500 and constant diameter.

4.2 3D pipe elbow

For the 3D pipe geometry, only reshaping with the constant diameter constraint was performed. Interestingly, the achieved relative pressure drop does not follow a clear trend. Figure 13 shows that Δp only converges to a minimum during the reshaping at Re = 40. To eliminate the influence of potential spatial limitations, a longer outlet was simulated by doubling the domain in the x-direction. 3D renderings of the derived shapes for Re = 100 and Re = 400 are visible in Figure 10. The shape derived at Re = 40 closely resembles those found in 2D at this Reynolds number (see Figure 11a): inlet and outlet are connected by an arc like shape. However, the shapes derived at Re ≥ 100 behave differently: the elbow dips below the level of the outlet, as visible in Figure 12a. Furthermore, the pipe cross section deviates from a circular profile: Figure 11b shows that for Re = 40, the cross section is still close to a circular profile. Above that, deviations increase with the Reynolds number, as visible in Figure 12b. The shear stress fields of the flow at Re = 200 in the initial and optimised 3D pipe are plotted in Figure 14 and Figure 15.

Figure 10. 3D visualisation of pipe optimised at two Reynolds numbers.

(a) Re = 100. (b) Re = 400. The pipe shows significant deformation.

Figure 11. 3D elbow: projections of shape derived at Re = 40.

(a) XY-plane projection of shape derived at Re = 40. The cross section in plane AA is depicted in figure b. The shape optimised at Re = 400 and the domain boundaries of the shorter 2D elbow are depicted for reference. (b) The cross section (plane AA) of the shape derived at Re = 40 retains an almost circular profile.

Figure 12. 3D elbow: projections of shapes derived at Re ≥ 100.

(a) XY-plane projection of shapes. The cross section in plane BB is depicted in figure (b). The domain boundaries of the shorter 2D elbow are depicted for reference. (b) The cross sections of shapes in plane BB (corresponds to YZ-plane at x̃ = 3, equivalent to n_x = 120). With increasing Reynolds number, the deviation from a circular profile becomes more pronounced.

Figure 13. 3D elbow: evolution of relative pressure drop during reshaping at different Reynolds numbers.

Figure 14. Shear stress in initial 3D pipe elbow at Re = 200.

(a) Full geometry. Grey box indicates area of detailed view in figure below. (b) Detail of geometry.

Figure 15. Shear stress in optimised 3D pipe elbow at Re = 200.

(a) Full geometry. Grey box indicates area of detailed view in figure below. (b) Detail of geometry.

4.3 Convergence of shapes

For this study, the direct dependence of the shape optimisation on a variation of the initial geometry was not investigated, but some conclusions can still be drawn from the observed convergence behaviour. In any 3D simulation with Re ≥ 100, the shapes did not converge to reach a minimum pressure drop. Instead, an intermediate iteration with a minimum pressure drop could be identified by the algorithm. This is visible in Figure 13. Consequently, the intermediate shape with the lowest pressure drop could not have been returned if a subsequent iteration of the reshaping had been used as the initial shape. Therefore, there seems to be a dependency of the optimisation algorithm on the initial geometry for 3D geometries at Re ≥ 100 which demands to be explored in more detail. While Wang et al.⁶ reported that the optimisation of 2D shapes at Re = 40 was largely unaffected by the initial geometry, this could not be verified for the reshaping with the new constant diameter constraint at Re = 40: here the smallest pressure drop was found after s = 379 iterations (see Table 5 and Figure 9a). For all other 2D cases, the geometry converged to a minimum. Overall, this indicates that the algorithm is indeed dependent on the starting geometry.

4.4 Quality of derived shapes

For both constraints and all setups (2D and 3D), shapes could be derived which yielded a reduced pressure drop compared to the initial geometries. The shapes derived at Re = 40 exhibit the shortest connection between in- and outlet. To rule out that these shapes perform better than the longer connections derived at Re ≥ 100, additional flow simulations were performed: for each tested Reynolds number Re ≥ 100, the flow through the reference shape derived at Re = 40 was analysed and Δp_Re=40 calculated. The values are summarised in Table 6. Surprisingly, the reference shape derived at Re = 40 resulted in the lowest pressure drop for flows, across all tested Reynolds numbers. This behaviour was found both among the shapes derived with the constant diameter and the constant volume constraint. These findings indicate that the algorithm as adapted from Wang et al. and Park et al. is not capable to find an intended optimum at Re ≥ 100, neither in 2D nor in 3D.

Table 6. Comparing pressure drop ∆p at a given Reynolds number through the dedicatedly derived shape and the shape derived at Re = 40.

	Flow at Re	Shape optimised at Re					∆p_opt/ ∆p_Re=40
2D channel		40		200		500
const. V	40	0.322		−−		−−	−−
	200	0.122		0.130		−−	1.065
	500	0.027		−−		0.036	1.335
const. d	40	0.739		−−		−−	−−
	200	0.452		0.472		−−	1.045
	500	0.206		−−		0.255	1.237
3D pipe		40	100	200	400
const. d	40	0.902	−−	−−	−−		−−
	100	0.776	0.972	−−	−−		1.253
	200	0.746	−−	0.943	−−		1.265
	400	0.672	−−	−−	0.832		1.237

5 Conclusion

The aim of this research was to assess the potential pressure drop reduction by changing the fluid domain path without affecting its diameter at Re ≥ 40. This was partly achieved by combining the existing heuristic optimality criterion with a further constraint of keeping the average diameter constant. Simulations conducted using this procedure revealed some limitations of the algorithm: Firstly, the cross section of 3D pipes deviated from the initial circular profile at Re ≥ 100. Secondly, and more importantly, for any flow with Re ≥ 100, 2D and 3D shapes optimised specifically for this Reynolds number resulted in a higher pressure drop than a reference shape derived at Re = 40. This behaviour was observed not only for shapes derived with a fixed average diameter, but also when the fluid volume was kept constant, as proposed by Wang et al.⁶ and Park et al.⁹. We concluded that the capability of this heuristic criterion to derive improved shapes with a minimum pressure drop is limited at Re ≥ 100 to achieving improvements which do not necessarily represent the optimal solution. The method can still be successfully applied to find an improved channel or pipe shape for complex layouts where other approaches may fail or be computationally prohibitive.

Data and software availability

Underlying data

Zenodo: Data for SEROS publication. https://doi.org/10.5281/zenodo.7788202¹⁶.

This project contains the following underlying data:

rhobarInletAvg_analysis (Summary of achieved Pressure Drop Reduction (PDR) of all cases.)
pre-processing-cases.zip (Directories containing the following files directly generated with SEROS)
- – bestSeros (Iteration at PDR was highest.)
- – config2D_L.xml (Case configuration for creation of data.)
- – measurementsSeros0 (Fluid measurements taken at each reshaping iteration.)
- – plbLog.dat (Log file generated during simulation execution.)
- – rho*.ppm (Visualisation of fluid density fields at given iterations.)
- – vtk*.vti (Fluid field at highest PDR available for some cases.)
2D_L_Re*_PDR.txt (Identical with measurementsSeros0 for each 2D case.)
3D_L_Re*_PDR.txt (Identical with measurementsSeros0 for each 3D case.)

Data are available under the terms of the Creative Commons Attribution 4.0 International license (CC-BY 4.0)

Software availability

Source code available from: https://github.com/AIT-LKR/SEROS

David Blacher. (2023). AIT-LKR/SEROS: First Release with full source code (v1.0-full-source-code).

Archived source code at time of publication: https://doi.org/10.5281/zenodo.7788368¹⁷.

License: The software is also available under a AGPL license.

Acknowledgements

The authors would like to thank Stephan Jäger for his conceptual contributions towards the method and support in development of the software.

Faculty Opinions recommended

References

1. Gersborg-Hansen A, Sigmund O, Haber RB: Topology optimization of channel flow problems. Struct Multidiscipl Optim. 2005; 30(3): 181–192. Publisher Full Text
2. Villanueva CH, Maute K: CutFEM topology optimization of 3d laminar incompressible flow problems. Comput Methods Appl Mech Eng. 2017; 320: 444–473. Publisher Full Text
3. Borrvall T, Petersson J: Topology optimization of fluids in stokes flow. Int J Numer Methods Fluids. 2003; 41(1): 77–107. Publisher Full Text
4. Dilgen CB, Dilgen SB, Fuhrman DR, et al.: Topology optimization of turbulent flows. Comput Methods Appl Mech Eng. 2018; 331: 363–393. Publisher Full Text
5. Behrou R, Ranjan R, Guest JK: Adaptive topology optimization for incompressible laminar flow problems with mass flow constraints. Comput Methods Appl Mech Eng. 2019; 346: 612–641. Publisher Full Text
6. Wang L, Fan Y, Luo L: Heuristic optimality criterion algorithm for shape design of fluid flow. J Comput Phys. 2010; 229(20): 8031–8044. Publisher Full Text
7. Wang L, Fan Y, Luo L: Lattice boltzmann method for shape optimization of fluid distributor. Comput Fluids. 2014; 94: 49–57. Publisher Full Text
8. Tarlet D, Fan Y, Luo L: Design and mixing performance characterization of a mini-channel mixer with nature-inspired geometries. Chem Eng Res Des. 2020; 164: 230–239. Publisher Full Text
9. Park J, Safdari A, Chun Kim K: Shape optimization of flow channels based on lattice boltzmann method. J Mech Sci Technol. 2018; 32(6): 2619–2627. Publisher Full Text
10. Tao M, Guo K, Huang Z, et al.: A hybrid optimization method to design shapes of three-dimensional flow channels. Chem Eng Res Des. 2016; 114: 190–201. Publisher Full Text
11. Latt J, Malaspinas O, Kontaxakis D, et al.: Palabos: Parallel lattice boltzmann solver. Comput Math Appl. Development and Application of Open-source Software for Problems with Numerical PDEs. 2021; 81: 334–350.Publisher Full Text
12. Succi S: The lattice Boltzmann equation: for fluid dynamics and beyond. Oxford university press, 2001. Reference Source
13. Latt J, Chopard B: Lattice boltzmann method with regularized pre-collision distribution functions. Math Comput Simul. 2006; 72(2–6): 165–168. Publisher Full Text
14. Goodman AW, Goodman G: Generalizations of the theorems of pappus. Am Math Mon. 1969; 76(4): 355–366. Publisher Full Text
15. Fuhrmann S, Kazhdan M, Goesele M: Accurate isosurface interpolation with hermite data. In: 2015 International Conference on 3D Vision. IEEE, 2015; 256–263. Publisher Full Text
16. Blacher D: Data for SEROS publication. Zenodo. [Data set]. 2023. http://www.doi.org/10.5281/zenodo.7788202
17. Blacher D: AIT-LKR/SEROS: First Release with full source code (v1.0-full-source-code). [Code] Zenodo. 2023. http://www.doi.org/10.5281/zenodo.7788368

Comments

Publication information

Section

Research Article

Published

September 21, 2023

DOI

https://doi.org/10.12688/openreseurope.15943.1

Versions

April 9, 2026 (Version of Record 1.0)
September 21, 2023 (Published Manuscript Under Review 1.0)

Included in

Funding data

European Commission

Grant ID(s): 831851

Metrics

Views

180

Downloads

Citations

See full metrics

Keywords

Lattice Boltzmann method, shape optimization, heuristic optimality criterion, laminar internal fluid flow, pressure drop minimization, Cellular Automaton

Comments

David Blacher

LKR Light Metals Technologies, AIT Austrian Institute of Technology, Lamprechtshausener Straße 61, Ranshofen, 5282, Austria

Roles

Corresponding Author, Co-Author, Conceptualization, Data Curation, Formal Analysis, Investigation, Methodology, Project Administration, Software, Validation, Visualization, Writing – Original Draft Preparation, Writing – Review & Editing

0000-0003-3617-5374

david.blacher@ait.ac.at

Michael Harasek

Institute of Chemical, Environmental and Bioscience Engineering, TU Wien, Getreidemarkt 9/166, Vienna, 1060, Austria

Roles

Co-Author, Conceptualization, Project Administration, Supervision, Validation, Writing – Original Draft Preparation, Writing – Review & Editing

0000-0002-6490-5840

michael.harasek@tuwien.ac.at

Downloads over the last 12 months:

Cited by

via Crossref

View citing articles

Optimisation of pipes with constant diameter using the heuristic optimality criterion

List of Authors

Peer review

Publication information

Included in

Funding data

Metrics

Keywords

Comments

Abstract

Comments

Publication information

Included in

Funding data

Metrics

Keywords

Comments

Authors & roles

David Blacher

Michael Harasek

Downloads over the last 12 months:

Cited by

via Crossref