COMSOL 6.4 - Minimizing the Flow Velocity in a Microchannel

Minimizing the Flow Velocity in a Microchannel

Topology optimization of the Navier–Stokes equations is encountered in different branches and applications, such as in the design of ventilation systems for cars and optimal reactors. A common technique applicable to such problems is to let the distribution of porous material vary continuously. In this model, the objective is to find the optimal distribution of a porous material in a microchannel such that the horizontal velocity component at the center of the channel is minimized.

The model is inspired by Ref. 1.

Model Definition

The model geometry (Figure 1) consists of three regions: the inlet channel, the design domain, and the outlet channel. A prescribed pressure drop between the inlet and the outlet drives the flow.

Figure 1: The model geometry.

The fluid flow in the channel is described by the Navier–Stokes equations

where −α(γ) u is a force term in which the coefficient

(1)

characterizes the flow in a porous medium. In Equation 1, α is interpreted as a continuous mapping determined by the function γ: Ω → [0, 1], which in the limit of decreasing Darcy number and decreasing mesh size should be a discrete-valued function. When γ equals 1, α is zero, corresponding to free flow. Conversely, for γ = 0, α = αmax, where αmax is related to the dimensionless Darcy number, Da, according to

The convergence of the optimization process depends on three important factors: the Darcy number, the mesh size, and the coefficient q. Rewriting Equation 1 it is easily seen that α/αmax → 1 − γ in the limit q → ∞. In this limit, γ can be interpreted as the local porosity, ranging between 0 (filled) and 1 (open channel).

Results and Discussion

Figure 2 displays the design variable, γ, which represents the distribution of porous material. As the plot shows, γ is either 0 or 1 in most of the domain, with a narrow transition zone in between. The width of this transition zone is mesh dependent; you can reduce it by decreasing the mesh-element size. Alternatively, decreasing the Darcy number also gives harder boundaries at the interface between porous and open domains.

Figure 2: Distribution of porous material; the red areas represent open channel.

A question that naturally arises in this type of problems concerns uniqueness of the solution. In this case, there is at least one more solution that gives exactly the same result; because the channel has no upside or downside, a solution mirrored around the axis y = 0.5 mm would give exactly the same flow.

Figure 3 contains a surface plot of the horizontal velocity component and a streamline plot of the velocity field resulting from the optimization process. In addition, the contour γ = 0.5 indicates the border between the open channel and filling material. The plots reveal how the flow turns around, with a negative horizontal velocity at the center of the channel. Note also that the x-velocity has a minimum of roughly −14 mm/s at the design point.

Figure 3: The horizontal velocity (surface plot) and velocity field (streamlines) after optimization. In addition, the contour γ = 0.5 indicates the border between open channel and filling material.

If you were to increase the streamline density in the above plot, some streamlines passing through the barriers would appear. This effect is due to a small amount of leakage, which can be reduced further by increasing the mesh resolution.

Figure 4 shows the pressure distribution, verifying the prescribed pressure drop through the channel length. The pressure drop is naturally concentrated to the region with porous material.

Figure 4: Pressure distribution in the channel, the pressure drop is concentrated to the porous domain.

Reference

1. L. Højgaard Olesen, F. Okkels, and H. Bruus, “A High-level Programming-language Implementation of Topology Optimization Applied to Steady-state Navier–Stokes Flow,” Int. J. Numer. Methods Eng., vol. 65, pp. 975–1001, 2005.

Notes About the COMSOL Implementation

This model combines a feature with a Laminar Flow interface. First, you calculate the solution for the flow in an empty channel (that is, with no porous material). You then solve the optimization problem. In each iteration, the software calculates a solution for the flow problem and feeds it to the optimization routine, which updates the design variables.

If the specified convergence criterion is fulfilled, the solution process terminates; otherwise the new design-variable values are used in the next calculation of the flow problem. However, for design problems such as this one, there is a tradeoff between computation time and convergence. The solution may be sufficiently improved long before convergence is reached in the mathematical sense. Therefore, it is useful to limit the number of steps taken by the optimization algorithm after which the solution can be evaluated and restarted if not yet satisfactory.

Setting up this kind of model with a general optimization routine requires quite a bit of work, but as you will discover, solving this problem with the built-in tools for optimization in COMSOL Multiphysics is easy.

Optimization_Module/Topology_Optimization/reversed_flow

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select > > .

Click .

Click

In the tree, select > .

Click

Global Definitions

Parameters 1

In the window, under click .

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Value	Description
q	1	1	Optimization parameter
meshsz	20[um]	2E-5 m	Center point mesh size
L	1[mm]	0.001 m	Inlet height

The purpose of the constant u0 is to introduce a rough velocity scale that you can use to make the objective function dimensionless with a value of the order 1.

Geometry 1

In the window, under click .

In the window for , locate the section.

From the list, choose .

Rectangle 1 (r1)

In the toolbar, click

In the window for , locate the section.

In the text field, type 2.5.

Rectangle 2 (r2)

In the toolbar, click

In the window for , locate the section.

In the text field, type 5.

Locate the section. In the text field, type 2.5.

Rectangle 3 (r3)

In the toolbar, click

In the window for , locate the section.

In the text field, type 2.5.

Locate the section. In the text field, type 7.5.

Point 1 (pt1)

In the toolbar, click

In the window for , locate the section.

In the text field, type 5.

In the text field, type 0.5.

Locate the section. Select the checkbox.

Click

Click the

button in the toolbar.

The geometry should now look like that in Figure 1.

Add Material

In the toolbar, click

to open the window.

Go to the window.

In the tree, select > .

Click the button in the window toolbar.

In the toolbar, click

to close the window.

Materials

Water, liquid (mat1)

Next use the feature to set up the design variable.

Component 1 (comp1)

Density Model 1 (dtopo1)

In the toolbar, click

and choose .

Only the center part of the channel geometry is needed in the optimization, so you only have to define the control variable there.

Select Domain 2 only.

In the window for , locate the section.

From the list, choose .

Locate the section. From the list, choose .

In the q text field, type q.

Locate the section. In the θ0 text field, type 1.

This defines the shape functions for the design variable used in the optimization. The initial value 1 corresponds to a channel free from porous material.

Definitions

Variables 1

In the toolbar, click

In the window for , locate the section.

From the list, choose .

Select Domain 2 only.

The maximum damping is limited by the mesh size, because the mesh has to resolve the transition between marginal flow in the porous material and free flow in fluid region.

Locate the section. In the table, enter the following settings:


Name	Expression	Unit	Description
alpha_max	4*spf.mu/meshsz^2	Pa·s/m²	Volume force coefficient, max value
alpha	alpha_max*dtopo1.theta_p	Pa·s/m²	Volume force coefficient