Search in book...
Toggle Font Controls
Create new playlist

Name your new playlist

Playlist description (optional)
Sign In

Email address

Password

Forgot Password?

or

Continue with Facebook

Continue with Google
Sign Up

Full Name

Email address

Confirm Email Address

Password

or

Continue with Facebook

Continue with Google

Chapter 33

Numerical Solutions to Transient Flow Problems

33.1 Introduction

We have already looked at analytical solutions to transient flow problems in section 18. Most of the solutions discussed in this section have been solved by using transient flow solutions, i.e., solutions that are composed of a steady-state solution and a transient solution. We distinguished two types of transient flows.

Time-Dependent Transient Flows. The first type discussed was time-dependent transient flow problems. In these problems, the solution consisted of a time-independent steady-state solution and a time-dependent transient solution (see Eq. 18.2). As we have seen, these problems are commonly used for discussing accelerating and decelerating flow problems.

Entrance Effects. The second type of problem discussed was entrance effects where the solution consisted of a steady-state solution and a space-dependent transient solution that decays as the flow travels (see Eq. 18.64). This type of flow problem occurs whenever the flow experiences a cross-section change and the momentum produced by the new boundary conditions will start diffusing into the fluids thus modifying the flow profiles.

Numerical Approaches. For both examples we were able to derive solutions analytically. We found time-dependent flow problems to be solved conveniently even though we first required the steady-state solution and had to perform a series expansion in order to obtain a solution. Entrance effects were somewhat more difficult to treat analytically since they usually involve the convection term on the left-hand side of the Navier-Stokes equation (see Eq. 11.40), which makes the problem very difficult to solve analytically. We could only derive a solution because we neglected convection, which resulted in a somewhat over-optimistic solution (see section 18.5.5).

As can be imagined, if the flow conditions or the cross-sections become more complex, it becomes more and more difficult to find analytical solutions. Usually, we deviate to numerical methods for solving such problems. We have already used a numerical scheme for solving a two-dimensional Poisson equation in section 28.3.2, which allowed us to derive the flow profiles in rectangular channels conveniently. Obviously, we can extend the function SolvePoissonSOR such that we can also use it to solve transient flow problems. As a matter of fact, this is not very complicated.

Finite Difference Method. In this section, we will extend SolvePoissonSOR to include both time-dependent flow effects and space-dependent effects such as the ones encountered when discussing the entrance flows. As you can see, we are slowly but surely expanding the simple function we wrote initially in a way that it can be used to solve pretty complex flow patterns. The solver uses a scheme that classifies it as a finite difference method solver, i.e., a solver that calculates the finite difference both in space coordinates and time on a finite mesh. In this case, we will limit ourselves to the Cartesian coordinate system, i.e., we will solve the fluid mechanical problems on a two-dimensional mesh using square cells. The edge length of the squares is given by the stepwidth h of our numerical scheme.

Navier-Stokes Equation. If you recall the Poisson equation in two dimensions given by Eq. 16.40 we find

$\frac{\partial^{2} v_{x}}{\partial z^{2}} + \frac{\partial^{2} v_{x}}{\partial y^{2}} = \frac{1}{η} \frac{d p}{d x}$ $\frac{\partial^{2} v_{x}}{\partial z^{2}} + \frac{\partial^{2} v_{x}}{\partial y^{2}} = \frac{1}{η} \frac{d p}{d x}$

si5_e

During derivation, we removed the time-dependent term of the Navier-Stokes equation because we limited the solutions to stationary flows. Obviously, we need to reiterate this derivative because we need the time-dependent term now. However, we still assume the flow to be parallel and driven by a pressure drop along the x-axis only, which is why we can neglect the terms v_y and v_z in which case we can simplify Eq. 11.40 to

$ρ (\frac{\partial v_{x}}{\partial t} + v_{x} \frac{\partial v_{x}}{\partial x}) = k_{x} - \frac{d p}{d x} + η (\frac{\partial^{2} v_{x}}{\partial x^{2}} + \frac{\partial^{2} v_{x}}{\partial y^{2}} + \frac{\partial^{2} v_{x}}{\partial z^{2}})$ $ρ (\frac{\partial v_{x}}{\partial t} + v_{x} \frac{\partial v_{x}}{\partial x}) = k_{x} - \frac{d p}{d x} + η (\frac{\partial^{2} v_{x}}{\partial x^{2}} + \frac{\partial^{2} v_{x}}{\partial y^{2}} + \frac{\partial^{2} v_{x}}{\partial z^{2}})$

si6_e (Eq. 33.1)

where we have left the relevant terms for one dimension (along the x-axis) in the equation. Eq. 33.1 is the unsimplified form of the Navier-Stokes equation for parallel flows, i.e., flows that are described only by the velocity term v_x. Please note that this equation contains no additional simplification, i.e., the Navier-Stokes equation is given entirely at this point.

Burgers’ Equation. You may occasionally come across a version of Eq. 33.1 that is referred to as Burgers’ equation. Burger¹ studied the mathematics of combined convection and diffusion effects that are described by the general formula

$\frac{\partial f}{\partial t} + f \frac{\partial f}{\partial x} = c \frac{\partial^{2} f}{\partial x^{2}}$ $\frac{\partial f}{\partial t} + f \frac{\partial f}{\partial x} = c \frac{\partial^{2} f}{\partial x^{2}}$

where c is a generic diffusion constant, e.g., momentum diffusivity η, mass diffusivity D, or thermal diffusivity α. In case of the Navier-Stokes equation, the dependent variable f is the velocity v_x. Burgers-type equations are found commonly in physics and mathematics: From the study of shock waves to traffic flow, there is a surprisingly large collection of phenomena that can be described by Burgers’ equation.

33.2 A Numerical Solver for Two-dimensional Time-Dependent Flow Problems

33.2.1 Introduction

The first modification we will make to our solver is implementing time-dependency. For this, we use Eq. 33.1. However, for this first step, we will assume the flow to be stationary along the x-axis and not subjected to volume forces. In this case Eq. 33.1 simplifies to

$\begin{array}{l} ρ \frac{\partial v_{x}}{\partial t} = - \frac{d p}{d x} + η (\frac{\partial^{2} v_{x}}{\partial y^{2}} + \frac{\partial^{2} v_{x}}{\partial z^{2}}) \\ \frac{\partial v_{x}}{\partial t} = - \frac{1}{ρ} \frac{d p}{d x} + \frac{η}{ρ} (\frac{\partial^{2} v_{x}}{\partial y^{2}} + \frac{\partial^{2} v_{x}}{\partial z^{2}}) \end{array}$ $\begin{array}{l} ρ \frac{\partial v_{x}}{\partial t} = - \frac{d p}{d x} + η (\frac{\partial^{2} v_{x}}{\partial y^{2}} + \frac{\partial^{2} v_{x}}{\partial z^{2}}) \\ \frac{\partial v_{x}}{\partial t} = - \frac{1}{ρ} \frac{d p}{d x} + \frac{η}{ρ} (\frac{\partial^{2} v_{x}}{\partial y^{2}} + \frac{\partial^{2} v_{x}}{\partial z^{2}}) \end{array}$

si8_e (Eq. 33.2)

In analogy to the implementation of SolvePoissonSOR from listing 28.4 we use the following numerical approximation for the partial differentials on the right-hand side of Eq. 33.2

$\frac{\partial^{2} v_{x}}{\partial y^{2}} + \frac{\partial^{2} v_{x}}{\partial z^{2}} = \frac{1}{h_{y z}^{2}} (v_{x}^{(y + 1, z)} + v_{x}^{(y - 1, z)} + v_{x}^{(y, z + 1)} + v_{x}^{(y, z - 1)} - 4 v_{x}^{(y, z)})$ $\frac{\partial^{2} v_{x}}{\partial y^{2}} + \frac{\partial^{2} v_{x}}{\partial z^{2}} = \frac{1}{h_{y z}^{2}} (v_{x}^{(y + 1, z)} + v_{x}^{(y - 1, z)} + v_{x}^{(y, z + 1)} + v_{x}^{(y, z - 1)} - 4 v_{x}^{(y, z)})$

si9_e (Eq. 33.3)

This is a second-order method. Obviously, we also need an implementation for the time-dependency. For this, we can use, e.g., a forward Euler scheme (see section 27.2.2) according to Eq. 27.8 or a backward Euler scheme (see section 27.2.3) according to Eq. 27.14.

33.2.2 Forward Euler Scheme

The simplest implementation is a forward Euler scheme that is given by

$v_{x}^{(i + 1)} = v_{x}^{(i)} + h_{t} X^{(i)}$ $v_{x}^{(i + 1)} = v_{x}^{(i)} + h_{t} X^{(i)}$

(Eq. 33.4)

The numerical scheme can be assembled from Eq. 33.2, Eq. 33.3, and Eq. 33.4. For this, we use Eq. 33.2 to calculate X⁽ⁱ⁾ in Eq. 33.4. This in terms will give us the value $v_{x}^{(i + 1)}$ $v_{x}^{(i + 1)}$ , i.e., the next value of v_x in time. Our iterative scheme is thus given by

$v_{x}^{(i + 1)} = v_{x}^{(i)} + h_{t} (- \frac{1}{ρ} \frac{d p}{d x} + \frac{η}{ρ h_{y z}^{2}} (v_{x}^{(i, y + 1, z)} + v_{x}^{(i, y - 1, z)} + v_{x}^{(i, y, z + 1)} + v_{x}^{(i, y, z - 1)} - 4 v_{x}^{(i, y, z)}))$ $v_{x}^{(i + 1)} = v_{x}^{(i)} + h_{t} (- \frac{1}{ρ} \frac{d p}{d x} + \frac{η}{ρ h_{y z}^{2}} (v_{x}^{(i, y + 1, z)} + v_{x}^{(i, y - 1, z)} + v_{x}^{(i, y, z + 1)} + v_{x}^{(i, y, z - 1)} - 4 v_{x}^{(i, y, z)}))$

