Particle Tracing in a Micromixer

Micromixers can either be static or dynamic depending on the required mixing time and length scale. For static mixers, the Reynolds number has to be suitably high to induce turbulence-enhanced mixing. Often micromixers operate in the laminar flow regime due to their small characteristic size. The diffusivity of a solute in the flowing fluid may also be extremely small, on the order of 10−10 m2/s. This results in mixing length scales on the order of meters — clearly unacceptable for a microscale device. One way to alleviate this problem is to add mixing elements to induce vorticity into the flow. A dynamic mixer uses rotating blades to enhance the mixing process, allowing for smaller-scale devices. The one big disadvantage of a dynamic mixer is that moving parts are required.

Model Definition

This example examines how mixing between microscopic particles occurs in a micromixer. Particles enter the mixer through 3 features and exit through the feature. The particles enter the modeling domain through the inlets in a continuous stream. A new set of particles is released every 50 milliseconds for a total duration of one second. After this, no more particles are released but the model is solved for an additional second. For each release inlet and each release time, 50 particles are released with an initial velocity equal to the fluid velocity, so a total of 3 × 21 × 50 = 3150 particles are released.

The geometry is an assembly containing stationary and rotating domains. The particles are free to cross the pair boundary between the stationary and moving domains as it if were invisible, provided that the feature is used in the Particle Tracing for Fluid Flow interface.

The blades are rotating at a constant angular velocity of 1 revolution per second in the anti-clockwise direction.

Figure 1: Plot of the model geometry. The geometry length unit is millimeters.

The particles obey Newton’s second law:

where

•

mp (SI unit: kg) is the particle mass,

•

v (SI unit: m/s) is the particle velocity, and

•

Ft (SI unit: N) is the total force on the particle.

In this example, the total force is dominated by the FD, for which Stokes’s law is

where

•

u (SI unit: m/s) is the fluid velocity,

•

μ (SI unit: Pa s) is the fluid dynamic viscosity, and

•

dp (SI unit: m) is the particle diameter.

In addition to the drag force, the optional virtual mass force Fvm and pressure gradient force Fp on the particle can also be considered. These forces are defined as

(1)

where mf (SI unit: kg) is the mass of the fluid displaced by the particle volume,

and ρ (SI unit: kg/m3) is the density of the fluid. In Equation 1 the derivative d/dt is a material (or total) derivative in the direction of the particle velocity, and D/Dt is a material derivative in the direction of the fluid velocity. That is, for an arbitrary vector field f,

The virtual mass and pressure gradient forces can usually be neglected when the density of the particle phase is much greater than the density of the fluid phase, as is true for solid particles in a gas. However, these forces might approach the same order of magnitude as the drag force if the particles are in a liquid. Particular attention should be paid to these forces when the flow is not stationary, since they each depend on both the spatial and time derivatives of the fluid velocity field.

The flow field is computed using the Laminar Flow interface. The force exerted on the fluid from the particles is neglected in this model. So, it is possible to solve for the flow field only in one study, then use a separate study to compute the particle trajectories based on that flow field. This is usually the recommended approach, if the field is computed from a study. In this case, there are very strong transients in the model, meaning that a huge number of time steps have to be stored if the model is to be solved sequentially. It is more attractive to solve for the particle trajectories and flow field in a single study step.

This geometry sequence is treated as an assembly rather than a union, so that the mesh in the inner domain containing the mixing blades can rotate freely. For the fluid flow, the feature must be added on the pair boundaries outside the rotating domain. For the particle tracing, the feature must be used on pairs. The mesh needs to be quite fine on the stationary/sliding interface so that the fluid motion remains continuous. The mesh used in this model is plotted in Figure 2.

Figure 2: The mesh is quite fine on the pair boundary to accurately resolve the flow field.

Results and Discussion

