Working on my undergraduate degree in Mechanical Engineering, my favorite class was heat transfer. For a class project, I built a simple simulation that showed convection for a randomly generated temperature field. At the time, I didn’t have enough programming skills to build a quality simulation. The purpose of this is to revisit that project with my improved science and programming skills. Along the way, I will teach the math behind building the model. The first step will be to research the fundamental equations for this system.

Problem Definition

In this project, we will model a compressible fluid. Natural convection forces will affect the fluid. The simulation will model the transient behavior of a starting temperature or input heat. The simulation will be 2-Dimensional to display the results easier. Lastly, the model will use constant Cartesian cell structure to simplify the problem.

Governing Equations

The Navier-Stokes momentum equations for 2D Cartesian flow will be used to determine the fluid velocities

$\rho \left( \frac { du }{ dt } +u\frac { du }{ dx } +v\frac { du }{ dy } \right) =-\frac { dp }{ dx } +\mu \left( \frac { { d }^{ 2 }u }{ { dx }^{ 2 } } +\frac { { d }^{ 2 }u }{ { dy }^{ 2 } } \right) +\rho { g }_{ x } \\ \rho \left( \frac { dv }{ dt } +u\frac { dv }{ dx } +v\frac { dv }{ dy } \right) =-\frac { dp }{ dy } +\mu \left( \frac { { d }^{ 2 }v }{ { dx }^{ 2 } } +\frac { { d }^{ 2 }v }{ { dy }^{ 2 } } \right) +\rho { g }_{ y }$

Conservation of mass will be used to determine fluid densities

$\frac { d\rho }{ dt } +\rho \left( \frac { du }{ dx } +\frac { dv }{ dy } \right) =0$

A conservation of energy equation was developed by including conduction, convection and external heat to the first law of thermodynamics

$\rho { C }_{ v }\left( \frac { dT }{ dt } +u\frac { dT }{ dx } +v\frac { dT }{ dy } \right) =k\left( \frac { { d }^{ 2 }T }{ { dx }^{ 2 } } +\frac { { d }^{ 2 }T }{ { dy }^{ 2 } } \right) +S$

The ideal gas law will be used to determine pressure

$p=\rho { R }_{ specific }T$

Finite Difference Approximation

A finite difference approximation will be used to solve the fluid state at future time intervals. A small time step and cell size will be desired to maximize accuracy. Central difference approximations will be used to reduce any directional biases.

When using a finite difference, we can encounter a problem. The grid below uses the same cells to represent both pressure and velocity.

Since we are using a central difference approximation, the velocities of center square is based on the pressure difference between it’s neighbors. We can see how the equation does not compare the squares current pressure between it’s neighbor squares.

${ u }_{ i,j }^{ t+1 }={ u }_{ i,j }^{ t }+\Delta t\left( -\frac { { p }_{ i+1,j }^{ t }-{ p }_{ i-1,j }^{ t } }{ \rho 2\Delta x } +\frac { \mu }{ \rho } \left( \frac { { u }_{ i+1,j }^{ t }+{ u }_{ i-1,j }^{ t }-2u_{ i,j }^{ t } }{ { \Delta x }^{ 2 } } +\frac { { u }_{ i,j+1 }^{ t }+{ u }_{ i,j-1 }^{ t }-2u_{ i,j }^{ t } }{ { \Delta y }^{ 2 } } \right) +{ g }_{ x }-{ u }_{ i,j }^{ t }\frac { { u }_{ i+1,j }^{ t }-{ u }_{ i-1,j }^{ t } }{ 2\Delta x } -{ v }_{ i,j }^{ t }\frac { { u }_{ i,j+1 }^{ t }-{ u }_{ i,j-1 }^{ t } }{ 2\Delta y } \right)$

If there is not a pressure differential between the 2 opposite neighbor cells then we do not get flow. As a result, a checkerboard pattern will emerge. The image below shows a possible solution.

This solution solves the governing equations but creates a system that is unrealistic. In a real system, we would see flow from the 4-pressure elements into the 2-pressure elements.

A better solution uses different, overlapping areas for the pressures and velocities. In the improved solution, the pressures exist in each area of the grid and the velocities represent the flows between the pressure cells. This solution is represented below.

Solved Equations

With our improved grid model, we can finally solve the governing equations using a finite difference approximation. The following equations can calculate the next simulation conditions based on the current conditions.

Navier-Stokes equations for velocity

