COMPUTATIONAL ANALYSIS OF GAS SOLID FLOW IN A PNEUMATIC DRYER

Master's Degree Thesis ISRN: BTH-AMT-EX--2008/D-04--SE

Supervisors: Sharon Kao-Walter, Ph.D Mech Eng, BTH Junruo Chen, Professor, KUST, Kunming, China

Department of Mechanical Engineering Blekinge Institute of Technology

Karlskrona, Sweden 2008

Gao Feng

Adelere Oluwafemi. A.

Computation and Analysis of Gas-Solid Flow Multiphysics

in a Pneumatic Dryer

Computation and

Analysis of Gas-Solid Flow Multiphysics in a Pneumatic

Dryer

Gao Feng

Adelere Oluwafemi. A.

Department of Mechanical Engineering Blekinge Institute of Technology

Karskrona, Sweden 2008

Thesis submitted for completion of Master of Science in Mechanical Engineering with emphasis on Structural Mechanics at the Department of Mechanical Engineering, Blekinge Institute of Technology, Karlskrona, Sweden.

Abstract:

This work presents a design of drying small food particles like tea, tobacco with a reduced height of about fifty percent compared to the convectional straight type. The geometry is designed such that, the particle leave the dryer at an average constant time, to prevent variation in the dryness properties,

Simulation result of straight type dryer in FLUENT, the commercial software is validated by industrial data implemented in MATLAB, which shows reasonable agreement. A reduced height type with a swirling or screw type is designed to reduce the height and also maintain the drying residence time.

Keywords:

Finite volume method ( FVM ) , SIMPLE method, Straight pipe dryer, Screw type dryer, MATLAB^®, GAMBIT^®, FLUENT^®

2

Acknowledgements

This work was carried out under the cooperation between the Department of Mechanical Engineering, Blekinge Institute of Technology and the Faculty of Electro-Mechanical Engineering, Kunming University of Science and Technology China under the supervision of Professor Junruo Chen and Dr. Sharon Kao-Walter

We would like to express sincere gratitude to our supervisors Professor Junruo Chen and Dr.Sharon Kao-Walter, who made this collaborative research effort. The Lineaus-Palme International exchange between Blekinge Institute of Technology Sweden and Kunming University of Science and Technology China an incredible reality.

We are also indebted to Mr. Xianxi Liu (PhD. student) for his valuable suggestions, and all members of staffs in Mechanical Engineering Departments of Blekinge Institute of Technology (BTH), Sweden and Kunming University of Science and Technology (KUST), China. For their guidance and professional engagement on the course of our study.

Kunming, March 2008.

Gao Feng.

Adelere Oluwafemi. A.

3

Contents

1 Notation 7

2 Introduction 9

2.1 Background 9

2.2 Project Description 9

2.3Method 10

2.4 Limitation 10

3 Overview of Drying Process 11

3.1 Drying Process 12

3.2 Available Convective types 12

3.3 Existing Work 13

4 Theoretical Analysis 14

4.1 Derivations of Equations 14

5 Available Experimental Data 19

6 Numerical Analysis 20

6.1 Control-Volume Formulation (F V M) or (FEM) 20

6.2 The SIMPLE algorithm 21

6.3 Implementations Of Boundary Conditions 29

6.4 Drying In a Straight Pipe Model 30

6.5 Design and Working principle of the Swirling Dryer 32

7 Simulation 34

7.1 Choice of Commercial Softwear 34

7.2 Problem Formulation in Commercial Softwear 34

7.3 Boundary Condition 34

7.4 Assumptions 35

7.5 Postprocessing 35

8 Result And Analysis 36

8.1 Single Phase 36

8.2 Double Phase 40

8.3 Swirling type double phase 43

9 Conclusion and Future Work 46

9.1 Conclusion 46

9.2 Future Work 46

10 References 47

11 Appendices 49

A. Other Work In FLUENT 49

B. MATLAB Code 53

4

1 Notations

c specific heat of the body

k Thermal conductivity of the body

T Temperature

A Area

t time

F Force

τ Shear Force

d distance

L Length

u Velocity of x component v Velocity of y component

p pressure

R Gas constant

Г Diffusivity p’ Pressure correction

φ Generalized variables

V Volume

S Source term

ρ Density of the body α Relaxation factor

b Tolerance

i,j Referring to nodes i,j

a Coefficient of SIMPLE method equation Min Mass flux at the inlet

Mout Mass flux at the outlet

5

new New value after correction * Assumption value

e East direction w West direction

s South direction n North direction

6

2 Introduction

2.1 Background

Gas flow drying is an approach of convective desiccation. Hot gas convey the material to be dried into the drying equipment; As a result of the high temperature difference, convective heat transfer takes place between the hot gas and material during the transportation, thereby, removing moisture, inform of vapor from the material. The evaporated moisture from the product becomes part of the drying medium, and does not need to be exhausted, since the hot gas is superheated steam, unless the pressure increases beyond a set point, at which the excess steam may be released.

This allows for the recycling of the drying medium,

This approach is widely used in food, tobacco and some chemical, medical, fertilizer, dyestuff, inorganic compound, industry.

