MATLAB code that solves for the velocity, pressure, and temperature for laminar flow between parallel plates using the Finite Volume Method (FVM) and the SIMPLE algorithm. Details of the equations implemented are included on the attached pdf file. The inlet condition is uniform flow and temperature. The boundary condition for the top wall can be either fixed temperature or fixed heat flux. The boundary condition for the bottom wall is symmetry (for half way between plates). All gradients except pressure are zero at the outlet.
Example On Using This Code
% Geometry H = 0.01; % Height of channel in y-direction(m) L = 10*H; % Length of cavity in x-direction (m) % Grid geometry Nx = 200; % Number of main grid points in x-direction within domain dx = L/Nx; % Grid spacing in x-direction Ny = 40; % Number of main grid points in y-direction within domain dy = H/Ny; % Grid spacing in y-direction dz = 0.01; % Width in z-direction for flux calculations (m)
% Properties (air at STP) rho = 1.2; % Density (kg/m^3) mu = 1.8e-5; % Absolute viscosity (N-s/m^2) nu = mu/rho; % Kinematic viscosity (m^2/s) kt = 0.025; % Thermal conductivity (W/m-K) cp = 1006; % Specific heat (J/kg-K) alpha = kt/(rho*cp); % Thermal diffusivity (m^2/s) Pr = nu/alpha; % Prandtl number % Boundary conditions Re = 100; % Reynolds number U = Re*nu/(2*H); % Average velocity(m/s) Ti = 20; % Inlet temperature (deg. C) Tw = 100; % Wall temperature (deg. C) qw = 100; % Wall heat flux (W/m^2) BC_N = 1; % BC_N = 0 for Tw, BC_N = 1 for qw BC_S = 1; % BC_S = 0 for Tw, BC_S = 1 for symmetry % Solution controls alphaU = 0.3; % Velocity relaxation (under) alphaP = 0.2; % Pressure relaxation (under) NmaxSIM = 1e+4; % Iteration max for SIMPLE algorithm (-) NmaxGSI = 1e+1; % Iteration max for numerical method (-) err = 1e-5; % Convergence criteria (-) div = 1e+1; % Divergence criteria (-)
Value of x-Velocity, y-Velocity and Pressure at each Iteration: Contours of Temperature, Axial Temperature Variation and Dimensionless Temperature Profiles: Contour Lines of X-Velocity, Y-Velocity, Pressure and Dimensionless Velocity Profiles: Residuals of X-Velocity, Y-Velocity and Pressure:
In computational fluid dynamics (CFD), the SIMPLE algorithm is a widely used numerical procedure to solve the Navier-Stokes equations. SIMPLE is an acronym for Semi-Implicit Method for Pressure Linked Equations. The SIMPLE algorithm was developed by Prof. Brian Spalding and his student Suhas Patankar at Imperial College, London in the early 1970s. Since then it has been extensively used by many researchers to solve different kinds of fluid flow and heat transfer problems.
Many popular books on computational fluid dynamics discuss the SIMPLE algorithm in detail. A modified variant is the SIMPLER algorithm (SIMPLE Revised), that was introduced by Patankar in 1979.
The algorithm is iterative. The basic steps in the solution update are as follows:
1. Set the boundary conditions.
2. Compute the gradients of velocity and pressure.
3. Solve the discretized momentum equation to compute the intermediate velocity field.
4. Compute the uncorrected mass fluxes at faces.
5. Solve the pressure correction equation to produce cell values of the pressure correction.
6. Update the pressure field:
where urf is the under-relaxation factor for pressure.
7. Update the boundary pressure corrections p’_(b).
8. Correct the face mass fluxes:
9. Correct the cell velocities:
; where nabla_p’ is the gradient of the pressure corrections, a_(P)^(v) is the vector of central coefficients for the discretized linear system representing the velocity equation and Vol is the cell volume.
10. Update density due to pressure changes.
About Finite Volume Method (FVM)
The finite volume method (FVM) is a method for representing and evaluating partial differential equations in the form of algebraic equations [LeVeque, 2002; Toro, 1999]. Similar to the finite difference method or finite element method, values are calculated at discrete places on a meshed geometry. “Finite volume” refers to the small volume surrounding each node point on a mesh. In the finite volume method, volume integrals in a partial differential equation that contain a divergence term are converted to surface integrals, using the divergence theorem. These terms are then evaluated as fluxes at the surfaces of each finite volume. Because the flux entering a given volume is identical to that leaving the adjacent volume, these methods are conservative. Another advantage of the finite volume method is that it is easily formulated to allow for unstructured meshes. The method is used in many computational fluid dynamics packages.
Consider a simple 1D advection problem defined by the following partial differential equation
Here, rho represents the state variable and f represents the flux or flow of rho . Conventionally, positive f represents flow to the right while negative f represents flow to the left. If we assume that equation (1) represents a flowing medium of constant area, we can sub-divide the spatial domain, x, into finite volumes or cellswith cell centres indexed as i. For a particular cell, i, we can define the volume average value of rho_(i) at time t = t_(1) and
and at time t = t_(2) as,
where x_(i – 0.5) and x_(i + 0.5) represent locations of the upstream and downstream faces or edges respectively of the ith cell. Integrating equation (1) in time, we have:
To obtain the volume average of rho(x, t) at time t = t_(2), we integrate rho(x, t2) over the cell volume,
and divide the result by
We assume that f is well behaved and that we can reverse the order of integration. Also, recall that flow is normal to the unit area of the cell. Now, since in one dimension
, we can apply the divergence theorem, i.e.
, and substitute for the volume integral of the divergence with the values of f(x) evaluated at the cell surface (edges x_(i – 0.5) and x_(i + 0.5) of the finite volume as follows:
.We can therefore derive a semi-discrete numerical scheme for the above problem with cell centres indexed as i, and with cell edge fluxes indexed as i +- 0.5, by differentiating (6) with respect to time to obtain:
where values for the edge fluxes,f_(i +- 0.5), can be reconstructed by interpolation or extrapolation of the cell averages. Equation (7) is exact for the volume averages; i.e., no approximations have been made during its derivation.
This method can also be applied to a 2D situation by considering the north and south faces along with the east and west faces around a node.
General Conservation Law
We can also consider the general conservation law problem, represented by the following PDE,
Here,u represents a vector of states and f represents the corresponding flux tensor. Again we can sub-divide the spatial domain into finite volumes or cells. For a particular cell,i, we take the volume integral over the total volume of the cell, v_(i), which gives,
On integrating the first term to get the volume average and applying the divergence theorem to the second, this yields
where S_(i) represents the total surface area of the cell and n is a unit vector normal to the surface and pointing outward. So, finally, we are able to present the general result equivalent to (8), i.e.
Again, values for the edge fluxes can be reconstructed by interpolation or extrapolation of the cell averages. The actual numerical scheme will depend upon problem geometry and mesh construction. MUSCL reconstruction is often used in high resolution schemes where shocks or discontinuities are present in the solution.
Finite volume schemes are conservative as cell averages change through the edge fluxes. In other words, one cell’s loss is another cell’s gain!
- “SIMPLE solver for driven cavity flow problem” (PDF). Retrieved 2011-08-21.
- Mangani, L.; Bianchini, C. (2007). Heat transfer applications in turbomachinery (PDF). Proceedings of the OpenFOAM International Conference 2007. Retrieved 2016-03-16.
- Patankar, S. V. (1980). Numerical Heat Transfer and Fluid Flow. Taylor & Francis. ISBN 978-0-89116-522-4.
- Ferziger, J. H.; Peric, M. (2001). Computational Methods for Fluid Dynamics. Springer-Verlag. ISBN 978-3-540-42074-3.
- Tannehill, J. C.; Anderson, D. A.; Pletcher, R. H. (1997). Computational Fluid Mechanics and Heat Transfer. Taylor & Francis.
- Eymard, R. Gallouët, T. R. Herbin, R. (2000) The finite volume method Handbook of Numerical Analysis, Vol. VII, 2000, p. 713–1020. Editors: P.G. Ciarlet and J.L. Lions.
- LeVeque, Randall (2002), Finite Volume Methods for Hyperbolic Problems, Cambridge University Press.
- Toro, E. F. (1999), Riemann Solvers and Numerical Methods for Fluid Dynamics, Springer-Verlag.
- Patankar, Suhas V. (1980), Numerical Heat Transfer and Fluid Flow, Hemisphere.
- Hirsch, C. (1990), Numerical Computation of Internal and External Flows, Volume 2: Computational Methods for Inviscid and Viscous Flows, Wiley.
- Laney, Culbert B. (1998), Computational Gas Dynamics, Cambridge University Press.
- LeVeque, Randall (1990), Numerical Methods for Conservation Laws, ETH Lectures in Mathematics Series, Birkhauser-Verlag.
- Tannehill, John C., et al., (1997), Computational Fluid mechanics and Heat Transfer, 2nd Ed., Taylor and Francis.
- Wesseling, Pieter (2001), Principles of Computational Fluid Dynamics, Springer-Verlag.
- Finite volume methods by R. Eymard, T Gallouët and R. Herbin, update of the article published in Handbook of Numerical Analysis, 2000
- Rübenkönig, Oliver. “The Finite Volume Method (FVM) – An introduction”. Archived from the original on 2009-10-02., available under the GFDL.
- FiPy: A Finite Volume PDE Solver Using Python from NIST.
- CLAWPACK: a software package designed to compute numerical solutions to hyperbolic partial differential equations using a wave propagation approach