{ u }_{ i,j }^{ t+1 }={ u }_{ i,j }^{ t }+\Delta t\left( -\frac { { p }_{ i,j }^{ t }-{ p }_{ i-1,j }^{ t } }{ \rho \Delta x } +\frac { \mu }{ \rho } \left( \frac { { u }_{ i+1,j }^{ t }+{ u }_{ i-1,j }^{ t }-2u_{ i,j }^{ t } }{ { \Delta x }^{ 2 } } +\frac { { u }_{ i,j+1 }^{ t }+{ u }_{ i,j-1 }^{ t }-2u_{ i,j }^{ t } }{ { \Delta y }^{ 2 } } \right) +{ g }_{ x }-{ u }_{ i,j }^{ t }\frac { { u }_{ i+1,j }^{ t }-{ u }_{ i-1,j }^{ t } }{ 2\Delta x } -\left( { v }_{ i,j }^{ t }+{ v }_{ i-1,j }^{ t } \right) \frac { \left( { u }_{ i,j }^{ t }-{ u }_{ i,j-1 }^{ t } \right) }{ 4\Delta y } -\left( { v }_{ i,j+1 }^{ t }+{ v }_{ i-1,j+1 }^{ t } \right) \frac { \left( { u }_{ i,j+1 }^{ t }-{ u }_{ i,j }^{ t } \right) }{ 4\Delta y } \right) \\ { v }_{ i,j }^{ t+1 }=v_{ i,j }^{ t }+\Delta t\left( -\frac { { p }_{ i,j }^{ t }-{ p }_{ i,j-1 }^{ t } }{ \rho \Delta y } +\frac { \mu }{ \rho } \left( \frac { { v }_{ i+1,j }^{ t }+{ v }_{ i-1,j }^{ t }-2v_{ i,j }^{ t } }{ { \Delta x }^{ 2 } } +\frac { { v }_{ i,j+1 }^{ t }+{ v }_{ i,j-1 }^{ t }-2v_{ i,j }^{ t } }{ { \Delta y }^{ 2 } } \right) +{ g }_{ y }-{ v }_{ i,j }^{ t }\frac { { v }_{ i,j+1 }^{ t }-{ v }_{ i,j-1 }^{ t } }{ 2\Delta y } -\left( { u }_{ i,j }^{ t }+{ u }_{ i,j-1 }^{ t } \right) \frac { \left( { v }_{ i,j }^{ t }-{ v }_{ i-1,j }^{ t } \right) }{ 4\Delta x } -\left( { u }_{ i+1,j }^{ t }+{ u }_{ i+1,j-1 }^{ t } \right) \frac { \left( { v }_{ i+1,j }^{ t }-{ v }_{ i,j }^{ t } \right) }{ 4\Delta x } \right)

Conservation of mass for density

$\rho _{ i,j }^{ t+1 }=\rho _{ i,j }^{ t }+\frac { \Delta t }{ 2 } \left( \frac { u_{ i,j }^{ t } }{ \Delta x } \left( \rho _{ i,j }^{ t }+\rho _{ i-1,j }^{ t } \right) +\frac { u_{ i+1,j }^{ t } }{ \Delta x } \left( \rho _{ i,j }^{ t }+\rho _{ i,j-1 }^{ t } \right) +\frac { v_{ i,j }^{ t } }{ \Delta y } \left( \rho _{ i,j }^{ t }+\rho _{ i+1,j }^{ t } \right) -\frac { v_{ i,j+1 }^{ t } }{ \Delta y } \left( \rho _{ i,j }^{ t }+\rho _{ i,j+1 }^{ t } \right) \right)$

Conservation of energy for temperature

$T_{ i,j }^{ t+1 }=\frac { \Delta t }{ \rho _{ i,j }^{ t } } \left( \frac { u_{ i,j }^{ t } }{ 2\Delta x } \left( T_{ i-1,j }^{ t }-T_{ i,j }^{ t } \right) \left( \rho _{ i-1,j }^{ t }+\rho _{ i,j }^{ t } \right) +\frac { u_{ i+1,j }^{ t } }{ 2\Delta x } \left( T_{ i+1,j }^{ t }-T_{ i,j }^{ t } \right) \left( \rho _{ i+1,j }^{ t }+\rho _{ i,j }^{ t } \right) +\frac { v_{ i,j }^{ t } }{ 2\Delta y } \left( T_{ i,j-1 }^{ t }-T_{ i,j }^{ t } \right) \left( \rho _{ i,j-1 }^{ t }+\rho _{ i,j }^{ t } \right) +\frac { v_{ i,j+1 }^{ t } }{ 2\Delta y } \left( T_{ i,j+1 }^{ t }-T_{ i,j }^{ t } \right) \left( \rho _{ i,j+1 }^{ t }+\rho _{ i,j }^{ t } \right) +\frac { k }{ { C }_{ v } } \left( \frac { T_{ i-1,j }^{ t }+T_{ i+1,j }^{ t }-2T_{ i,j }^{ t } }{ { \Delta x }^{ 2 } } +\frac { T_{ i,j-1 }^{ t }+T_{ i,j+1 }^{ t }-2T_{ i,j }^{ t } }{ { \Delta y }^{ 2 } } \right) +{ S }_{ i,j } \right)$

Ideal gas law for pressure

$p_{ i,j }^{ t+1 }=\rho _{ i,j }^{ t+1 }{ R }_{ specific }T_{ i,j }^{ t+1 }$

These equations can directly be used in our simulation. Beyond these equations, we will have to consider the boundary conditions. The boundaries can be modeled in many different ways so we will review them later. With all that in mind, we can begin the design and development of our simulation.

Problem Definition

Governing Equations

Finite Difference Approximation

Solved Equations

Symbol Definitions

Leave a Reply Cancel reply