The location of the particles at different snapshots in time is plotted in Figure 3. The particle color is different for each particle , which conveniently allows the effect of the mixing to be visualized. The particles make their way normally inward from the inlets and, like the fluid velocity, begin to assume a parabolic profile. The particles entering from the left (the blue particles) are then swept downward due to the presence of the rotating blades. However, a few of these particles released at later times actually go clockwise, depending on the exact position of the blades when they first enter the mixer.

The particles entering from the right (red) are swept upward, but some of the particles go past the outlet because of the momentum they gained from the surrounding fluid. At about 0.6 seconds, particles from the bottom inlet (green) begin to reach the outlet as well. Mixing of the three particle streams continues even after new particles stop entering the domain at 1 s. This is because liquid continues to flow in from all of the inlets after the particle stream is terminated, and because the mixing blades continue to rotate.

Figure 3: Plot of the particle coordinates at different stages of the mixing process.

Reference

1. G. Karniadakis, A. Beskok, and N. Aluru, Microflows and Nanoflows, Springer, 2005.

Particle_Tracing_Module/Fluid_Flow/micromixer_particle_tracing

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select .

Click .

In the tree, select .

Click .

Click

In the tree, select .

Click

Geometry 1

The micromixer is only a few millimeters in size, so change the geometry length unit to millimeters:

In the window, under click .

In the window for , locate the section.

From the list, choose .

Circle 1 (c1)

In the toolbar, click

In the window for , locate the section.

In the text field, type 3.

Click

Circle 2 (c2)

In the toolbar, click

In the window for , locate the section.

In the text field, type 2.75.

Click

Difference 1 (dif1)

In the toolbar, click

and choose .

Select the object only.

In the window for , locate the section.

Find the subsection. Click to select the

toggle button.

Select the object only.

Select the check box.

Click

Rectangle 1 (r1)

In the toolbar, click

In the window for , locate the section.

In the text field, type 0.2.

In the text field, type 5.25.

Locate the section. From the list, choose .

Click

Rectangle 2 (r2)

In the toolbar, click

In the window for , locate the section.

In the text field, type 5.25.

In the text field, type 0.2.

Locate the section. From the list, choose .

Click

Rectangle 3 (r3)

In the toolbar, click

In the window for , locate the section.

In the text field, type 0.5.

Locate the section. In the text field, type -3.4.

From the list, choose .

Click

Rotate 1 (rot1)

In the toolbar, click

and choose .

Select the object only.

In the window for , locate the section.

In the text field, type 90 180 270.

Locate the section. Select the check box.

Click

Union 1 (uni1)

In the toolbar, click

and choose .

Select the objects , , , , and only.

In the window for , locate the section.

Clear the check box.

Click

Difference 2 (dif2)

In the toolbar, click

and choose .

Select the object only.

In the window for , locate the section.

Find the subsection. Click to select the

toggle button.

Select the objects and only.

Click

Form Union (fin)

The interface requires that a pair is present between the stationary and rotating domains. In order to do this, use the option. This will automatically create boundaries between the stationary and rotating domains.

In the window, under click .

In the window for , locate the section.

From the list, choose .

Click

Click the

button in the toolbar. The geometry should look like Figure 1.

Definitions

It is usually convenient to define an explicit selection for the pair boundaries.

Pair boundaries

In the toolbar, click

In the window for , type Pair boundaries in the text field.

Locate the section. From the list, choose .

Select Boundaries 15–18 and 33–36 only.

The easiest way to select these boundaries is to copy the text ’15-18, 33-36’, click in the selection box, and then press Ctrl+V. Alternatively, click the button and type or paste the boundary numbers in the dialog box that appears.

Now define a function for the inlet velocity. The boundary condition for the inlet velocity must be consistent with the initial condition for the velocity. The initial velocity in this model will be zero so the inlet velocity must be ramped up from zero to its maximum value over a certain period of time. In this case the ramp time is 0.01 seconds. To achieve this, the function is used with a of 100, meaning that the ramp function reaches its maximum value after 0.01 seconds.

Ramp 1 (rm1)

In the toolbar, click

and choose .

In the window for , locate the section.