si12_e

$v_{x}^{(t + 1, y, z)} = v_{x}^{(t, y, z)} - \frac{h_{t}}{ρ} \frac{d p}{d x} + \frac{η h_{t}}{ρ h_{y z}^{2}} (v_{x}^{(t, y + 1, z)} + v_{x}^{(t, y - 1, z)} + v_{x}^{(t, y, z + 1)} + v_{x}^{(t, y, z - 1)} - 4 v_{x}^{(t, y, z)})$ $v_{x}^{(t + 1, y, z)} = v_{x}^{(t, y, z)} - \frac{h_{t}}{ρ} \frac{d p}{d x} + \frac{η h_{t}}{ρ h_{y z}^{2}} (v_{x}^{(t, y + 1, z)} + v_{x}^{(t, y - 1, z)} + v_{x}^{(t, y, z + 1)} + v_{x}^{(t, y, z - 1)} - 4 v_{x}^{(t, y, z)})$

si13_e (Eq. 33.5)

As you can see, Eq. 33.5 steps forward in time and we only require values at time step t to calculate the values for time step t + 1. Therefore, the method does not require iteration. Unfortunately, this scheme is a first-order scheme since we used the forward Euler method for the time step. This means that we have to be careful about the size of the steps in time h_t since large values quickly lead to numerical instability. The stepwidth in space h_yz is significantly less problematic since this is a second-order scheme.

33.2.3 Backward Euler Scheme

The second option is using a backward Euler scheme. As discussed in section 27.2.3 the backward Euler method is a good alternative if higher numerical stability is required. The general formula is given by Eq. 27.14 as

$v_{x}^{(i + 1)} = v_{x}^{(i)} + h X^{(i + 1)}$ $v_{x}^{(i + 1)} = v_{x}^{(i)} + h X^{(i + 1)}$

(Eq. 33.6)

If you compare Eq. 33.6 (backward Euler) with Eq. 33.4 (forward Euler) you will see the subtle difference. The backward Euler method requires the gradient at time step i + 1 in order to calculate the value at i + 1. Obviously, the gradient cannot be determined if the value is not known. Therefore, we have to implement an iterative solver. Using Eq. 33.6 to approximate the step in time we can write our iterative scheme Eq. 33.5 to

$\begin{array}{l} v_{x}^{(i + 1)} = v_{x}^{(i)} + h_{t} (- \frac{1}{ρ} \frac{d p}{d x} + \frac{η}{ρ h_{y z}^{2}} (v_{x}^{(i + 1, y + 1, z)} + v_{x}^{(i + 1, y - 1, z)} + v_{x}^{(i + 1, y, z + 1)} + v_{x}^{(i + 1, y, z - 1)} - 4 v_{x}^{(i + 1, y, z)})) \\ v_{x}^{(t + 1, y, z)} = v_{x}^{(t, y, z)} - \frac{h_{t}}{ρ} \frac{d p}{d x} + \frac{η h_{t}}{ρ h_{y z}^{2}} (v_{x}^{(t + 1, y + 1, z)} + v_{x}^{(t + 1, y - 1, z)} + v_{x}^{(t + 1, y, z + 1)} + v_{x}^{(t + 1, y, z - 1)} - 4 v_{x}^{(t + 1, y, z)}) \end{array}$ $\begin{array}{l} v_{x}^{(i + 1)} = v_{x}^{(i)} + h_{t} (- \frac{1}{ρ} \frac{d p}{d x} + \frac{η}{ρ h_{y z}^{2}} (v_{x}^{(i + 1, y + 1, z)} + v_{x}^{(i + 1, y - 1, z)} + v_{x}^{(i + 1, y, z + 1)} + v_{x}^{(i + 1, y, z - 1)} - 4 v_{x}^{(i + 1, y, z)})) \\ v_{x}^{(t + 1, y, z)} = v_{x}^{(t, y, z)} - \frac{h_{t}}{ρ} \frac{d p}{d x} + \frac{η h_{t}}{ρ h_{y z}^{2}} (v_{x}^{(t + 1, y + 1, z)} + v_{x}^{(t + 1, y - 1, z)} + v_{x}^{(t + 1, y, z + 1)} + v_{x}^{(t + 1, y, z - 1)} - 4 v_{x}^{(t + 1, y, z)}) \end{array}$

si15_e (Eq. 33.7)

As we have already discussed, the backward Euler scheme ensures numerical stability. However, due to the fact that it is numerically significantly more expensive (as it requires iteration) we will restrict ourselves with implementing Eq. 33.5, i.e., the forward Euler scheme.

33.2.4 Preparing the Grid

Before implementing the actual function for calculating the velocity profiles we will implement two functions in C for preparing the grid. The first function C_SolveFlow2DPrepareGrid_RectIV will prepare a rectangular grid, and the second function C_SolveFlow2DPrepareGrid_CircIV will prepare a circular grid. As we will see, the function we are about to implement can be used on any geometry. We simply have to correctly set up the boundary and the initial values.

Rectangular Grid. The first function we will implement sets up a rectangular grid (see listing 33.1). The function is called SolveFlow2DPrepareGrid_RectIV and expects the number of time steps T, the width W, the height H, as well as the initial value IV at the first time step. The data to operate on are passed in the form of the pointer data. The pointer BC is expected to contain the boundary values.

Listing 33.1

[C] Listing for preparing a rectangular grid. This listing is part of the NeptunLib and can be found in the file NeptunLib.c in the supplementary material.


 1 void DLLEXPORTIMPORT SolveFlow2DPrepareGrid_RectIV(int T, int W, int H, double IV, double* data, char* BC){.
 2  int currenty=0.
 3  for(int y=0;y<H;y++){.
 4   for(int z=0;z<W;z++){.
 5     BC[currenty+z]=1;
 6     if(y==0|| y==H−1|| z==0|| z==W−1){
 7        data[currenty+z]=0;
 8     }else{
 9       data[currenty+z]=IV;
10     }
11    }
12    currenty+=W;
13   }
14   for(int t=1;t<T;t++){
15     for(int y=0;y<H;y++){
16       for(int z=0;z<W;z++){
17        if(y==0|| y==H−1|| z==0|| z==W−1){
18          data[currenty+z]=0;
19          BC[currenty+z] = 1;
20        }
21       }
22       currenty+=W;
23     }
24   }
25 }

The main variable used is the variable currenty (see line 2). It is used to access the data in both arrays data and BC. The first loop (see line 3) sets up the boundary conditions. These are only required for the first time step, i.e., for T = 0. All of these values are assumed to be boundary values and thus the respective values in the array BC are set to 1 (see line 5). On the boundary of the domain, i.e., for y = 0, y = H − 1, z = 0, and z = W − 1, the velocity, i.e., the values in data, is set to be 0 as they are boundary values (see line 7). For all other values, we assume the given initial value IV (see line 9). As you can see, the two for loops are set up such that the counting variable currenty is incremented by W after finishing with each inner loop (see line 12). This makes the code just a little bit faster as we do not need to recalculate the current position in the arrays every time we pass the inner loop.

The second half of the listing sets up the boundary values for all time points T > 0. Again, for all values y = 0, y = H − 1, z = 0, and z = W − 1, the array BC is set to 1 (see line 19), thus indicating that this is a boundary value in which case we can set the velocity in the array data to zero due to the no-slip boundary condition (see line 18).

Calling the function SolveFlow2DPrepareGrid_RectIV effectively configures the array BC to contain an entry of 1 for all boundary values (at the edge of the grid) and the initial values (all values for T = 0), whereas the array data is set up such that it is set to 0 for all boundary values and to IV for all initial values.

As we will see in just a moment, setting up a grid for a different geometry to be used with our solver only requires a modification of this setup function. The code of the solver is identical.

Circular Grid. We are now going to implement the function SolveFlow2DPrepareGrid_CircIV (see listing 33.2). As you can see, the implementation of this function is almost identical. However, instead of a width and height, this function expects an edge length WH and a radius R. From these it calculates the midpoint mid (see line 3). It then assigns boundary values to all points of the grid that are outside of the circle with radius R around the midpoint mid (see line 7). The rest of the implementation is identical to SolveFlow2DPrepareGrid_RectIV of listing 33.1.

Listing 33.2

[C] Listing for preparing a circular grid. This listing is part of the NeptunLib and can be found in the file NeptunLib.c in the supplementary material.

 1 void DLLEXPORTIMPORT SolveFlow2DPrepareGrid_CircIV(int T, int WH, int R, double IV, double* data, char* BC){
 2  int currenty=0;
 3  double mid=WH/20;
 4  for(int y=0;y<WH;y++){
 5   for(int z=0;z<WH;z++){
 6    BC[currenty+z]=1;
 7    if(sqrt(pow(mid−y,2)+pow(mid−z,2))>=R){
 8     data[currenty+z]=0;
 9    }else{
10     data[currenty+z]=IV;
11    }
12   }
13   currenty+=WH;
14  }
15  for(int t=1;t<T;t++){
16   for(int y=0;y<WH;y++){
17    for(int z=0;z<WH;z++){
18     if(sqrt(pow(mid−y,2)+pow(mid−z,2))>=R){
19        data[currenty+z]=0;
20        BC[currenty+z] = 1;
21      }
22    }
23    currenty+=WH;
24   }
25  }
26 }

33.2.5 Implementation