Since the velocity of the airflow is quite high, the material granule is in suspending status, the contact area between the gas-solid two phases is large, so the heat transfer is intensified. The material drying time depends on the temperature of the hot gas, diffusibility of the material, initial and final moisture content. The temperature of the hot gas can adjust to control the temperature of the material and avoid quality change as result of overheating.

2.2 Project Description

However, it is not easy to design a proper drying system, factors mentioned earlier has to be taken into consideration. Such as the dryness quality, moisture content of the product. The disadvantages of the gas-flow drying;

are the resistance of system flow is high and the effective distance of drying is rather long. The straight pipe type dry system commonly used for material drying is designed at a height of about 20 meters,which is too high to build the workshop. So ,the need to consider a new approach to lower the equipment height and also maintain the drying distance arise, with the design the screw type dryer system and make some numerical calculations to simulate the model for optimum design.

7

2.3 Method

Since the process of drying involves movement of particle from one point to another by hot gas, it is a mixture flow of solid (material) and gas (hot gas).

There are research works on solid-gas, two-phase flows, and there are detail results, nevertheless, there are slight variations between the numerical simulation and the real situation due to the complexity and the limitation.

An exploring investigation carried out in this paper is as follows:

• Finding the suitable mathematic model of coupled physic field;

• Simplifying the model and making the reasonable assumptions;

• Using numerical method to solve and implement in MATLAB;

• Validating the results from MATLAB by commercial software;

• Verifying the results from experimental test (optional);

Numerical simulation technology on the computer, to analysis the flow, and the heat transfer process in the project in question is an extremely economical and effective method, which save time and resources. The research results provided the theoretical basis for optimum product development and design, simultaneously, the method has provided useful means for the solid-gas two-phase flows research, but also built the foundation for the practical application.

2.4 Limitations

Effect of diffusibility, and porosity of the material media, compartment of the materials is not considered, only the drying residence time or drying time of the material is considered. Diffusion coefficient is a material property which is affected by the composition, shape, size of the product and process temperature.

Also the mechanisms of separating the dry product from the superheated steam are not included in this work.

8

3 Overview of Drying Process

3.1 Drying Process

The primary objective of food drying is preservation of the product, and lowering the moisture content, which helps to prevent, or reduce microbial and enzymatic reactions; however, drying can have an adverse effect on chemical, physical, and nutritional value of food products [1]. Superheated steam is generally used gas for drying food substance to prevent contaminations as a result of chemical reactions; and the temperature is moderately high to prevent saturation of the drying gas and burning of the food particles as result of heat transfer.

Goals of research is improving current drying technologies, and to improve on economics of operations ,that is reduce energy consumption, increase capacity, reduce size of equipment, increase ease of control, environmental considerations, reduce emissions, increase safety, and improve product quality, minimize chemical, physical, and nutritional degradation. Some of these goals must be proven to provide benefits that exceed additional costs and risks before industry is willing to adopt a new technology. [2]

Drying is process of removing moisture from a material, this can be achieved or done by various methods and technologies, such as chemical decomposition of water in substance, that is, desiccants [3], which could be in form of absorption in gases and capillary action in solid [4].Compression, centrifugal and gravitational forces are employed as mechanical method of drying .Thermal drying is the most commonly used method in practical applications, it is an important means of reducing the moisture content of fruit, vegetables, and enhance resistance to degradation due to decrease in water activity.

And the widely use type, is convective drying ,in which the gas temperature decrease as result of uptake of water from the material[5].Hot gases are made to flow over the surface of material to be dried and thermal energy is transferred from the hot gas to the material, by convection. The vapor as result of evaporation of the moisture is transported away by the air flow.

In solids drying, there has often been a great difference between academic theory and industrial design practice. Traditionally, practical dryer design has tended to be base on simple correlations and scale-up from pilot-plant tests, rather than rigorous theoretical model [6].

However, drying processes involving solids are much more difficult to model than fluid-phase (liquid and gas) processes. Physical properties of

9

fluids can be obtained easily from databanks, and are uniquely defined for given temperatures and pressures, while the system is controlled by equilibrium thermodynamics. In contrast, physical properties of solids vary considerably with solids structure. For example, drying kinetics can differ by orders of magnitude for the same chemical substance depending on particle size, porosity, polymorph, etc. Drying is also inherently a non- equilibrium process; the moisture within a particle does not flash off instantly when exposed to hot dry air [7].

Drying of solid is generally understood to follow two distinct drying zones, namely the constant rate and the falling rate periods, with critical moisture content demarcating these zones. The constant rate period is understood to have maximum drying rate, which remain constant until the critical moisture content is reach. The extent of the drying zone is based on the type of the material. [8]

However, the drying time of a particle depends on its initial moisture content, mass flow rate, temperature at inlet, bulk density, and final moisture content [9].

3.2 Available convective types

Fluidized bed can by classify as gas-solid transport system which are characterized by a continuous convective drying process very high particle resident time.

The spouted bed type system for feeding solid into pneumatic transport tube provide a stable flow of particles, and offer operations of wide range of solid flow rate by changing air flow rate or the distance between the air inlet and the outlet in the tube transport ,unlike the normal feeder i.e.