In the text field, type 100.

Select the check box.

Click to expand the section.

Select the check box. In the associated text field, type 0.001.

Now that the function is defined, create an expression for the inlet velocity which will ramp up over 0.01 seconds.

Variables 1

In the toolbar, click

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Unit	Description
uin	0.02[m/s]*rm1(t[1/s])	m/s	Inlet velocity

Materials

Material 1 (mat1)

In the window, under right-click and choose .

In the window for , locate the section.

In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Density	rho	1E3	kg/m³	Basic
Dynamic viscosity	mu	1E-3	Pa·s	Basic

Add a feature which designates the rotating domain. The speed of revolution is also specified, in this case one revolution per unit time. This means the blade system will undergo one complete revolution (360 degrees) per second.

Moving Mesh

Rotating Domain 1

In the window, under click .

In the window for , locate the section.

Click

Select Domain 2 only.

Locate the section. In the f text field, type 1.

Laminar Flow (spf)

Inlet 1

In the window, under right-click and choose .

Select Boundaries 1, 5, and 12 only.

In the window for , locate the section.

In the U0 text field, type uin.

Outlet 1

In the toolbar, click

and choose .

Select Boundary 7 only.

Particle Tracing for Fluid Flow (fpt)

Wall 1

In the window, under click .

In the window for , locate the section.

From the list, choose .

Start by adding the drag force on the particles. This requires input of the fluid velocity and viscosity.

Drag Force 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. From the u list, choose .

From the μ list, choose .

Locate the section. Select the check box.

Now define a stream of particles over the first second for each inlet, with 50 particles per inlet and a new release every 50 milliseconds. Defining 3 separate inlet features will allow for improved visualization during results processing.

Inlet 1

In the toolbar, click

and choose .

Select Boundary 1 only.

In the window for , locate the section.

From the list, choose .

In the N text field, type 50.

Locate the section. From the u list, choose .

Locate the section. Click

In the dialog box, type 0 in the text field.

In the text field, type 1.

In the text field, type 0.05.

Click .

Inlet 2

Right-click and choose .

Select Boundary 5 only.

Inlet 3

Right-click and choose .

Select Boundary 12 only.

Outlet 1

In the toolbar, click

and choose .

Select Boundary 7 only.

Particle Properties 1

In the window, click .

In the window for , locate the section.

From the ρp list, choose . In the dp text field, type 10[um].

Mesh 1

A reasonably fine mesh is needed on the interface between the stationary and rotating domains.

Edge 1

In the toolbar, click

In the window for , locate the section.

From the list, choose .

Size 1

Right-click and choose .

In the window for , locate the section.

From the list, choose .

Free Triangular 1

In the toolbar, click

Size

In the window, click .

In the window for , locate the section.

From the list, choose .

Click

Click the

button in the toolbar. The mesh should look like Figure 2.

Study 1

Step 1: Time Dependent

In the window, under click .

In the window for , locate the section.

In the text field, type range(0,0.02,2).

From the list, choose .

In the toolbar, click

Results

Particle Trajectories (fpt)

The predefined variable fpt.prf can be used to place colors on a particle based on the inlet where it appeared. This allows you to visualize the effect of the mixing between the three inlets.

Click the

button in the toolbar.

Particle Trajectories 1

In the window, expand the node.

Color Expression 1

In the window, expand the node, then click .

In the window for , click in the upper-right corner of the section. From the menu, choose .

In the toolbar, click

Hide the pair boundary using the option in the node.

Definitions

Hide for Physics 1

In the window, expand the node.

Right-click and choose .

In the window for , locate the section.

From the list, choose .

Results

Particle Trajectories (fpt)

You can reproduce the results in Figure 3 by selecting different values for Time. A better way of visualizing the results is to click the button, in which case the following instructions can be skipped.

In the window, under click .

In the window for , locate the section.

From the list, choose .

In the toolbar, click

Repeat last two steps for the time values 0.4, 0.6, 0.8, 1, and 2 s.