We are now ready to proceed to the function that implements Eq. 33.5. The function is termed SolveFlow2DTimeDependent and is shown in listing 33.3. It accepts the number of time steps T, the width W, and height H of the grid, as well as the stepwidth in time ht (in μs) and space hyz (in μm). The density density (in gcm^− 3) and the viscosity viscosity (in mPas) as well as the pressure drop dpdx (in mbarmm^− 1) must also be provided. Furthermore, the two arrays data and BC must also be given. They are expected to have been prepared using a function such as SolveFlow2DPrepareGrid_RectIV or SolveFlow2DPrepareGrid_CircIV.

Listing 33.3

[C] Listing for solving time-dependent flow problems in rectangular channel cross-sections. This listing is part of the NeptunLib and can be found in the file NeptunLib.c in the supplementary material.

 1 void DLLEXPORTIMPORT SolveFlow2DTimeDependent(int T, int W, int H, double ht, double hyz, double density, double viscosity, double dpdx, double* data, char* BC){
 2   int currenty=0;
 3   int sizeT=W*H;
 4   double pressurefactor=−01*ht*dpdx/density;
 5   double laplacefactor=ht*viscosity/(density*hyz*hyz);
 6   for(int t=0;t<T−1;t++){
 7     for(int y=0;y<H;y++){
 8       for(int z=0;z<W;z++){
 9        if(BC[currenty+sizeT+z]! = 1){
10          data[currenty+z+sizeT]=data[currenty+z]+pressurefactor+laplacefactor*(data[currenty+z 1]+data[currenty+z+1]+data[currenty+z−W]+data[currenty+z+W] − 4*data[currenty+z]);
11        }
12       }
13       currenty+=W;
14     }
15   }
16   return;

The function first sets up the variable currenty, which is again used as the index variable to access the data in the two arrays. Furthermore, we set up the variable sizeT, which is set to the size of one time frame in the arrays. We can use this value as an offset to jump from one time step to the next in the arrays data and BC (see line 3). The function then calculates two prefactors that we can also find in Eq. 33.5. The prefactor pressurefactor is the prefactor that we find preceding the pressure term in Eq. 33.5. It is corrected with the factor 0.1 to account for the units (see line 4). The second factor laplacefactor precedes the (rewritten) Laplace operator in Eq. 33.5. pressurefactor has the unit mms^− 1, whereas laplacefactor has no unit.

The function then iterates over all time steps. For each time step the inner loops iterate over each cell u_x^(T,y,z) applying Eq. 33.5 to calculate the value of $u_{x}^{(T + 1, y, z)}$ $u_{x}^{(T + 1, y, z)}$ , i.e., the velocity of the next time step at the respective position (see line 10). The value of cell (y, z) is located at the position currenty + z in the arrays data and BC. The subsequent value in z-direction is located at currenty + z + 1 whereas the previous value in z-direction is located at currenty + z − 1. The subsequent value in y-direction is located at currenty + z + W whereas the previous value in y-direction is located at currenty + z − W. In order to access the next time frame of the cell (y, z) the location currenty + z + sizeT must be used.

As you can see, the value is only calculated for cells that are not boundary values, i.e., for which the value in the array BC is not 1 (see line 9). For boundary values, the velocity values are already set in the array data, which is why they are simply skipped. Once all values are processed, the function returns.

33.2.6 Calculating Time-Dependent Transient Flows in Rectangular Channel Profiles

We will now use Flow2DTimeDependent to calculate time-dependent transient flows in rectangular channel profiles. For this we repeat the example we already performed in section 18.4 where we calculated the velocity profiles of an accelerating flow in a rectangular channel cross-section. For this we use the following values that we also used for the example in Fig. 18.9a:

$\begin{array}{l} w = 100 μm \\ h = 50 μm \\ \frac{d p}{d x} = - 0.001 mba {mm}^{- 1} \\ η = 1 mPa s \\ ρ = 1 g {cm}^{- 3} \end{array}$ $\begin{array}{l} w = 100 μm \\ h = 50 μm \\ \frac{d p}{d x} = - 0.001 mba {mm}^{- 1} \\ η = 1 mPa s \\ ρ = 1 g {cm}^{- 3} \end{array}$

si17_e

The time points chosen are t = 10 μs, t = 100 μs, and t = 500 μs, which are also displayed in Fig. 18.9a. The profiles can be calculated using Eq. 18.61. In addition, we will also calculate the profiles for t = 1000 μs simply to verify that the flow will be fully developed and not changing significantly.

In order to properly output our results we implement a Maple worksheet (see listing 33.4). As usual, we first restart the Maple server and input the file Core.txt (see line 1). We then import the functions SolveFlow2DPrepareGrid_RectIV, SolveFlow2DPrepareGrid_CircIV, and SolveFlow2DTimeDependent from NetpunLib.dll, which we implemented in listing 33.1, listing 33.2, and listing 33.3 (see line 2). We pass the arguments as the expected types (integer or float) and announce the size of the arrays data and BC, which are passed by reference in order to allow the C-code to modify the data.

Listing 33.4

[Maple] Listing for solving time–transient flow problems numerically (part 1). This listing is part of the listing TransientFlowNumerical.mw, which can be found in the supplementary material.

 1 restart:read"Core.txt":
 2 C_SolveFlow2DPrepareGrid_RectIV:=define_external('SolveFlow2DPrepareGrid_RectIV',LIB="NeptunLib. dll",T::(integer[4]),W::(integer[4]),H::(integer[4]),IV::(float[8]), 'data' ::(ARRAY(1..W*H*T, float[8])), ' BC' ::(ARRAY(1..W*H*T,boolean[1])));
 3 C_SolveFlow2DPrepareGrid_CircIV:=define_external('SolveFlow2DPrepareGrid_CircIV',LIB="NeptunLib.
dll",T::(integer[4]),WH::(integer[4]),R::(integer[4]),IV::(float[8]), 'data' ::(ARRAY(1..W*H*T, float[8])), ' BC' ::(ARRAY(1..W*H*T,boolean[1])));
 4 C_SolveFlow2DTimeDependent:=define_external('SolveFlow2DTimeDependent',LIB="NeptunLib.dll",T::(integer[4]),W::(integer[4]),H::(integer[4]),ht::(float[8]),hxyz::(float[8]),density::(float[8]),viscosity::(float[8]),dpdx::(float[8]),'data'::(ARRAY(1..W*H*T,float[8])),'BC'::(ARRAY (1..W*H*T, boolean [1])));
 5
 6 SolveFlow2DTimeDependent_Rect:=proc(T,W,H,ht,hyz,density,viscosity,dpdx) local data, BC, z, y:
 7   data:=Array(1..T,1..H,1..W,datatype=float[8],order=C_order):
 8   BC:=Array(1..T,1..H,1..W,datatype=boolean [1],order=C_order):
 9   C_SolveFlow2DPrepareGrid_Rect IV(T,W,H,0,data,BC):
10   C_SolveFlow2DTimeDependent(T,W,H,ht,hyz,density,viscosity,dpdx,data,BC):
11   return data*1000;
12 end proc:
13
14 sol:=SolveFlow2DTimeDependent_Rect(10000,50,25,0.1,2,1,1,−0.001):
15 matrixplot(sol[100]);
16 matrixplot(sol[1000]);
17 matrixplot(sol[5000]);
18 matrixplot(sol[10000]);

In the next step, we define the function SolveFlow2DTimeDependent_Rect in Maple, which simply wraps the call to the respective C-functions and correctly outputs the result (see line 6). The function expects the total number of time steps T, the width W, the height H, the stepwidth in time ht, the stepwidth in space hyz, the density density (in g cm^− 3), the viscosity viscosity (in mPas), and the pressure drop dpdx (in mbar mm^− 1). It then creates the arrays data and BC in the required size. The data array is of type float, and the BC is passed as boolean array (see line 7).

We then call the function C_SolveFlow2DPrepareGrid_RectIV to set up the initial values and boundary conditions whereupon we call the function C_SolveFlow2DTimeDependent to perform the actual calculation (see line 9). The function simply returns the array data multiplied by 1000 after the calculation has completed (see line 11). The reason for this multiplication is due to the fact that we want the output velocity to be displayed in μms^− 1 rather than in mms^− 1. We then use this function to calculate our example for a total of 10000 time steps with a stepwidth of 0.1 μs, which equals a total duration of 1 ms (see line 14). The chosen stepwidth is 2 μm, which is about as large as possible without risking instability of the numerical scheme.

We then extract and plot the results as the time points 100 (10 μs), 1000 (100 μs), 5000 (500 μs), and 10000 (1 ms). The results are shown in Fig. 33.1. If you compare these plots with the analytical solutions of Fig. 18.9a you can see that they match very well. For 10 μs the analytical solution (given by Eq. 18.58) predicts a center velocity of 1.1 18.ms^− 1 where our numerical scheme finds 0.99 μms^− 1. For 100 μs the analytical solution predicts a center velocity of 9.88 μms^− 1 where our numerical scheme finds 9.44 μms^− 1. For the time step 500 μs the analytical solution predicts a velocity of 25.3 μms^− 1 whereas the numerical scheme finds 24.2 μms^− 1. Finally for 1 ms the analytical solution predicts a velocity of 28.56 μms^− 1 whereas the numerical solution finds 26.2 μms^− 1. As you can see, the numerical scheme finds the close-to-correct solution for each time point. You can also see that the change of velocity between 500 μm and 1000 μm is reasonably small and that the profile can be expected to be fully developed at this point.

f33-01-9788535284430 — Fig. 33.1 Numerical results for the accelerating flow in rectangular channel profiles using water as fluid. The channel was 100 μm wide and 50 μm high and the flow accelerated at t = 0 by applying a pressure drop of − 0.001 mbar mm^− 1. The results correspond to the analytical solutions that we obtained analytically for the same geometry in Fig. 18.9a. Please note that the axes denote the number of the cell and not the width or height of the channel.

As you can see, our fairly simple numerical scheme delivers accurate results for the rectangular cross-section. In the next step, we will employ the same solver to solve a different cross-section, i.e., a circular channel cross-section.

33.2.7 Calculating Time-Dependent Transient Flows in Circular Channel Profiles

In the next step, we repeat the example from section 18.3 in which we solved the time-dependent transient flow in circular cross-section channels. Again, we use the same values for the variables in order to be able to compare the results obtained. The values are

$\begin{array}{l} R = 50 μ m \\ \frac{d p}{d x} = - 0.01 m bar {mm}^{- 1} \\ η = 1 m Pas \\ ρ = 1 g {cm}^{- 3} \\ υ_{x, max} = \frac{R^{2}}{4 η} \frac{d p}{d x} = 0.625 {mms}^{- 1} \\ t_{char} = R^{2} \frac{ρ}{η} = 2.5 ms \end{array}$ $\begin{array}{l} R = 50 μ m \\ \frac{d p}{d x} = - 0.01 m bar {mm}^{- 1} \\ η = 1 m Pas \\ ρ = 1 g {cm}^{- 3} \\ υ_{x, max} = \frac{R^{2}}{4 η} \frac{d p}{d x} = 0.625 {mms}^{- 1} \\ t_{char} = R^{2} \frac{ρ}{η} = 2.5 ms \end{array}$

si18_e (Eq. 33.8)

where we have calculated the maximum velocity v_x,_max according to Eq. 16.19 and the characteristic time t_char according to Eq. 18.36. We can now use our implemented function to calculate the velocity profiles during acceleration of the Hagen-Poiseuille flow and compare the results with the dimensionless velocity profiles that we have derived analytically in section 18.3.5. Fig. 18.6a shows these profiles for t = 0, t = 0.02t_char = 50 μs, t = 0.05t_char = 125 μs, t = 0.1t_char = 250 μs, and t = 0.5t_char = 1250 μs.

We then implement the function SolveFlow2DTimeDependent_Circ shown in listing 33.5, which is a slightly modified version of the function SolveFlow2DTimeDependent_Rect implemented in listing 33.4. The only real difference is the fact that we call the function SolveFlow2DTimeDependent_Circ instead of the function C_SolveFlow2DPrepareGrid_RectIV to set up the initial and the boundary values for a circular cross-section. In order to better compare the numerical results we normalize the output data to the maximum velocity given by Eq. 33.8. We then call the function SolveFlow2DTimeDependent_Circ with the given values and request calculation of a total of 12 500 time steps, which corresponds to 1.25 ms (see line 10).

Listing 33.5

[Maple] Listing for solving time-transient flow problems in circular channel cross-sections numerically (part 2). This listing is part of the listing TransientFlowNumerical.mw, which can be found in the supplementary material.

 1 SolveFlow2DTimeDependent_Circ:=proc(T,WH,R,ht,hxyz,density,viscosity,dpdx) local data, BC, z, y:
 2   data:=Array(1..T,1..WH,1..WH,datatype=float[8],order=C_order):
 3   BC:=Array(1..T,1..WH,1..WH,datatype=boolean[1],order=C_order):
 4   C_SolveFlow2DPrepareGrid_CircIV(T,WH,R,0,data,BC):
 5   C_SolveFlow2DTimeDependent(T,WH,WH,ht,hxyz,density,viscosity,dpdx,data,BC):
 6   data:=data/(− 10*R^2/(4*viscosity)*dpdx):
 7   return data;
 8 end proc:
 9
10 solcirc:=SolveFlow2DTimeDependent_Circ(12500,25,10,0.1,5,1,1,-0.01):
11 matrixplot(solcirc[500]);
12 matrixplot(solcirc[1250]);
13 matrixplot(solcirc[2500]);
14 matrixplot(solcirc[12500]);

From the result we plot the normalized velocity profiles at time step 500 (50 μs), 1250 (125 μs), 2500 (250 μs), and 12 500 (1250 μs). The results are shown in Fig. 33.2. If you compare these normalized profiles with the analytical results shown in Fig. 18.6a you can see that they match very well.

f33-02-9788535284430 — Fig. 33.2 Numerical results for the accelerating flow in circular channel profiles using water as fluid. The channel has a radius of 50 μm and the flow is accelerated at t = 0 by applying a pressure drop of − 0.01mbarmm^− 1. The results correspond to the analytical solutions that we obtained analytically for the same geometry in Fig. 18.6a. Please note that the axes denote the number of the cell and not the width or height of of the channel.

For 50 μs, 125 μs, and 250 μs the analytical solution predicts normalized center velocities of 0.08, 0.2, and 0.39. respectively, all of which are also found by the numerical solution. For 1250 μs the analytical solution expects a slightly smaller value (0.95) than the value predicted by the numerical solution (0.99). Overall the results are very close to the correct values.

Overall, it is fairly easy to use the function we implemented with various channel geometries. It is merely required to correctly set up the boundary and initial values differently.

33.2.8 Summary

As you can see, it is reasonably simple to transfer the somewhat more complicated case of time-transient flows to a numerical scheme. Using a forward Euler method we obtain a fast and pretty accurate solver for these types of problems for any type of flow channel geometry. We have demonstrated the solver of listing 33.3 both in rectangular and in circular cross-sections and found that it delivers results that are fully comparable to the results we obtained analytically. This shows how conveniently usable these numerical schemes are. Obviously, there is more than one way to implement the time dependency. We have opted for a forward Euler scheme due to the fast that is numerically cheap and can be implemented conveniently, yielding very fast code. In fact, on most modern computers calculating the 10 000 time steps we used in listing 33.4 will hardly be noticeable.

33.3 A Numerical Solver for Two-Dimensional Entrance Flow Problems

33.3.1 Introduction

In the next problem we will adapt our solver for are entrance flows. We already discussed this case in section 18.5 for Hagen-Poiseuille flow analytically. Entrance flows are the transient regime a flow experiences upon a change of the cross-section. As an example, consider a moving fluid that enters a pipe. The uniform velocity distribution will gradually convert to the characteristic steady-state flow profiles. In case of Hagen-Poiseuille flow, this profile is parabolic (see Eq. 16.20). The length along x at which these entrance effects have vanished and the flow can be considered to be in steady state, is referred to as x_char.

In order to achieve an analytical solution, we had to neglect the convection term on the left-hand side of the Navier-Stokes equation (see Eq. 33.1). We reasoned that neglecting convection would result in a not-quite-correct estimate of the total length x_char (see section 18.5.5). In cases where convection must be considered, i.e., at high Reynolds numbers, the momentum diffusion proceeds significantly slower since it is overlaid by the convective contribution. We failed to proof this analytically, but we will do so in a moment using our numerical solver. In all cases, we assume the flow to be stationary, i.e., the flow changes as it moves along x but it does not change in time. This allows us to drop the partial derivative with respect to time on the left-hand side of Eq. 33.1. In addition we also neglect influences of volume forces in which case Eq. 33.1 simplifies to

$ρ v_{x} \frac{\partial v_{x}}{\partial x} = - \frac{d p}{d x} + η (\frac{\partial^{2} v_{x}}{\partial x^{2}} + \frac{\partial^{2} v_{x}}{\partial y^{2}} + \frac{\partial^{2} v_{x}}{\partial z^{2}})$ $ρ v_{x} \frac{\partial v_{x}}{\partial x} = - \frac{d p}{d x} + η (\frac{\partial^{2} v_{x}}{\partial x^{2}} + \frac{\partial^{2} v_{x}}{\partial y^{2}} + \frac{\partial^{2} v_{x}}{\partial z^{2}})$

si19_e (Eq. 33.9)

33.3.2 Neglecting Convection

33.3.2.1 Simplified Version of the Navier-Stokes Equation

In a first step, we will implement the solver in a way that it yields the same results as previously obtained from our analytical solution. Therefore, we neglect convection, i.e., we drop the left-hand side of Eq. 33.9, which simplifies to

$η (\frac{\partial^{2} v_{x}}{\partial x^{2}} + \frac{\partial^{2} v_{x}}{\partial y^{2}} + \frac{\partial^{2} v_{x}}{\partial z^{2}}) = \frac{d p}{d x}$ $η (\frac{\partial^{2} v_{x}}{\partial x^{2}} + \frac{\partial^{2} v_{x}}{\partial y^{2}} + \frac{\partial^{2} v_{x}}{\partial z^{2}}) = \frac{d p}{d x}$

si20_e (Eq. 33.10)

33.3.2.2 Numerical Scheme

The numerical scheme is the three-dimensional version of Eq. 33.3 that reads

$\frac{\partial^{2} v_{x}}{\partial x^{2}} + \frac{\partial^{2} v_{x}}{\partial y^{2}} + \frac{\partial^{2} v_{x}}{\partial z^{2}} = \frac{1}{h_{x y z}^{2}} (v_{x}^{(x + 1, y, z)} + v_{x}^{(x - 1, y, z)} + v_{x}^{(x, y + 1, z)} + v_{x}^{(x, y - 1, z)} + v_{x}^{(x, y, z + 1)} + v_{x}^{(x, y, z - 1)} - 6 v_{x}^{(x, y, z)})$ $\frac{\partial^{2} v_{x}}{\partial x^{2}} + \frac{\partial^{2} v_{x}}{\partial y^{2}} + \frac{\partial^{2} v_{x}}{\partial z^{2}} = \frac{1}{h_{x y z}^{2}} (v_{x}^{(x + 1, y, z)} + v_{x}^{(x - 1, y, z)} + v_{x}^{(x, y + 1, z)} + v_{x}^{(x, y - 1, z)} + v_{x}^{(x, y, z + 1)} + v_{x}^{(x, y, z - 1)} - 6 v_{x}^{(x, y, z)})$

si21_e (Eq. 33.11)

Obviously, this is the second-order method we have been using so far to approximate the second-order partial derivatives. Inserting Eq. 33.11 into Eq. 33.10 yields

$v_{x}^{(x + 1, y, z)} + v_{x}^{(x - 1, y, z)} + v_{x}^{(x, y + 1, z)} + v_{x}^{(x, y - 1, z)} + v_{x}^{(x, y, z + 1)} + v_{x}^{(x, y, z - 1)} - 6 v_{x}^{(x, y, z)} = \frac{h_{x y z}^{2}}{η} \frac{d p}{d x}$ $v_{x}^{(x + 1, y, z)} + v_{x}^{(x - 1, y, z)} + v_{x}^{(x, y + 1, z)} + v_{x}^{(x, y - 1, z)} + v_{x}^{(x, y, z + 1)} + v_{x}^{(x, y, z - 1)} - 6 v_{x}^{(x, y, z)} = \frac{h_{x y z}^{2}}{η} \frac{d p}{d x}$

si22_e

from which we find

$v_{x}^{(x, y, z)} = - \frac{h_{x y z}^{2}}{6 η} \frac{d p}{d x} + \frac{1}{6} (v_{x}^{(x + 1, y, z)} + v_{x}^{(x - 1, y, z)} + v_{x}^{(x, y + 1, z)} + v_{x}^{(x, y - 1, z)} + v_{x}^{(x, y, z + 1)} + v_{x}^{(x, y, z - 1)})$ $v_{x}^{(x, y, z)} = - \frac{h_{x y z}^{2}}{6 η} \frac{d p}{d x} + \frac{1}{6} (v_{x}^{(x + 1, y, z)} + v_{x}^{(x - 1, y, z)} + v_{x}^{(x, y + 1, z)} + v_{x}^{(x, y - 1, z)} + v_{x}^{(x, y, z + 1)} + v_{x}^{(x, y, z - 1)})$

si23_e (Eq. 33.12)

In comparison to the scheme we used for time-transient flow problems, where we could simply apply a forward Euler scheme, Eq. 33.12 requires recursion since we require $v_{x}^{(x + 1, y)}$ $v_{x}^{(x + 1, y)}$ to calculate v_x^(x,y). This calls for an SOR scheme that we have used previously to solve Poisson equations (see listing 28.4).

33.3.2.3 Implementation

The implementation of Eq. 33.12 will be done in C again because we want to have fairly fast code for these types of problems. It is shown in listing 33.6. As you can see, the function SolveFlow2DEntranceFlowNoConvection accepts the length of the grid X as well as the width W and the height H. Furthermore it requires the stepwidth in space hxyz, the viscosity viscosity (in mPas), the pressure drop dpdx (in mbarmm^− 1), the required precision epsilon, the relaxation parameter ω omega, and the two arrays BC (for the boundary conditions) and data.

Listing 33.6

[C] Listing for solving a space-transient flow problem while neglecting convection. This listing is part of the NeptunLib and can be found in the file NeptunLib.c in the supplementary material.

 1 void DLLEXPORTIMPORT SolveFlow2DEntranceFlowNoConvection(int X, int W, int H, double hxyz, double viscosity, double dpdx, double epsilon, double omega, double* data, char* BC, char* log){
 2  double store;
 3  int count=0;
 4  int currenty=0;
 5  int sizeX=W*H;
 6  char isconverged;
 7  double pressurefactor=−0.1*hxyz*hxyz*dpdx/(6.0*viscosity);
 8  double laplacefactor=1.0/6.0;
 9  double prefactor1=1−omega;
10  double prefactor2=omega;
11  do{
12    currenty=sizeX;
13    isconverged=1;
14    for(int x=1;x<X-1;x++){
15     for(int y=0;y<H;y++){
16      for(int z=0;z<W;z++){
17       if(BC[currenty+z]!=1){
18        store=prefactor1*data[currenty+z]+prefactor2*(pressurefactor+laplacefactor*(data[ currenty+z−sizeX]+data[currenty+z+sizeX]+data[currenty+z−1]+data[currenty+z+1]+data[currenty−W+z]+data[currenty+W+z]));
19        if(isconverged && fabs(store-data[currenty+z])>epsilon){isconverged=0;}
20        data[currenty+z]=store;
21       }
22      }
23      currenty+=W;
24     }
25    }
26    for(int y=0;y<H;y++){
27     for(int z=0;z<W;z++){
28      data[currenty+z]=data[currenty+z−sizeX];
29     }
30     currenty+=W;
31    }
32    count++;
33    if (count==MAXITERATIONS) {
34     sprintf(log, "Reached maximum number of iterations");
35     return;
36    }
37   }while(!isconverged);
38   sprintf(log, "Converged! Number of iterations: %i",count);
39   return;
40 }

As you can see, the implementation is pretty similar to the implementation we used in listing 28.4. Again, in order to quickly step through the array we define the index variable currenty and set sizeX to quickly access the offsets of the individual layers of x in the array (see line 4). We also use the boolean variable isconverged to indicate if we have a converged solution. We then calculate a couple of prefactors: pressurefactor is the prefactor of the pressure term (corrected for the units) in Eq. 33.12, and laplacefactor is the prefactor in front of the (rewritten) Laplace operator in Eq. 33.12 (see line 7). We also defined prefactor1 and prefactor2, which are the two prefactors we require for the relaxation.

We then iterate over the array data starting with the second set (see line 12). As we expect the values for x = 0, i.e., the entry of the domain, to be boundary values we do not have to check these for convergence. We then iterate across the whole domain skipping boundary values. For all values that are not boundary values, we apply Eq. 33.12 (see line 18), i.e., using SOR. Before assigning the new values (which we store in store) we check if the difference between the current and the last value is bigger than the threshold epsilon. If this is the case, we reset the flag isconverged thus forcing another iteration. At the beginning of each iteration we assume the solution to be converged (see line 13) and check if this flag is still 1 at the end of the loop (see line 37). In order to avoid deadlocks, we only iterate for MAXITERATIONS number of iterations (see line 33). One thing to note is the fact that we have no boundary value for x = X, i.e., for the right-hand side of our domain. Therefore we simply assume Neumann boundary conditions here, i.e., a zero-gradient. We can achieve this by simply copying all values from the second to last x to the last x (see line 28).

In the next step, we can now use SolveFlow2DEntranceFlowNoConvection to solve space-transient flow problems. We will begin with the circular cross-section for which we already derived the analytical solutions in section 18.5.

33.3.2.4 Calculating Space-Dependent Transient Flows in Circular Channel Profiles

As before, we will quickly implement a function in Maple that allows us to visualize the results of our function (see listing 33.7). The function C_SolveFlow2DEntranceFlowNoConvection wraps the call to the C-function C_SolveFlow2DEntranceFlowNoConvection. Its implementation is almost identical to the implementation of SolveFlow2DTimeDependent_Circ in listing 33.5. We first call C_SolveFlow2DPrepareGrid_CircIV to prepare the grid with the boundary and initial values for a circular flow profile whereupon we call C_SolveFlow2DEntrance FlowNoConvection to perform the actual calculation. We then normalize the results to the velocity at x = 0 (v0) in order to make the results easier to compare (see line 8). As we have a logging variable in this case, we output log before returning the pointer to the array data.

Listing 33.7

[Maple] Listing for solving space-transient flow problems in circular channel cross-sections numerically neglecting convection (part 3). This listing is part of the listing TransientFlowNumerical.mw, which can be found in the supplementary material.

 1  C_SolveFlow2DEntranceFlowNoConvection:=define_external('SolveFlow2DEntranceFlowNoConvection',LIB=" NeptunLib.dll",X::(integer[4]),W::(integer[4]),H::(integer[4]),hxyz::(float[8]),viscosity::(float[8]),dpdx::(float[8]),epsilon::(float[8]),omega::(float[8]), 'data' ::(ARRAY(1..W*H*X,float [8])),'BC'::(ARRAY(1..W*H*X,boolean[1])),'log'::(string[256]));
 2
 3 SolveFlow2DEntranceFlowNoConvection_Circ:=proc(X,WH,R,hxyz,v0,viscosity,dpdx,epsilon,omega) local data,BC,z,y,log:=" ":
 4   data:=Array(1..X,1..WH,1..WH,datatype=float[8],order=C_order):
 5   BC:=Array(1..X,1..WH,1..WH,datatype=boolean[1],order=C_order):
 6   C_SolveFlow2DPrepareGrid_CircIV(X,WH,R,v0,data,BC):
 7   C_SolveFlow2DEntranceFlowNoConvection(X,WH,WH,hxyz,viscosity,dpdx,epsilon,omega,data,BC,'log'):
 8   data:=data/v0:
 9   print(log);
10   return data;
11 end proc:
12
13 solnoconv_circ:=SolveFlow2DEntranceFlowNoConvection_Circ(200, 25, 10, 5,0.3125,1,−0.01,1E-5,1.8):
14 matrixplot(solnoconv_circ[1]);
15 matrixplot(solnoconv_circ[2]);
16 matrixplot(solnoconv_circ[5]);
17 matrixplot(solnoconv_circ[10]);

We then call the function SolveFlow2DEntranceFlowNoConvection_Circ using the same values we have already used during our example with the accelerating flow in Hagen-Poiseuille flow, which are

$\begin{array}{l} R & = & 50 μm \\ \frac{d p}{d x} & = & - 0.01 mbar m m^{- 1} \\ η & = & 1 m Pa s \\ ρ & = & 1 g c m^{- 3} \\ v_{x, \max} & = & - \frac{R^{2}}{4 η} \frac{d p}{d x} = 0.625 mm s^{- 1} \\ v_{a v} & = & \frac{v_{x, \max}}{2} = 0.3125 mm s^{- 1} \\ t_{char} & = & R^{2} \frac{ρ}{η} = 2.5 ms \end{array}$ $\begin{array}{l} R & = & 50 μm \\ \frac{d p}{d x} & = & - 0.01 mbar m m^{- 1} \\ η & = & 1 m Pa s \\ ρ & = & 1 g c m^{- 3} \\ v_{x, \max} & = & - \frac{R^{2}}{4 η} \frac{d p}{d x} = 0.625 mm s^{- 1} \\ v_{a v} & = & \frac{v_{x, \max}}{2} = 0.3125 mm s^{- 1} \\ t_{char} & = & R^{2} \frac{ρ}{η} = 2.5 ms \end{array}$

si25_e

where we have calculated the average velocity from Eq. 16.35. We then plot the results for x = 1 (0 μm), x = 2 (5 μm), x = 5 (25 μm), and x = 10 (50 μm) in line 13. These are the same values of x that we plotted the analytical solution for in Fig. 18.11. Using a relaxation parameter of ω = 1.8 results in only 48 iterations for this example which is a reasonably fast calculation. The numerical output is shown in Fig. 33.3.

f33-03-9788535284430 — Fig. 33.3 Numerical results for the space-transient flow in circular channel profiles during entry flow using water as fluid. This approach ignores effects caused by convection. The channel has a radius of 50 μm. The flow is expected to be uniform at x = 0 with the Hagen-Poiseuille flow profile gradually forming as x increases. The pressure drop across the channel is − 0.01 mbarmm^− 1. The initial flow velocity is assumed to be the average flow velocity. The results correspond to the analytical solutions that we obtained for the same geometry in Fig. 18.11. Please note that the axes denote the number of the cell and not the width or height of the channel.

As you can see, the results are very close to the analytical results. The analytical solution expects a normalized center velocity of 1.00, 1.26, 1.77, and 1.94 for $\frac{x}{R} = 0, \frac{x}{R} = 0.1, \frac{x}{R} = 0.5$ $\frac{x}{R} = 0, \frac{x}{R} = 0.1, \frac{x}{R} = 0.5$ , and $\frac{x}{R} = 1$ $\frac{x}{R} = 1$ , respectively. The numerical values found are 1.00, 1.28, 1.78, and 2.05, which are obviously close-to-exact.

As you can see, we can conveniently achieve the same solutions we obtained analytically using the numerical scheme of Eq. 33.12. One interesting thing to note is the fact that Eq. 33.12 is a second-order scheme, i.e., it has a decent numerical stability. We can therefore choose a significantly larger stepwidth as we could do for SolveFlow2DTimeDependent, which used a (first-order) forward Euler scheme and therefore lacked numerical stability.

As stated, we were only able to obtain an analytical solution by neglecting convection. In the next step, we will modify our numerical scheme such that convection is accounted for.

33.3.3 Considering Convection

33.3.3.1 Simplified Version of the Navier-Stokes Equation

In the next step, we now want to include convection, which means that the Navier-Stokes equation of Eq. 33.9 is kept as it is. We require a scheme for the partial differential $\frac{\partial v_{x}}{\partial x}$ $\frac{\partial v_{x}}{\partial x}$ and for the second-code partial differentials $\frac{\partial^{2}}{\partial x^{2}}, \frac{\partial^{2} v_{x}}{\partial y^{2}}$ $\frac{\partial^{2}}{\partial x^{2}}, \frac{\partial^{2} v_{x}}{\partial y^{2}}$ and $\frac{\partial^{2} v_{x}}{\partial z^{2}}$ $\frac{\partial^{2} v_{x}}{\partial z^{2}}$ .

First Partial Derivatives. We will use forward distance schemes for all first-order partial derivatives with respect to space (see Eq. 4.7). This means that in general we find

$\frac{\partial v_{x}}{\partial x} \approx \frac{{v_{x}}^{(x + 1, y, z)} - {v_{x}}^{(x, y, z)}}{h}$ $\frac{\partial v_{x}}{\partial x} \approx \frac{{v_{x}}^{(x + 1, y, z)} - {v_{x}}^{(x, y, z)}}{h}$

(Eq. 33.13)

Second Partial Derivatives. Likewise we will use the following schemes for all second-order partial derivatives with respect to space:

$\begin{array}{c} \frac{\partial^{2} v_{x}}{\partial x^{2}} & \approx & \frac{{v_{x}}^{(x + 1, y, z)} + {v_{x}}^{(x - 1, y, z)} - 2 {v_{x}}^{(x, y, z)}}{h^{2}} \\ \frac{\partial^{2} v_{x}}{\partial y^{2}} & \approx & \frac{{v_{x}}^{(x, y + 1, z)} + {v_{x}}^{(x, y - 1, z)} - 2 {v_{x}}^{(x, y, z)}}{h^{2}} \\ \frac{\partial^{2} v_{x}}{\partial z^{2}} & \approx & \frac{{v_{x}}^{(x, y, z + 1)} + {v_{x}}^{(x, y, z - 1)} - 2 {v_{x}}^{(x, y, z)}}{h^{2}} \end{array}$ $\begin{array}{c} \frac{\partial^{2} v_{x}}{\partial x^{2}} & \approx & \frac{{v_{x}}^{(x + 1, y, z)} + {v_{x}}^{(x - 1, y, z)} - 2 {v_{x}}^{(x, y, z)}}{h^{2}} \\ \frac{\partial^{2} v_{x}}{\partial y^{2}} & \approx & \frac{{v_{x}}^{(x, y + 1, z)} + {v_{x}}^{(x, y - 1, z)} - 2 {v_{x}}^{(x, y, z)}}{h^{2}} \\ \frac{\partial^{2} v_{x}}{\partial z^{2}} & \approx & \frac{{v_{x}}^{(x, y, z + 1)} + {v_{x}}^{(x, y, z - 1)} - 2 {v_{x}}^{(x, y, z)}}{h^{2}} \end{array}$

si32_e

33.3.3.2 Numerical Scheme

The numerical scheme we find from Eq. 33.9 when inserted into Eq. 33.14 is

$\begin{array}{l} ρ {v_{x}}^{(x, y, z)} \frac{{v_{x}}^{(x + 1, y, z)} - {v_{x}}^{(x, y, z)}}{h} & = & - \frac{d p}{d x} + \frac{η}{h^{2}} (\begin{array}{l} {v_{x}}^{(x + 1, y, z)} + {v_{x}}^{(x - 1, y, z)} + {v_{x}}^{(x, y + 1, z)} \\ - {v_{x}}^{(x, y - 1, z)} + {v_{x}}^{(x, y, z + 1)} + {v_{x}}^{(x, y, z - 1)} - 6 {v_{x}}^{(x, y, z)} \end{array}) \\ {v_{x}}^{(x, y, z)} (\frac{{v_{x}}^{(x + 1, y, z)} - {v_{x}}^{(x, y, z)}}{h} + \frac{6 η}{ρ h^{2}}) & = & - \frac{1}{ρ} \frac{d p}{d x} + \frac{η}{ρ h^{2}} (\begin{array}{l} {v_{x}}^{(x + 1, y, z)} + {v_{x}}^{(x - 1, y, z)} + {v_{x}}^{(x, y + 1, z)} \\ + {v_{x}}^{(x, y - 1, z)} + {v_{x}}^{(x, y, z + 1)} - {v_{x}}^{(x, y, z - 1)} \end{array}) \end{array}$ $\begin{array}{l} ρ {v_{x}}^{(x, y, z)} \frac{{v_{x}}^{(x + 1, y, z)} - {v_{x}}^{(x, y, z)}}{h} & = & - \frac{d p}{d x} + \frac{η}{h^{2}} (\begin{array}{l} {v_{x}}^{(x + 1, y, z)} + {v_{x}}^{(x - 1, y, z)} + {v_{x}}^{(x, y + 1, z)} \\ - {v_{x}}^{(x, y - 1, z)} + {v_{x}}^{(x, y, z + 1)} + {v_{x}}^{(x, y, z - 1)} - 6 {v_{x}}^{(x, y, z)} \end{array}) \\ {v_{x}}^{(x, y, z)} (\frac{{v_{x}}^{(x + 1, y, z)} - {v_{x}}^{(x, y, z)}}{h} + \frac{6 η}{ρ h^{2}}) & = & - \frac{1}{ρ} \frac{d p}{d x} + \frac{η}{ρ h^{2}} (\begin{array}{l} {v_{x}}^{(x + 1, y, z)} + {v_{x}}^{(x - 1, y, z)} + {v_{x}}^{(x, y + 1, z)} \\ + {v_{x}}^{(x, y - 1, z)} + {v_{x}}^{(x, y, z + 1)} - {v_{x}}^{(x, y, z - 1)} \end{array}) \end{array}$

si33_e (Eq. 33.14)

If you look at Eq. 33.14 you will see that the left-hand side of this equation is nonlinear since ${v_{x}}^{(x, y, z)}$ ${v_{x}}^{(x, y, z)}$ occurs both before and within the bracket. However, as we need to implement an iterative scheme (due to the fact that we again have to calculate ${v_{x}}^{(x, y, z)}$ ${v_{x}}^{(x, y, z)}$ without knowing ${v_{x}}^{(x + 1, y, z)}$ ${v_{x}}^{(x + 1, y, z)}$ we can as well extend the iteration to the left-hand side and use a previous value of ${v_{x}}^{(i, x, y, z)}$ ${v_{x}}^{(i, x, y, z)}$ to calculate ${v_{x}}^{(i + 1 x, y, z)}$ ${v_{x}}^{(i + 1 x, y, z)}$ . We therefore rewrite the scheme to

$\begin{array}{l} {v_{x}}^{(i + 1, x, y, z)} (ρ h ({v_{x}}^{(i, x + 1, y, z)} - {v_{x}}^{(i, x, y, z)}) + 6 η) & = & - h^{2} \frac{d p}{d x} + η (\begin{array}{l} {v_{x}}^{(i, x + 1, y, z)} + {v_{x}}^{(i, x - 1, y, z)} + {v_{x}}^{(i, x, y + 1, z)} \\ + {v_{x}}^{(i, x, y - 1, z)} + {v_{x}}^{(i, x, y, z + 1)} - {v_{x}}^{(i, x, y, z - 1)} \end{array}) \\ {v_{x}}^{(i + 1, x, y, z)} ({v_{x}}^{(i, x + 1, y, z)} - {v_{x}}^{(i, x, y, z)} + 6 \frac{η}{ρ h}) & = & - \frac{h^{2}}{ρ h} \frac{d p}{d x} + \frac{η}{ρ h} (\begin{array}{l} {v_{x}}^{(i, x + 1, y, z)} + {v_{x}}^{(i, x - 1, y, z)} + {v_{x}}^{(i, x, y + 1, z)} \\ + {v_{x}}^{(i, x, y - 1, z)} + {v_{x}}^{(i, x, y, z + 1)} - {v_{x}}^{(i, x, y, z - 1)} \end{array}) \end{array}$ $\begin{array}{l} {v_{x}}^{(i + 1, x, y, z)} (ρ h ({v_{x}}^{(i, x + 1, y, z)} - {v_{x}}^{(i, x, y, z)}) + 6 η) & = & - h^{2} \frac{d p}{d x} + η (\begin{array}{l} {v_{x}}^{(i, x + 1, y, z)} + {v_{x}}^{(i, x - 1, y, z)} + {v_{x}}^{(i, x, y + 1, z)} \\ + {v_{x}}^{(i, x, y - 1, z)} + {v_{x}}^{(i, x, y, z + 1)} - {v_{x}}^{(i, x, y, z - 1)} \end{array}) \\ {v_{x}}^{(i + 1, x, y, z)} ({v_{x}}^{(i, x + 1, y, z)} - {v_{x}}^{(i, x, y, z)} + 6 \frac{η}{ρ h}) & = & - \frac{h^{2}}{ρ h} \frac{d p}{d x} + \frac{η}{ρ h} (\begin{array}{l} {v_{x}}^{(i, x + 1, y, z)} + {v_{x}}^{(i, x - 1, y, z)} + {v_{x}}^{(i, x, y + 1, z)} \\ + {v_{x}}^{(i, x, y - 1, z)} + {v_{x}}^{(i, x, y, z + 1)} - {v_{x}}^{(i, x, y, z - 1)} \end{array}) \end{array}$

si39_e

from which we find

${v_{x}}^{(i + 1, x, y, z)} = \frac{1}{{v_{x}}^{(i, x + 1, y, z)} - {v_{x}}^{(i, x, y, z)} + 6 \frac{η}{ρ h}} (- \frac{h}{ρ} \frac{d p}{d x} + \frac{η}{ρ h} (\begin{array}{l} {v_{x}}^{(i, x + 1, y, z)} + {v_{x}}^{(i, x - 1, y, z)} + {v_{x}}^{(i, x, y + 1, z)} \\ + {v_{x}}^{(i, x, y - 1, z)} + {v_{x}}^{(i, x, y, z + 1)} - {v_{x}}^{(i, x, y, z - 1)} \end{array}))$ ${v_{x}}^{(i + 1, x, y, z)} = \frac{1}{{v_{x}}^{(i, x + 1, y, z)} - {v_{x}}^{(i, x, y, z)} + 6 \frac{η}{ρ h}} (- \frac{h}{ρ} \frac{d p}{d x} + \frac{η}{ρ h} (\begin{array}{l} {v_{x}}^{(i, x + 1, y, z)} + {v_{x}}^{(i, x - 1, y, z)} + {v_{x}}^{(i, x, y + 1, z)} \\ + {v_{x}}^{(i, x, y - 1, z)} + {v_{x}}^{(i, x, y, z + 1)} - {v_{x}}^{(i, x, y, z - 1)} \end{array}))$

si40_e (Eq. 33.15)

Eq. 33.15 is an iterative scheme that works around the problem of the nonlinearity of the convection term. It uses the velocity values ${v_{x}}^{(i)}$ ${v_{x}}^{(i)}$ at i to calculate a better approximation ${v_{x}}^{(i + 1)}$ ${v_{x}}^{(i + 1)}$ .

A Note to the Choice of Scheme for the First Derivative. Obviously, it is important which of the two velocity terms on the left-hand side of Eq. 33.14 is chosen to be the actual variable at step i + 1. This will be the value we will solve for iteration i + 1. In contrast the second value can be substituted with the value of v_xof iteration i. Here, we chose the first derivative, i.e., $\frac{\partial v_{x}}{\partial x}$ $\frac{\partial v_{x}}{\partial x}$ to be known. This is the common approach of doing it. We also could have chosen a different numerical scheme for this derivation, e.g., the central finite difference given by Eq. 4.6. This does not even include the value ${v_{x}}^{(x, y, z)}$ ${v_{x}}^{(x, y, z)}$ but only the values ${v_{x}}^{(x - 1, y, z)}$ ${v_{x}}^{(x - 1, y, z)}$ and ${v_{x}}^{(x + 1, y, z)}$ ${v_{x}}^{(x + 1, y, z)}$ . In this case, we work around this problem of having the value ${v_{x}}^{(x, y, z)}$ ${v_{x}}^{(x, y, z)}$ twice.

Given the fact that the choice of numeric scheme for the first derivative is more or less random, we consider the prefactor of $v_{x} \frac{\partial v_{x}}{\partial x}$ $v_{x} \frac{\partial v_{x}}{\partial x}$ to be the variable and the derivative to be calculated from the known values of the last iteration, if required.

33.3.3.3 Implementation

Before implementing the scheme in C we will quickly write out the prefactors and correct them for the units. We assume all velocity values to be given in mm s^− 1, the stepwidth h in μm, the density ρ in g cm^− 3, the viscosity η in mPa s, and the pressure drop in mbar mm^− 1. The unit-corrected prefactors are

$\begin{array}{l} pressure factor & = & - 100 \frac{h}{ρ} \frac{d p}{d x} [{mm}^{2} s^{- 2}] \\ laplace factor & = & 1000 \frac{η}{ρ h} [{mm s}^{- 1}] \end{array}$ $\begin{array}{l} pressure factor & = & - 100 \frac{h}{ρ} \frac{d p}{d x} [{mm}^{2} s^{- 2}] \\ laplace factor & = & 1000 \frac{η}{ρ h} [{mm s}^{- 1}] \end{array}$

si49_e

The implementation of Eq. 33.15 is SolveFlow2DEntranceFlowWithConvection shown in listing 33.8. As you can see, the listing is almost identical to the implementation of SolveFlow2DEntranceFlowNoConvection from listing 33.6. The main difference is line 18, which implements Eq. 33.15, i.e., the equation with convection.

Listing 33.8

[C] Listing for solving a space-transient flow problem taking convection into account. This listing is part of the NeptunLib and can be found in the file NeptunLib.c in the supplementary material.

 1 void SolveFlow2DEntranceFlowWithConvection(int T, int W, int H, double hxyz, double density, double viscosity, double dpdx, double epsilon, double omega, double* data, char* BC, char* log ){
 2   double store;
 3   int count=0;
 4   int currenty=0;
 5   int sizeX=W*H;
 6   char isconverged;
 7   double pressurefactor=−100*hxyz*dpdx/density;
 8   double laplacefactor=1000*viscosity/(density*hxyz);
 9   double prefactor1=1−omega;
10   double prefactor2=omega;
11   do{
12     currenty=sizeX;
13     isconverged=1;
14     for(int t=1;t<T−1;t++){
15      for(int y=0;y<H;y++){
16       for(int z=0;z<W;z++){
17        if(BC[currenty+z]!=1){
18          store=prefactor1*data[currenty+z]+prefactor2*(pressurefactor+laplacefactor*(data[ currenty+z−sizeX]+data[currenty+z+sizeX]+data[currenty+z−1]+data[currenty+z+1]+data[currenty−W + z]+data[currenty+W+z]))/(6*laplacefactor+data[currenty+z+sizeX] − data[currenty+z]);
19          if(isconverged && fabs(store-data[currenty+z])>epsilon){isconverged=0;}
20          data[currenty+z]=store;
21        }
22       }
23       currenty+=W;
24      }
25     }
26     for(int y=0;y<H;y++){
27      for(int z=0;z<W;z++){
28       data[currenty+z]=data[currenty+z−sizeX];
29      }
30      currenty+=W;
31     }
32     count++;
33     if (count==MAXITERATIONS) {
34      sprintf(log, "Reached maximum number of iterations");
35      return;
36     }
37    }while(!isconverged);
38    sprintf(log, "Converged! Number of iterations: %i", count);
39    return;
40 }

33.3.3.4 Calculating Space-Dependent Transient Flows in Circular Channel Profiles

Having implemented SolveFlow2DEntranceFlowWithConvection we are now able to compare the two cases of including and neglecting convection side by side. However, before doing these comparisons we first want to rerun our example using listing 33.7, taking convection into account this time. For this we implement the function SolveFlow2DEntranceFlowWithConvection_Circ in Maple, which is shown in listing 33.9. The implementation is identical to the implementation of SolveFlow2DEntranceFlowNoConvection_Circ from listing 33.7 with the sole exception of calling the C-function C_SolveFlow2DEntranceFlowWithConvection instead of C_SolveFlow2DEntranceFlowNoConvection (see line 7). We then run the same example as before (see line 13). Again, we plot the calculated velocity profiles at $\frac{x}{R} = 0, \frac{x}{R} = 0.1, \frac{x}{R} = 0.5, and \frac{x}{R} = 1$ $\frac{x}{R} = 0, \frac{x}{R} = 0.1, \frac{x}{R} = 0.5, and \frac{x}{R} = 1$ (see Fig. 33.4). As you can see, the output looks exactly like the plots in Fig. 33.3, which were calculated for the same example but without taking convection into account. This shows that convection does not play a significant role in this example and can therefore be ignored. As discussed the dimensionless number that allows us to state whether or not we are allowed to ignore convection is the Reynolds number (see section 9.9.8). In our example it is

Listing 33.9

[Maple] Listing for solving space-transient flow problems in circular channel cross-sections numerically taking convection into account (part 4). This listing is part of the listing TransientFlowNumerical.mw that can be found in the supplementary material.

 1  C_SolveFlow2DEntranceFlowWithConvection:=define_external('SolveFlow2DEntranceFlowWithConvection',
LIB="NeptunLib.dll",X::(integer[4]),W::(integer[4]),H::(integer[4]),h::(float[8]),density::( float[8]),visosity::(float[8]),dpdx::(float[8]),epsilon::(float[8]),omega::(float[8]), ' data'::(ARRAY(1..W*H*X,float[8])),'BC'::(ARRAY(1..W*H*X,boolean[1])),'log'::(string[256]));
 2
 3 SolveFlow2DEntranceFlowWithConvection_Circ:=proc(X,WH,R,hxyz,v0,density,viscosity,dpdx,epsilon,omega) local data,BC,z,y,log:="":
 4   data:=Array(1..X,1..WH,1..WH,datatype=float[8],order=C_order):
 5   BC:=Array(1..X,1..WH,1..WH,datatype=boolean[1],order=C_order):
 6   C_SolveFlow2DPrepareGrid_CircIV(X,WH,R,v0,data,BC):
 7   C_SolveFlow2DEntranceFlowWithConvection(X,WH,WH,hxyz,density,viscosity,dpdx,epsilon,omega,data,BC,'log'):
 8   data:=data/v0:
 9   print(log);
10   return data;
11 end proc:
12
13 solwithconv_circ:=SolveFlow2DEntranceFlowWithConvection_Circ(200,25,10,5,0.3125,1, − 0.01,1E-5,1.8):
14 matrixplot(solwithconv_circ[1]);
15 matrixplot(solwithconv_circ[2]);
16 matrixplot(solwithconv_circ[5]);
17 matrixplot(solwithconv_circ[10]);

f33-04-9788535284430 — Fig. 33.4 Numerical results for the space-transient flow in circular channel profiles during entry flow using water as fluid. This approach takes convection into account. The channel has a radius of 50 μm. The flow is expected to be uniform at x = 0 with the Hagen-Poiseuille flow profile gradually forming as x increases. The pressure drop across the channel is − 0.01 mbar mm^− 1. The initial flow velocity is assumed to be the average flow velocity. The results are identical to the results shown in Fig. 33.3 where the same example was used but without taking convection into account. This suggests that convection plays no significant role in this example, which is to be expected due to the low Reynolds number Re = 0.03125. Please note that the axes denote the number of the cell and not the width or height of of the channel.

$\begin{array}{l} R e & = & \frac{ρ v L_{char}}{η} \\ = & \frac{ρ v d}{η} = 0.03125 ≪ 1 \end{array}$ $\begin{array}{l} R e & = & \frac{ρ v L_{char}}{η} \\ = & \frac{ρ v d}{η} = 0.03125 ≪ 1 \end{array}$

si51_e

which states that we are perfectly allowed to ignore effects of convection. Please recall that during the analytical derivation of the flow profiles in section 15 we have made this simplification, without which analytical solutions are almost impossible to make.

33.3.3.5 An Example at Higher Reynolds Numbers

Obviously, we require an example with higher Reynolds numbers in order to effectively see the difference convection can make in such cases. For this we take the same geometry but increase the pressure drop substantially. If we assume a value of − 10 mbar mm^− 1 we end up with an average velocity of 312.55 mm s^− 1 and a Reynolds number of 32.25. This value shows us that convection cannot be neglected and we should see a difference in the outputs of our two functions SolveFlow2DEntranceFlowNoConvection_Circ and C_SolveFlow2DEntranceFlowWithConvection.

The Maple listing we will use is shown in listing 33.10. It calculates the given example with and without taking convection into account and compares the velocity profiles obtained in both cases. As discussed in section 18.5.5 the momentum transport due to convection along the channel axis will overlay the momentum diffusion along the radius. Momentum diffusing from the boundary of the channel will therefore be transported continuously out of the control volume, and we therefore expect the velocity profiles to take significantly longer to establish fully.

Listing 33.10

[Maple] Listing for solving space-transient flow problems in circular channel cross-sections numerically. This example studies the influence of convection using a case at high Reynolds numbers (part 5). This listing is part of the listing TransientFlowNumerical.mw, which can be found in the supplementary material.

1 solnoconv_circ:=SolveFlow2DEntranceFlowNoConvection_Circ(500,105,50,1,312.55,1,−10,1E−5,1.8):
2 solwithconv_circ:=SolveFlow2DEntranceFlowWithConvection_Circ(500,105,50,1,312.55,1,1, − 10,1E−5, 1.8) :
3
4 matrixplot(solnoconv_circ[10]);
5 matrixplot(solwithconv_circ[10]);
6 matrixplot(solnoconv_circ[100]);
7 matrixplot(solwithconv_circ[100]);
8 matrixplot(solnoconv_circ[500]);
9 matrixplot(solwithconv_circ[500]);

The Maple code we need for this example is shown in listing 33.10. As you can see, we decrease the stepwidth to 1 μm, which ensures numerical stability of SolveFlow2DEntranceFlowWithConvection_Circ, which tends not to be able to find a solution if the stepwidth is too large. This is indicated by the log message “Reached maximum number of iterations” from our C-code, which tells us that even after the given maximum number of iterations, the solution has not converged. The output of listing 33.10 is shown in Fig. 33.5 and Fig. 33.6. This example clearly shows the significance of convection in the given example. We know that the fully developed velocity profile will have a normalized center velocity of around 2. For the case in which we ignore convection, this case is achieved at X = 100 ( $\frac{x}{R} > 2$ $\frac{x}{R} > 2$ ). However, if we take convection into account, the flow profile will only have reached an normalized center velocity of around 1.59 at this point. The fully developed profile is only achieved for a value around X = 500 ( $\frac{x}{R} > 10$ $\frac{x}{R} > 10$ ), where it has reached a value of 1.97.

f33-05-9788535284430 — Fig. 33.5 Comparison of numerical results for the space-transient flow in circular channel profiles during entry flow studying the effects of convection using water as fluid (part 1). The channel has a radius of 50 μm. The flow is expected to be uniform at x = 0 with the Hagen-Poiseuille flow profile gradually forming as the x increases. The pressure drop across the channel is − 10mbarmm^− 1. The Reynolds number of this problem is 32.25. The initial flow velocity is assumed to be the average flow velocity. As you can see, convection plays a significant effect. The momentum in- and outflux along the channel axis leads to a significantly reduced momentum diffusion in radial direction. As a consequence, the fully developed flow profiles are only obtained for $\frac{x}{R} > 10$ $\frac{x}{R} > 10$ , which is significantly longer than the length we found in the case where we neglected convection. In the case where convection is neglected, the profile is fully established around $\frac{x}{R} > 2$ $\frac{x}{R} > 2$ . Please note that the axes denote the number of the cell and not the width or height of of the channel.

f33-06-9788535284430 — Fig. 33.6 Comparison of numerical results for the space-transient flow in circular channel profiles during entry flow studying the effects of convection using water as fluid (part 2). The channel has a radius of 50 μm. The flow is expected to be uniform at x = 0 with the Hagen-Poiseuille flow profile gradually forming as the x increases. The pressure drop across the channel is − 10mbarmm^− 1. The Reynolds number of this problem is 32.25. The initial flow velocity is assumed to be the average flow velocity. As you can see, convection plays a significant effect. The momentum in- and outflux along the channel axis leads to a significantly reduced momentum diffusion in radial direction. As a consequence, the fully developed flow profiles are only obtained for $\frac{x}{R} > 10$ $\frac{x}{R} > 10$ , which is significantly longer than the length we found in the case where we neglected convection. In the case where convection is neglected, the profile is fully established around $\frac{x}{R} > 2$ $\frac{x}{R} > 2$ . Please note that the axes denote the number of the cell and not the width or height of of the channel.

This finding is along the lines of what we expect. The momentum diffusing from the boundary of the channel is transported along the axis by the momentum influx due to convection. As a consequence, it takes diffusion significantly longer to propagate in the radial direction, thereby establishing a fully developed flow profile. This can be conveniently seen using the numerical solvers we just implemented.

33.4 Summary

In this chapter, we have studied numerical methods to solve transient-flow profiles both transient in space and time. For both cases, we implemented numerical solvers that we found to yield results that compare well to the analytical solutions we derived for these cases. In addition, implementing a numerical solver for space-transient problems allowed us, for the first time, to study the influence of convection in a Hagen-Poiseuille flow. Although we have only studied these effects on one example, it is very easy to modify the Maple worksheets in a way that arbitrary flow problems can be solved. The way we set up the solvers allows us to reuse the functions for arbitrary cross-sections, which yields very generic and usable solvers for a wide variety of flow problems. Although this chapter describes two-dimensional flow profiles, the solvers we wrote essentially are still one-dimensional since they only need to solve one of the three equations of the Navier-Stokes equation. In the following chapter, we will extend this code to account for fully three-dimensional flow problems.

..................Content has been hidden....................

You can't read the all page of ebook, please click here login for view all page.

Table of Contents for Chapter 33: Numerical Solutions to Transient Flow Problems

Create new playlist

Sign In

Sign Up