fluidized bed type[10]. Which is very important factor for product drying .Other types are fixed bed, pneumatic conveying, and impingement superheated steam dryer e.t.c

3.3 Existing work

Ferria et al [10] analyzed one dimensional fluid dynamic model based on continuity and momentum balance for dilute gas solid flow in a pneumatic dryer with a spouted bed type solid feeding system. And concluded that one –dimensional model does not provide a good description of the physical

10

phenomena and flow pattern involved in pneumatic vertical flow transport, due to inconsistence estimated value of solid-fluid interaction forces and slip velocities .

More so, Somkiat et al [11] compared the performances of pulsed and conventional fluidized-bed dryers, based on consumption of energy, outlet moisture content, and quality of the dried. Experimental results have shown that the variation of moisture content at the exits of both dryer types in test runs was very small. Calculated thermal and electrical energy consumptions indicated that the pulsed flow dryer was more economical than the conventional dryer.

Investigation of flow field and the heat transfer characteristics in a tangential inlet cyclone mainly used for the separation of the dens phase of a two phase flow Was carried out by Irfan and Fuat [12]. It was observed that heat transfer increases at all surfaces with the inlet velocity and decreases towards to cone apex. The maximum value of the Nusselt number occurs on the region opposite to the inlet section and this region displaces down slightly for high inlet velocities.

Results obtained from their computer modeling have demonstrated that CFD is suitable for simulating the flow and heat transfer characteristics in cyclone separators.

11

4 Theoretical Analysis

The solution of heat transfer, fluid flow, and other related processes can begin when the laws governing these processes have been expressed in mathematical form, generally in terms of differential equations. The governing equations of fluid flow represent mathematical statement of the conservation laws of physics. The statements are the following:

• The mass of a fluid is conserved.

• The rate of change of momentum equals the sum of the forces on a fluid particle, which is the Newton’s second law.

• The rate of change of energy is equal to the sum of the rate of heat addition to and the rate of work done on a fluid particle, which is the first law of thermodynamics.

4.1 Derivations and Descriptions Of Equations

4.1.1 Mass conservation equations

According to Versteeg H.K and Malalasekera. [17], any fluid flow problem must fulfill the mass conservation law. It is stated that the rate of the increase of mass in fluid element is equal to the net rate of flow of mass into fluid element. This yields the mass conservation equations:

( ) ( ) ( )

u v w 0

t x y z

ρ ρ ρ ρ

∂ ∂ ∂ ∂

+ + + =

∂ ∂ ∂ ∂ (4.1)

Or in more compact vector notation

div u( ) 0 t

ρ ρ

∂ + =

∂ (4.2)

The equation (4.2) is unsteady, three-dimensional mass conservation or continuity equation at a point in a compressible fluid. The first term on the left hand side is the rate of change in time of the density (mass per unit volume). The second term describes the net flow of mass out of the element across its boundaries and is called the convective term.

4.1.2 Momentum conservation equations

Newton’s second law states that the rate of change of momentum of a fluid

12

particle equals the sum of the forces on the particle [17]. According to this law, one can derive the momentum conservation equations in x, y and z three directions：

( )

( )

( )

( )

( )

( )

xx yx zx

x

xy yy zy

y

xz yz zz

z

u p

div uu F

t x x y z

v p

div vu F

t y x y z

w p

div wu F

t z x y z

τ τ τ

ρ ρ

τ τ τ

ρ ρ

τ τ τ

ρ ρ

∂ ∂ ∂

∂ ∂

+ = − + + + +

∂ ∂ ∂ ∂ ∂

∂ ∂ ∂

∂ + = −∂ + + + +

∂ ∂ ∂ ∂ ∂

∂ ∂ ∂

∂ + = −∂ + + + +

∂ ∂ ∂ ∂ ∂

(4.3)

Where p is the pressure on the fluid element; τxx, τxy and τxz etc. are the viscosity shear components on the fluid element due to the molecule viscosity. Fx, Fy and Fz are the body force of element.

For Newtonian fluid, the viscosity stresses τ is proportional to the rate of deformation [17]. Thus:

2 ( )

2 ( )

2 ( )

xx

yy

zz

xy yy

xz zy

yz zy

u div u x

v div u y

u div u z

u v

y x

u w

z x

v w

z y

τ μ λ

τ μ λ

τ μ λ

τ τ μ

τ τ μ

τ τ μ

= ∂ +

∂

= ∂ +

∂

= ∂ +

∂

⎛∂ ∂ ⎞

= = ⎜⎝∂ +∂ ⎟⎠

∂ ∂

⎛ ⎞

= = ⎜⎝∂ + ∂ ⎟⎠

⎛∂ ∂ ⎞

= = ⎜⎝∂ + ∂ ⎟⎠

(4.4)

Where μ is dynamic viscosity, λ is the second viscosity. Substitute equation (4.4) into (4.3):

13

( )

( ) ( )

( )

( ) ( )

( )

( ) ( )

u

v

w

u p

div uu div grad u S

t x

v p

div vu div grad v S

t y

w p

div wu div grad w S

t z

ρ ρ μ

ρ ρ μ

ρ ρ μ

∂ + = −∂ +

∂ ∂

∂ ∂

+ = − +

∂ ∂

∂ ∂

+ = − +

∂ ∂

(4.5)

Where Su, Sv and Sw are the generalized source term。Su=Fx,+Sx, Su=Fy,+Sy, Su=Fz,+Sz, and the expressions of Sx, Sy and Sz is the following:

( )

( )

( )

x

y

z

u v w

S div u

x x y x z x x

u v w

S div u

x y y y z y y

u v w

S div u

x z y z z z z

μ μ μ λ

μ μ μ λ

μ μ μ λ

∂ ⎛ ∂ ⎞ ∂ ⎛ ∂ ⎞ ∂ ⎛ ∂ ⎞ ∂

=∂ ⎜⎝ ∂ ⎟⎠+∂ ⎜⎝ ∂ ⎟⎠+∂ ⎜⎝ ∂ ⎟⎠+∂

⎛ ⎞ ⎛ ⎞ ⎛ ⎞

∂ ∂ ∂ ∂ ∂ ∂ ∂

=∂ ⎜⎝ ∂ ⎟⎠+∂ ⎜⎝ ∂ ⎟⎠+∂ ⎜⎝ ∂ ⎟⎠+∂

∂ ⎛ ∂ ⎞ ∂ ⎛ ∂ ⎞ ∂ ⎛ ∂ ⎞ ∂

=∂ ⎜⎝ ∂ ⎟⎠+∂ ⎜⎝ ∂ ⎟⎠+∂ ⎜⎝ ∂ ⎟⎠+∂

(4.6)

This momentum equation is also called Navier-Stokes equations [17].

4.1.3 Energy conservation equations

The energy equation is derived from the first law of the thermodynamics which states that the rate of the change of energy of a fluid particle is equal to the rate of heat addition to the fluid particle plus the rate of the work done in the particle. One can obtain the energy conservation equation with respect of temperature T

( )

( ) _T

p

T k

div uT div grad T S

t c

ρ ρ ^{⎛ ⎞}

∂ ∂ + = ⎜⎜⎝ ⎟⎟⎠+ (4.7)

The full form is

14

( ) ( ) ( ) ( )

T

p p p

T div uT div vT div wT

t x y z

k T k T k T

x c x y c y z c z S

ρ ρ ρ ρ

∂ + + +

∂ ∂ ∂ ∂

⎛ ⎞ ⎛ ⎞ ⎛ ⎞

∂ ∂ ∂ ∂ ∂ ∂

=∂ ⎜⎜⎝ ∂ ⎟⎟⎠+∂ ⎜⎜⎝ ∂ ⎟⎟⎠+∂ ⎜⎜⎝ ∂ ⎟⎟⎠+

(4.8)

Where cp is the specific heat, k is the heat conduction coefficient. ST is called viscous dissipation term. Among the above equations, there are six variables: u, v, w, p, T and ρ. The equations of state relate the variables to the two state variables T and ρ:

p=p (T, ρ) (4.9) For a perfect gas the following, well-known, equation of state is useful:

p=ρRT (4.10) Where R is gas constant.

Now the equation group (4.1) , (4.6), (4.8), (4.9) are closed and can be solved.

4.1.4 Conservative form of the governing equations of fluid flow

To summarize the findings thus far we quote in Table 4.1 the conservative or divergence form of the system of the equations which governs the time- dependent three dimensional fluid flow and heat transfer of a compressible Newtonian fluid.

Table 4.1 governing equations of the flow of a compressible Newtonian fluid

Mass

( ) 0

div u t

ρ ρ

∂ + =

∂

x-momentum ( )

( ) ( ) _u

u p

div uu div grad u S

t x

ρ ρ μ

∂ + = −∂ +

∂ ∂

y-momentum ( )

( ) ( ) _v

v p

div vu div grad v S

t y

ρ ρ μ

∂ + = −∂ +

∂ ∂

z-momentum ( )

( ) ( ) _w

w p

div wu div grad w S

t z

ρ ρ μ

∂ + = −∂ +

∂ ∂

Energy ( )

( ) _T

p

T k

div uT div grad T S

t c

ρ ρ ^{⎛ ⎞}

∂ ∂ + = ⎜⎜⎝ ⎟⎟⎠+

Equation of state p=p(T, ρ)

15

4.1.4 The general form of the governing equations of fluid flow

For the convenience of analysis of the governing equation, and possibility of numerically solving every equation by one program code, it is of very importance to establish the general form of the governing equations.

Via the observation the equation group (4.1) , (4.6) and (4.8), it can be found that they all represent the conservation law of physics quantity in unit time and volume, though there are different variables. If φ is used to express the generalized variables, the each governing equation can be written as the following format:

( )

( ) ( )+S

div u div grad t

ρφ ρ φ φ

∂ + = Γ

∂ (4.10) The full pattern is:

( ) div u( ) ( v ) ( w )

t x y z

x x y y z z S

ρφ ρ φ ρ φ ρ φ

φ φ φ

∂ ∂ ∂

+ + +

∂ ∂ ∂ ∂

⎛ ⎞

∂ ⎛ ∂ ⎞ ∂ ∂ ∂ ⎛ ∂ ⎞

=∂ ⎜⎝Γ∂ ⎟⎠+∂ ⎜⎝Γ∂ ⎟⎠+∂ ⎜⎝Γ∂ ⎟⎠+

(4.11)

Where φ is a generalized variable could be u, v, w, T or other variables; Γ is the generalized diffusion coefficient; S is generalized source term. In equation (4.10), the terms from left to right are transient term, convective term, diffusive tern and source term respectively. For specific equation,φ, Γ and S are the specific form, the table 4.2 gives the corresponding relationship of the three symbol and the specific equations.

Table 4.2 Symbol pattern in generalized governing equation

φ _{Γ S}

continuity 1 0 0

momentum ui μ -∂ ∂ + p x/ _i S_i

energy T k/c ST

Therefore, if the code can numerically solve the equation (4.10), it can solve the different type fluid flow and heat transfer problem. For differentφ, just call the program repeatedly, and give the suitable boundary and initial conditions, the result can be obtained.

16

5. Available Experimental Result

Data are obtained from factory; Kunming shipbuilding equipment co. Such as the temperature velocity of gas particle and these value are serves as input to the analysis and the final results are also compared.

The values are as follows:

Inlet velocity 18m/s Temperature 600K

Residence time of particle 5-12 seconds Length of pipe 18 meter.

The temperature measurement is obtained as standard value, used in drying of a specified food particle, Kunming shipbuilding equipment Co. design for their client base on demand specifications, so the design specification on which this studies is base is the type which the drying chamber 18 meter long pipe . And the value of the velocity is standardized as well, for this length.

But the drying time is seriously affected by this parameters or values, since high velocity give lower resident time and lower velocity give high residence time vice versa.

The drying time is the actual time the food particle spend in the equipment, and the time of transporting the particle to and away from the equipment is not measured.

The residence time is measured with the help of stop watch which is switched on when the food particle get inside the drying chamber, where drying and heat exchange take place, and stop immediately the particle leave the chamber and this is done repeatedly but the average is used.

However, the time to transport the particle to and away from the drying chamber is not included due to complexity of the whole structure.

17

6 Numerical Analysis

6.1 Control-Volume Formulation (CVF) or Finite Volume Method (FVM)

Some elementary textbooks on heat transfer, derive the finite-difference equation via the Taylor-series method, and then demonstrate that, the resulting equations are consistent with heat balance over a small region surrounding a grid point. Also, the control-volume formulation can be regarded as a special version of the method of weighted residuals. The basic idea of the control-volume formulation is easy to understand and lends itself to direct physical interpretation. The calculation domain is divided into a number of nonoverlapping control volumes such that there is one control volume surrounding each grid point. The differential equation is integrated over each control volume. Piecewise profiles expressing the variation of φ between the grid points are used to evaluate the required integrals. The result is the discretization equation containing the values of

φ for a group of grid points.

The discretization equation obtained in this manner expresses the conservation principle for φ for the finite control volume, just as the differential equation expresses it for an infinitesimal control volume.

Indeed, deriving the control-volume discretization equation by integrating the differential equation over a finite control volume is a rather roundabout process, of DISCRETIZATION METHODS.

The most attractive feature of the control-volume formulation is that the resulting solution would imply that the integral conservation of quantities such as mass, momentum, and energy is exactly satisfied over any group of control volumes and, of course, over the whole calculation domain. This characteristic exists for any number of grid points-not just in a limiting sense when the number of grid points becomes large. Thus, even the coarse- grid solution exhibits exact integral balances.

When the discretization equations are solved, to obtain the grid-point values of the dependent variable, the result can be viewed in two different ways. In the finite-element method and in most weighted-residual methods, the assumed variation of φ consisting of the grid-point values and the interpolation functions (or profiles) between the grid points is taken as the approximate solution. In the finite-difference method, however, only the

18

grid-point values of φ are considered to constitute the solution, without any explicit reference as to how φ varies between the grid points. This is akin to a laboratory experiment, where the distribution of a quantity is obtained in terms of the measured values at some discrete locations, without any statement about the variation between these locations. This view is adopted in the present control-volume approach. The solution is sought in the form of the grid-point values only. The interpolation formulas or the profile assumptions can be forgotten. This viewpoint permits complete freedom of choice in employing, base on wish, different profile assumptions for integrating different terms in the differential equation.

To make the foregoing discussion more concrete, below are derivations of the control-volume discretization equation for a simple situation.

6.2 The SIMPLE algorithm

The acronym SIMPLE stands for Semi-Implicit Method for Pressure- Linked Equations. The algorithm was originally put forward by Patanker and Spalding (1972) and is essentially a guess-and correct procedure for the calculation of pressure on the staggered grid arrangement introduced above.

The method is illustrated by considering the two-dimensional laminar steady flow equations in Cartesian Co-ordinates.

To initiate the SIMPLE calculation process a pressure field p^* is guessed.

Discretised momentum equations (6.8) and (6.10) are solved using the guessed pressure field to yield velocity components u^* and v^* as follows:

* * * *

, , ( 1, , ) , ,

i J i J nb nb I J I J i J i J

a u =

∑

a u + p₋ − p A +b (6.12) a v_{I j I j,} ^*_, =

∑

a v_{nb nb}^* +(p^*_{I J, −1}− p^*_{I J,} )A_{I j,} +b_{I j,} (6.13)

Now we define the correction p′ as the difference between the correct pressure field p and the guessed pressure field p^*, so that

p= p*+ p′ (6.14)

Similarly we define velocity corrections u^´ and v^´ to relate the correct velocities u and v to the guessed velocities u ^*and v^*

u u= *+u′ (6.15)

v v= +* v′ (6.16)

Substitution of the correct pressure field p into the momentum equations yields the correct velocity field (u,v). Discretized equations (6.8) and (6.10) link the correct velocity fields with correct pressure field.

19

Subtraction of equations (6.12) and (6.13) from (6.8) and (6.10), respectively, gives

* * * *

, , , 1, 1, , , ,

* * * *

, , , , 1 , 1 , , ,

, , 1, ,

, , , 1 ,

, ,

( ) ( ) [( ) ( )]

( ) ( ) [( ) ( )]

( )

( )

i J i J i J nb nb nb I J I J I J I J i J

I j I j I j nb nb nb I J I J I J I J I j

nb nb

i J i J I J I J

I j I j I J I J

i i J

a u u a u u p p p p A

a v v a v v p p p p A

a v

u d p p

v d p p

d A

− −

− −

−

−

− = − + − − −

− = − + − − −

′

′ = ′ − ′

′ = ′ − ′

=

∑ ∑

∑

,

, ,

,

*

, , , 1, ,

*

, , , , 1 ,

( )

( )

J i J I j I j

I j

i J i J i J I J I J

I j I j I j I J I J

a d A

a

u u d p p

v v d p p

−

−

=

′ ′

= + −

′ ′

= + −

(6.17)

a_{I j,} (v_{I j,} −v^*_{I j,} )=

∑

a v_nb( _nb−v^*_nb) [(+ p_{I J, −1}−p^*_{I J, −1}) (− p_{I J,} −p^*_{I J,} )]A_{I j,} ^(6.18) Using correction formulae (6.14-6.16) the equations (6.17-6.18) may be rewritten as follows:

, , ( 1, , ) ,

i J i J nb nb I J I J i J

a u′ =

∑

a u′ + p′₋ −p′ A (6.19)

, ( , 1 , ) ,

I j nb nb I J I J I j

v′ =

∑

a v′ + p′ ₋ − p′ A (6.20)

At this point an approximation is introduced:

∑

a u′_{nb nb}^and

∑

^{a v′nb nb are}

dropped to simplify equations (6.19) and (6.20) for the velocity corrections.

Omission of these terms is the main approximation of the SIMPLE algorithm. We obtain

, , ( 1, , )

i J i J I J I J

u′ =d p′₋ −p′ (6.21)

, , ( , 1 , )

I j I j I J I J

v′ =d p′ ₋ −p′ (6.22)

Where _, ^,

, i J i J

i J

d A

= a and _, ^,

, I j I j

I j

d A

= a (6.23.)

Equation (6.21) and (6.22) describe the corrections to be applied to velocities through formulae (6.15) and (6.16), which gives

20

*

, , , ( 1, , )

i J i J i J I J I J

u =u +d p′₋ − p′ (6.24)

*

, , , ( , 1 , )

I j I j I j I J I J

v =v +d p′ ₋ −p′ (6.25) Similar expressions exist for ui+1,J and υI,j+1:

) ( _{, 1,}

, , 1

* 1 ,

1J i J i J I J I J

i u d p P

u+ = + + ₊ ′ − ′₊ (6.26) )

( _{, , 1}

1 , 1 , 1

, ∗ + +

+

+ = _{I j} + _{I j I}′_J − _I′_J

J

I ν d P P

ν (6.27)

Where

J i

J i J

i a

d A

, 1

, 1 ,

1

+ +

+ = and

1 ,

1 , 1 ,

+ +

+ =

j I

j I j

I a

d A (6.28)

Thus far we have only considered the momentum equations but, as mentioned earlier, the velocity field is also subject to the constraint that it should satisfy continuity equation (6.3). Continuity is satisfied in discretised form for the scalar control volume show in Figure 6.1:

[(ρuA)_i+1,J −(ρuA)_i,J]+[(ρνA)_I,j+1 −(ρνA)_I,j]=0 (6.29)

Figure 6.1 The scalar control volume used for the discretisation of the continuity equation.

21

Substitution of the corrected velocities of equations (6.24-6.27) into discretised continuity equation (6.29) gives

1, 1, 1, 1 , 1,

J , , , I 1 J ,

, 1 I , 1 , 1 , , 1

, , , , , 1 ,

[ A ( ( ))

A ( ( ))]

[ A ( ( ))

( ( ))] 0

i J i J i J i I J I J

i i J i J i J I J

I j I j I j I J I J

I j i j i j I j I J I J

u d p p

u d p p

d p p

A d p p

ρ ρ

ρ ν

ρ ν

∗

+ + + + +

∗

−

+ ∗ + + +

∗ ∗

−

′ ′

+ −

′ ′

− + −

′ ′

+ + −

′ ′

− + − =

，J

， ，

，j +1

(6.30)

This may be re-arranged to give

_{, 1}

1, , , 1 ,

1, 1, , 1, , 1

, 1, , , 1

[( ) ( ) ( ) ( ) ]

(( ) ( ) ( )

[( ) ( ) ( ) ( ) ]

I j

i J i J I j I j

i J I J i J i J I J

i J i J I j I j

dA dA dA dA

dA p dA p dA p

u A u A A A

ρ ρ ρ ρ

ρ ρ ρ

ρ ρ ρν ρν

+

+ +

+ + − +

∗ ∗ ∗ ∗

+ +

+ + +

′ ′ ′

= + +

+ − + −

(6.31)

Identifying the coefficients of p′ this may be written as

^{, , 1, 1, 1, 1, , 1 , 1}

, 1 , 1 ,

I J I J I J I J I J I J I J I J

I J I J I J

a p a p a p a p

a p b

+ + − − + +

− −

′ = ′ + ′ + ′

′ ′

+ + (6.32)

Where a_{I J,} =a_I+1,J +a_I−1,J +a_{I J, +1}+a_{I J, −1} and the coefficients are given below:

1,

I J

a₊ a_I−1,J a_{I J, +1} a_{i J, −1} (ρdA)_i+1,_J (ρdA)_{i J,} (ρdA)_{I j, +1} (ρdA)_{I , j}

,

b′I J

, 1, , , 1

(ρu A^∗ )_{i J} −(ρu A^∗ )_{i+ J} +(ρν^∗A)_{I j}−(ρν^∗A)_{I j+}

Equation (6.32) represents the discretised continuity equation as an equation for pressure correction p′ . The source term b′ in the equation is the continuity imbalance arising from the incorrect velocity field u^∗,ν^∗, By solving equation (6.32), the pressure correction field p′ can be obtained at all points. Once the pressure correction field is known, the correct pressure field may be obtained using formula (6.14) and velocity components through correction formulae (6.24-6.27). The omission of terms such as

nb nb

∑

a u′ in the derivation does not affect the final solution because the pressure correction and velocity corrections will all be zero in a converged solution giving p^∗ = ,up ^∗ = andu ν^∗= . ν

The pressure correction equation is susceptible to divergence unless some under-relaxation is used during the iterative process and new, improved,

22

pressures p^new are obtained with

new

p = p^∗+a pp ′ (6.33) Where α_p is the pressure under-relaxation factor. If we select α_p equal to I the guessed pressure field p^∗ is corrected byp′ . However, the correctionsp′ , in particular when the guessed field p^∗ is far away from the final solution, is often too large for stable computations. A value of α_p equal to zero would apply no correction at all, which is also undesirable.

Taking α_p between 0 and 1 allows us to add to guessed field p^∗ a fraction of the correction field of the correction field p′ that is large enough to move the iterative improvement process forward, but small enough to ensure stable computation.

The velocities are also under-relaxed. The iteratively improved velocity components u^new and ν^new are obtained from

u^new=α_uu+ −(1 α_u)u⁽ⁿ⁻¹⁾ (6.34)

( 1)

(1 )

new n

ν ν

ν =α ν+ −α ν ⁻ (6.35) Where α_u and α_ν are the u- and ν -velocity under-relaxation factors with values between 0 and 1, u and ν are the corrected velocity components without relaxation and u⁽ⁿ⁻¹⁾ and ν⁽ⁿ⁻¹⁾ represent their values obtained in the previous iteration. After some algebra it can be shown that with under- relaxation the discretised u-momentum equation takes the form

, , ( 1)

, ( 1, , ) , (1 ) ,

i J i J n

i J nb nb I J I J i J i u i J

u u

a a

u a u p p A b u

a α

α

−

−

⎡ ⎤

= + − + +⎢ − ⎥

⎣ ⎦

∑

^，J (6.36)

and the discretised ν -momentum equation

a^{I j,} _{I j, nb nb} ( _{I J, 1 I J,} ) _{I j, I j,} (1 )a^{I j,} _{I j}⁽_,^{n 1)}

a p p A b

a_ν ν ν α^ν _ν ν

α ⁻

−

⎡ ⎤

= + − + +⎢ − ⎥

⎣ ⎦

∑

(6.37)

The pressure correction equation is also affected by velocity under- relaxation and it can be shown that d-terms of pressure correction equation become

, ,

, i J u i J

i J

d A

a

= α , _1, ^1,

1,J

i J u

i J

i

d A

a

+ α

+

+

= , _, ^,

, I j I j

I j

d A

a αν

= and

23

, 1 , 1

, 1 I j I j

I j

d A

a αν + +

+

=

Note that in these formulae a_{i J,} ,a_i+1,J,a_{I j,} and a_{I j, +1} are the central coefficients of discretised velocity equations at positions (i,J), (i+1,J), (I,j) and (I,j+1) of a scalar cell centered around P.

A correct choice of under-relaxation factors α is essential for cost-effective simulations. Too large a value of α may lead to oscillatory or even divergent iterative solutions and a value which is too small will cause extremely slow convergence. Unfortunately, the optimum values of under- relaxation factors are flow dependent and must be sought on a case-by-case basis.

The SIMPLE algorithm gives a method of calculating pressure and velocities, The method is iterative and when other scalars are coupled to the momentum equations, the calculation needs to be done sequentially. The sequence of operations in a CFD procedure which employs the SIMPLE algorithm is given in Figure 6.2

24

Initial guessp u^∗, , ,^∗ ν^∗ ∅^∗

, u^∗ν^∗

_u^∗_,ν^∗

p′

p u, , ,ν ∅^∗

No ∅

yes

Figure 6.2 the SIMPLE Algorithm Flow chart START

SETP 3: Correct pressure and velocities

, , ,

, , , 1, ,

, , , 1, ,

( )

( )

I J I J I J

I J i J i J I J I J

I j I j i J I J I J

p p p

u u d p p

d p p

ν ν

∗

∗

−

∗

−

= + ′

′ ′

= + −

′ ′

= + −

STEP 1：Solve discretised momentum equations

, , ( 1, , ) , ,

i J i J nb nb I J I J i J i J

a u^∗ =

∑

a u^∗ + p^∗₋ − p^∗ A +b

, , ( , 1 , ) , ,

I J I J nb nb I J I J I J I J

a ν^∗ =

∑

a ν^∗ + p^∗ ₋ − p^∗ A +b

STEP 2: Solve pressure correction equation

, , 1, 1, 1, 1, , 1 , 1 , 1 , 1 ,

I J I J I J I J I J I J I J I J I J I J I J

a p′ =a₋ p′₋ +a₊ p′₊ +a ₋ p′ ₋ +a ₊ p′ ₊ +b′

STEP 4：Solve all other discretised transport equations

, , 1, 1, 1, 1, , 1 , 1 , 1 , 1 ,

I J I J I J I J I J I J I J I J I J I J I J

a ∅ =a ₋ ∅ ₋ +a₊ ∅ ₊ +a ₋∅ ₋ +a ₊∅ ₊ + ∅b

Convergence?

STOP ,

, Set

P p u u

ν ν

∗ ∗

∗ ∗

= =

= ∅ = ∅

25

6.3 Implementation of boundary conditions

6.3.1 Grid arrangement

• The physical boundary coincide with the scalar control volume

• The grid extends outside the physical boundary.

6.3.2 Inlet boundary conditions

To accommodate the inlet boundary condition for these variables, it is unnecessary to make any modifications to their discretized equations. The link with the boundary side is cut in the discretized pressure correction equation by setting the boundary side coefficient aw (if the inlet is on the west) equal to zero. Since the velocity is known at inlet, it is also not necessary to make a velocity correction here and hence we have

u^*_w=u_w

6.3.3 Outlet boundary conditions

The outlet is such a region that the flow often teaches a fully developed state where no change occurs in the flow direction. It means all the variables (except pressure), the gradient is 0.

For the v- and scalar equations this implies setting

v_{NI j,} =v_NI−1,j ; φ_{NI J,} =φ_NI−1,J

Where NI is the last line of the nodes.

Special care should be taken in the case of the u-velocity. Calculation of u at the outlet plane i=NI by assuming a zero gradient gives

u_{NI j,} =u_NI−1,j

During the iteration cycles of the SIMPLE algorithm there is no guarantee that these velocities will conserve mass over the computational domain as a whole. To ensure that overall continuity is satisfied the total mass flux going out of the domain (Mout ) is firstly computed by summing all the extrapolated outlet velocities. To make the mass flux out equal to the mass flux Min coming into the domain all the outlet velocity components uNI, J are multiplied by the ratio Min/ Mout. Thus the outlet plane velocities with the continuity correction are given by

26 u_{NI j,} =u_NI−1,j×M_in/ M_out

The velocity at the outlet boundaries is not corrected by means of pressure corrections. Hence in the discretised p’-equation the link to the outlet (if it is on the east) boundary side is suppressed by setting aE=0. The contribution to the source term in this equation is calculated as normal, noting that

u^*_E =u_E

No addition modifications are required.

6.4.4 Outlet boundary conditions

The wall is the most common boundary encountered in confined fluid flow problems. The no-slip condition is the appropriate condition for the velocity components at solid walls. Which is

u=v=0

Since the wall velocity is known, it is also unnecessary to perform a pressure correction here. In the discretised p^’-equation for the cell nearest to the wall, the wall link (if it is on the south) is cut by setting aS = 0, and in the source term

v_S^* =v_S

6.4 Drying in Straight Pipe Model

The straight type is the common type which is too high, often difficult to fabricate in the workshop is shown in the diagram below

27

Figure 6.3 straight type dryer Kunming shipbuilding equipment Co

In the above equipment，only the cavity where drying occur; that is, the part labeled straight pipe dryer is modeled and solved in Fluent and hand written code is used to validate the solution. The results of the single and double phase are presented in chapter 8.

However, other components shown in the above picture are mechanisms for