INTRODUCTION TO FINITE VOLUME METHOD FOR COMPUTATIONAL FLUID DYNAMICS I - SOLVING  TWO DIMENSIONAL CONDUCTION PROBLEMS

This notebook has been written in Mathematica by 
Alexandre M. S. Costa
Paulo H. Mareze
Aluizio B. Drugovich
Maringa State University
Maringa PR
Brazil

desouzadacosta@yahoo.com
http://www.geocities.com/desouzadacosta


Version:  23/8/05
More recent versions of this notebook should be available at the web site:
http://www.geocities.com/desouzadacosta/wphd/stabilityfluidiz.zip

Making use of symbolic and numeric capabilities of Mathematica, in this notebook we explore the fundamentals of the finite volume method (FVM). This method is largely employed for solution of computational fluid dynamics (CFD) problems in engineering. Nowadays, There are many commercial CFD packages available. However, for efficient use of these packages, the knowledge of the physics and mathematics behind their solutions are important In this notebook, we discuss the following concepts related to the FVM: regular mesh construction, assembling coefficient matrix, implicit and explicit formulation, and treatment of boundary condition types. Finally, we apply all concepts for solving two typical conduction problems.

The application of mathematical modeling and computational simulation to processes and equipment has been almost fundamental in the development of these. Among the computational simulation tools, we detach the CFD (Computational Fluid Dynamics). This tool allows the quantitative prediction of the behavior of equipment, structures and natural phenomena. This prediction are obtained by the following steps:
·	Expressing the relevant physical and chemical laws, besides the required complementary models, in the form of equations which relates pressure, velocity, temperature and concentration values, among others, to points distributed in space and time.
·	Placing the points (i.e., the computational mesh) very close to each other in order to satisfy the continuum hypothesis.
·	Solving the equations systematically and iteratively using error reducing methods
·	Finishing the calculations when convergence is obtained, i.e., when the errors are sufficiently small.

Finally, CFD programs can manipulate a large variety of physical processes in circumstances such as:
·	Fluid flow (e.g. steam, combustion gases, oil),
·	Heat transfer,
·	Phase change,
·	Other simultaneous effects (e.g. thermal stresses, chemical reactions)

The finite volume method is a CFD method for solution of the set of differential equations that describes mass, momentum and energy conservation. The  FVM is related to other approximated numerical methods such as finite differences and finite elements (vide Chung, for a discussion).
The finite volume method is characterized by the partition of the spatial domain in a finite number of elementary volumes for which are applied the balances of mass, momentum and energy. The approximated equations for the method can be obtained by two approaches. The first consists in applying balances for the elementary volumes (finite volumes), and the second consists in the integration spatial-temporal of the conservation equations. In this work, we are going to follow the second approach.
The general form of the differential conservation equation written in a cartesian reference frame to be treated by the FVM is:

\(\(\(\(∂\/∂t\) \((ρϕ)\) + \(∂\/∂x\) \((ρuϕ)\) + \(∂\/∂y\) \((ρvϕ)\) + \(∂\/∂z\) \((ρwϕ)\) = \(∂\/∂x\) \((\(Γ\^ϕ\) ∂ϕ\/∂x)\) + \(∂\/∂y\) \((\(Γ\^ϕ\) ∂ϕ\/∂y)\) + \(∂\/∂z\) \((\(Γ\^ϕ\) ∂ϕ\/∂z)\) + S\^ϕ\)\(
\)
    \)\)

For two dimensional flow in the xy plane , incompressible fluid (density ρ constant), newtonian and constant properties, the variables for Eq. (1) are defined in the following table:

Values for ϕ,  (Γ^ϕ),  (S^ϕ)
                                             ({{\(\(\ \ \ \ \ \)\(\(Balance\)\(
\)
              \(\ \ \ \ \ \ \)\(Equation\)\(\ \ \)\)\), ϕ, (Γ^ϕ), (S^ϕ)}, {Mass, 1, 0, 0}, {x Momentum, u, μ, -∂(p/∂)x}, {y Momentum, v, μ, -∂(p/∂)y}, {Thermal Energy, T, (k/(Subscript[c, p])), (μ/(Subscript[c, p])) Φ}})

u, v, p and T are velocity components for x and y directions, pressure, and temperature respectively. μ , k and cp stands for viscosity, thermal conductivity and specific heat at constant pressure, respectively. We neglect the gravity force and thermal energy generation by viscosity Φ. Integrating the Eq. (1) over a cubic volume of dimensions  Δx, Δy, Δz (see Fig. (1)) and over time (interval between t and t + Δt):

\(\(\(∫\+\(V, t\)\(∂\/∂t\) \((ρϕ)\) dVdt + ∫\+\(V, t\)\(∂\/∂x\) \((ρuϕ)\) 
            dVdt \(∫\+\(V, t\)\(∂\/∂y\) \((ρvϕ)\) 
                dVdt\) = ∫\+\(V, t\)\(∂\/∂x\) \((\(Γ\^ϕ\) ∂ϕ\/∂x)\) 
          dVdt \(∫\+\(V, t\)\(∂\/∂y\) \((\(Γ\^ϕ\) ∂ϕ\/∂y)\) 
              dVdt \(∫\+\(V, t\)\(S\^ϕ\) dVdt\)\)\)\(
\)
    \)\)

Results

\(\(\(\(\(M\_P\) ϕ\_p - \(M\_p\%o\) ϕ\_p\%o\)\/Δt + \(M\& . \_e\) ϕ\_e\%θ - \(M\& . \_w\) ϕ\_w\%θ + \(M\& . \_n\) ϕ\_n\%θ - \(M\& . \_s\) ϕ\_s\%θ = 
\(D\_1[∂ϕ\/∂x]\)\_e\%θ - \(D\_1[∂ϕ\/∂x]\)\_w\%θ + \(D\_2[∂ϕ\/∂y]\)\_n\%θ - \(D\_2[∂ϕ\/∂y]\)\_s\%θ + 
        L[S\^ϕ]\^θ\)\(
\)
    \)\)

The process of integration described by Eqs. (2) and (3) is named discretization. The equation  (3) and all similar equations obtained by this type of process are known as discretized equations.
From Eq. (3)

\(M\_p = ρΔxΔyΔz\ \ \ \ \ \ \ ; \ \ \ \ \ M\_p = \(ρ\^o\) ΔxΔyΔz\)

\(M\& . \_e = \([ρuΔyΔz]\)\_e\ \ \ \ ; \ \ \ \ M\& . \_w = \([ρuΔyΔz]\)\_w\)

\(M\& . \_n = \([ρvΔxΔz]\)\_n\ \ \ \ ; \ \ \ \ M\& . \_s = \([ρvΔxΔz]\)\_s\)

\(D\_\(1  e\) = \([\(Γ\^ϕ\) ΔyΔz]\)\_e\ \ \ ; \ \ D\_\(1  w\) = \([\(Γ\^ϕ\) ΔyΔz]\)\_w\)

\(D\_\(2  n\) = \([\(Γ\^ϕ\) ΔxΔz]\)\_n\ \ \ ; \ \ D\_\(2  s\) = \([\(Γ\^ϕ\) ΔxΔz]\)\_s\)

From the previous equations the “o” superscript stands for the time t and no superscript stands for the time t + Δt.The superscript “θ” stands for a time between t and t+Δt.In this time the ϕ variable is weighted by θ value (between 0 and 1):

\(\(\(\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \)\(ϕ\^θ\  = \ θ\ ϕ\  + \ \((1\  - \ θ)\) ϕ\^o\)\)\)

Yet, referring to Eq. (3), the volume face values are weighted by the adjacent volume values:

\(ϕ\_e = \((1\/2 + α\_e)\) ϕ\_p + \((1\/2 - α\_e)\) ϕ\_E\)

\(ϕ\_w = \((1\/2 + α\_w)\) ϕ\_W + \((1\/2 - α\_w)\) ϕ\_p\)

\(ϕ\_n = \((1\/2 + α\_n)\) ϕ\_p + \((1\/2 - α\_n)\) ϕ\_N\)

\(\(\(ϕ\_s = \((1\/2 + α\_s)\) ϕ\_S + \((1\/2 - α\_s)\) ϕ\_p\)\(
\)
    \)\)

Similarly, the face gradient values for ϕ are weighted by the β variable:

\(\([∂ϕ\/∂x]\)\_e = \(β\_e\) \((\(ϕ\_E - ϕ\_p\)\/Δx\_e)\)\)

\(\([∂ϕ\/∂x]\)\_w = \(β\_w\) \((\(ϕ\_p - ϕ\_W\)\/Δx\_w)\)\)

\(\([∂ϕ\/∂x]\)\_n = \(β\_n\) \((\(ϕ\_N - ϕ\_p\)\/Δy\_n)\)\)

\(\([∂ϕ\/∂x]\)\_s = \(β\_s\) \((\(ϕ\_p - ϕ\_S\)\/Δy\_s)\)\)

The case for which α = 0 and β = 1 is named central difference scheme.This scheme is commonly employed for pure conduction problems.The case for which β = 0 and α = 0,5 or α = -0,5 is named upwind scheme for positive and negative velocities, respectively.
The last term from equation  corresponds to the last integral on equation . We are going to assume that:

L[(S^ϕ)] = [(Subscript[S, P]) ((Subscript[ϕ, P])^o)  +  
            (Subscript[S, c])] ΔV

Where 
(Subscript[S, P]) and  (Subscript[S, c]) are coefficients defined according to ϕ
ΔV is the volume Δx Δy Δz

Using eqs (4) to (17) and the discretized mass conservation equation (underlined below), equation (3) becomes:

\(\(\(M\_P\) ϕ\_p\)\/Δt + ϕ\_p\%θ[\(M\& . \_e\) \((1\/2 + α\_e)\) - \(M\& . \_w\) \((1\/2 - α\_w)\) + \(M\& . \_n\) \((1\/2 + α\_n)\) - \(M\& . \_s\) \((1\/2 - α\_s)\) - \(S\_p\) ΔV + \([\(\(D\_1\) β\)\/Δx]\)\_e + \([\(\(D\_1\) β\)\/Δx]\)\_w + \([\(\(D\_2\) β\)\/Δy]\)\_n + \([\(\(D\_2\) β\)\/Δy]\)\_s - \(\(M\_P - M\_p\%o\)\/Δt - M\& . \_e + M\& . \_w - M\& . \_n + M\& . \_s\)\+_] = ϕ\_E\%θ[\(-\((1\/2 - α\_e)\)\) 
              M\& . \_e + \([\(\(D\_1\) β\)\/Δx]\)\_e] + ϕ\_W\%θ[\((1\/2 + α\_w)\) 
              M\& . \_w + \([\(\(D\_1\) β\)\/Δx]\)\_w] + ϕ\_N\%θ[\(-\((1\/2 - α\_n)\)\) 
              M\& . \_n + \([\(\(D\_2\) β\)\/Δy]\)\_n] + ϕ\_S\%θ[\((1\/2 + α\_s)\) 
              M\& . \_s + \([\(\(D\_2\) β\)\/Δy]\)\_s] + \(\(M\_p\%o\) ϕ\_p\%o\)\/Δt + \(S\_c\) ΔV\)

Collecting the ϕ coefficients:

\(\(\(M\_P\) ϕ\_p\)\/Δt + \(A\_p\%*\) ϕ\_p\%θ = \(A\_e\) ϕ\_E\%θ + \(A\_w\) ϕ\_W\%θ + \(A\_n\) ϕ\_N\%θ + \(A\_s\) ϕ\_S\%θ + \(\(M\_p\%o\) ϕ\_p\%o\)\/Δt + \(S\_c\) ΔV\)

Where

\(\(A\_p\%*\) = 
      A\_e + A\_w + A\_n + A\_s - \(S\_p\) ΔV - 
        M\_P\/Δt + 
        M\_p\%o\/Δt\)

\(A\_e = \(-\((1\/2 - α\_e)\)\) 
          M\& . \_e + \([\(\(D\_1\) β\)\/Δx]\)\_e\)

\(A\_w = \((1\/2 + α\_w)\) 
          M\& . \_w + \([\(\(D\_1\) β\)\/Δx]\)\_w\)

\(A\_n = \(-\((1\/2 - α\_n)\)\) 
          M\& . \_n + \([\(\(D\_2\) β\)\/Δy]\)\_n\)

\(A\_s = \((1\/2 + α\_s)\) 
          M\& . \_s + \([\(\(D\_2\) β\)\/Δy]\)\_s\)

The equation (20) for θ = 0 (explicit formulation) becomes:

\(\(\(M\_P\) ϕ\_p\)\/Δt = \(ϕ\_p\%o\) \((M\_p\%o\/Δt - \(A\_p\%*\))\) + ∑\(A\_nb\) ϕ\_NB\%o + \(S\_c\) ΔV\)

And for  θ = 1 (fully implicit):

\(\(A\_p\) ϕ\_p = ∑\(A\_nb\) ϕ\_NB\ \  + \ \(\(M\_p\%o\) ϕ\_p\%o\)\/Δt + \(S\_c\) ΔV\)

Where

(Subscript[A, p]) = ∑(Subscript[A, nb]) - (Subscript[S, p]) ΔV + 
        ((Subscript[M, p]^o)/Δ)t

In this section we illustrate the FVM applied to the solution of a two dimensional transient conduction problem without thermal energy sources.
Referring to Table (1), Eq. (20) and approximating the face volume derivatives by central differences,  Eqs (26) and (27)  becomes:

For the explicit formulation:

\(\(\(\(A\_p\) 
        T\_p = \(A\_e\) T\_E\%o + \(A\_w\) T\_W\%o + \(A\_n\) 
          T\_N\%o + \(A\_s\) 
          T\_S\%o + \((A\_p\%o - A\_e - A\_w - A\_n - 
              A\_s + \(S\_p\) ΔxΔy)\) 
          T\_p\%o + \(S\_c\) ΔxΔy\)\(
\)
    \)\)

Where

(Subscript[A, p]) = (Subscript[M, p])/Δt   ;     (Subscript[A, p]^o) = 
      (Subscript[M, p]^o)/Δt

For the fully implicit formulation

(Subscript[A, p]) 
        (Subscript[T, p]) = (Subscript[A, e]) (Subscript[T, E]) + (Subscript[A, w]) (Subscript[T, W]) + (Subscript[A, n]) (Subscript[T, N]) + (Subscript[A, s]) 
          (Subscript[T, S]) + (Subscript[B, p]^T)

Where

\(A\_e = \(\([k\/\(\(c\_p\) Δx\)]\)\_e\) Δy\ \ \ ; \ A\_w = \(\([k\/\(\(c\_p\) Δx\)]\)\_w\) Δy\ \ \ ; \ \ A\_n = \(\([k\/\(\(c\_p\) Δy\)]\)\_n\) Δx\ \ \ ; 
    A\_s = \(\([k\/\(\(c\_p\) Δy\)]\)\_s\) Δx\)

\(B\_p\%T = \(A\_p\%o\) T\_p\%o + 
        ScΔxΔy\ \ \ ; \ \ \ A\_p\%o = \(M\_p\%o\/Δt = \(ρ\^o\/Δt\) ΔxΔy\)\)

\(A\_p = 
      A\_e + A\_w + A\_n + A\_s + 
        M\_P\/Δt - \(S\_p\) ΔxΔy\)

Writing the equation (28) for each center point of each volume results in a linear set of equations given by:

\(\(\(\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \)\(\(A\ T = B\)\(
\)
    \)\)\)

where:
 A is the coefficient matrix
 T is the vector of temperatures
  B is the independent vector.

In the next subsection we are going to show the construction of matrix A and vector B.

In this section, we illustrate the structure of the FVM coefficient matrix for conduction in a two dimensional retangular domain

(* The next cell loads the LinearAlgebra`MatrixManipulation` package . 
This package includes functions for manipulating and constructing a variety of  matrices. For evaluating, select it (click the mouse in the middle of the cell bracket on the right corner) and press SHIFT + ENTER *)

<< LinearAlgebra`MatrixManipulation`

(*The two dimensional domain is going to be partitioned into 16 volumes :  
              4 in the x direction and 4 in the y direction .  
                  In the next cell we create the matrix for representing the partitioned domain (named here matrixGrid) *)

MatrixForm[Array[matrixGrid, {4, 4}]]

(* Next cell sets the matrix matrixGrid with elements .  Each one of these elements corresponds to numbers from 1 to 16 from the upper right to the bottom right corner *)

k  =  0 ;
Do[  
 
        Do[
      k  =  k  + 1 ; 
          matrixGrid[i, j]  =  k; , {j, 1, 
            4}]
                               , {i, 1, 4}];
MatrixForm[Table[matrixGrid[i, j], {i, 1, 4}, {j, 1, 4}]]

(* In the next cell, we build the coefficients matrix A ( see equation 3  5). This matrix is named here matrixT. The matrix elements are obtained by sweeping of MatrixGrid .  The sweep is firstly by column and after that by row. *)

matrixT  =  ZeroMatrix[16]; (* First all the elements of the coefficients matrix are zeros .  The function ZeroMatrix creates an 16 by 16 matrix of zeros *)

Do[  
 Do[
      If[j + 1 ≤  4, 
                  matrixT[[matrixGrid[i, j], 
                        matrixGrid[i, j + 1]]]  =  -Ae];
      If[j - 1  ≥   1, 
                  matrixT[[matrixGrid[i, j], 
                        matrixGrid[i, j - 1]]]  =  -Aw];
      If[i + 1  ≤    4, 
                    matrixT[[matrixGrid[i, j], 
                          matrixGrid[i + 1, j]]]  =  -As]; 
      If[i - 1  ≥   1, 
                    matrixT[[matrixGrid[i, j], 
                          matrixGrid[i - 1, j]]]  =  -An];
      matrixT[[matrixGrid[i, j], 
                      matrixGrid[i, j]]]  =  Ap;
,{j, 1, 4}]
                               ,{i, 1, 4}];

MatrixForm[matrixT]

In this section we discuss one method for treating boundary conditions. This method is named balance for boundary volume. 
Initially, we use the Fourier law applied to the west face of a boundary volume in a conduction problem:

\(\([q\_f\%\(\*"""\)]\)\_w = \(k[∂T\/∂x]\)\_w\%o\)

This term is analogous to term “(Subscript[((Subscript[D, 1])[∂(ϕ/∂)x]), w]^θ) ” in equation  (3). The following table summarizes the three types of boundary conditions and corresponding expression for "\(\(\([\)\(q\_f\%\(\*"""\)\)\(]\)\)\_w\)":

\(Types\ of\ boundary\ conditions\ for\ Thermal\ Energy\ balance\)                            
              ({{□, \(\(q\_f\%''\)\( | \_w\)\)}, {Prescribed Temperature, \(\(k\_f\) \(T\_f - T\_p\)\/Δx\_f\)}, {Prescribed flux, \(\(q\_f\%''\)\( | \_w\)\)}, {Convection, (h/(1 + 
                      hΔ((Subscript[x, f])/(Subscript[k, f])))) ((Subscript[T, ∞]) - 
                    (Subscript[T, p]))}})

In the above table, (Subscript[T, f])  is the boundary face temperature and (Subscript[T, p]) is the center temperature for the boundary volume.
It should be noticed that the third line of this table is obtained by equaling the convection flux to the conduction flux:

\(q\_f\%\(\*"""\) = \(h \((T\_∞ - 
              T\_f)\) = \(k\_f\) \(T\_f - T\_p\)\/Δx\_f\)\)

Here,  h is the convection heat transfer coefficient and (Subscript[T, ∞]) is the outside temperature. Isolating (Subscript[T, f]) from the above equation we can obtain the value for \(\(\([\)\(q\_f\%\(\*"""\)\)\(]\)\)\_w\)

We illustrate the method using the next figure:

Applying the thermal energy balance for the volume centered on P, using equation (3) and Table 1, and the central difference scheme, results:

\(\(\(\(\(M\_P\) ϕ\_p - \(M\_p\%o\) ϕ\_p\%o\)\/Δt\ \  = \(\(\(2  k\_f\)\/c\_p\)[\((\(T\_f1 - T\_p\)\/Δx)\)]\)\_w - \ \ \([\(k\/c\_p\) \((\(T\_p - T\_E\)\/Δx)\)]\)\_e + \([\(\(h/c\_p\)\/\(1 + hΔy/2  k\_f\)\) \((T\_∞ - T\_p)\)]\)\_s - \([\ \ \(k\/c\_p\) \((\(T\_p - T\_N\)\/Δy)\)]\)\_n\)\(
\)
    \)\)

And the equation analogous to Eq. (31)  becomes:

(Subscript[A, p]) (Subscript[T, p]) = (Subscript[A, e]) (Subscript[T, E]) + (Subscript[A, n]) (Subscript[T, N]) + 
        (Subscript[B, p])

Where

\(B\_p = \(A\_p\%o\) 
          T\_p\%o + \(\([\(2  k\)\/\(\(c\_p\) Δx\)]\)\_w\) 
          T\_f1\  + \(\([\(h/c\_p\)\/\(1 + hΔy/2  k\_f\)]\)\_s\) 
          T\_∞\ \ ; \ \ A\_p\%o = \(M\_p\%o\/Δt = \(ρ\^o\/Δt\) ΔxΔy\)\)

\(A\_p = 
      A\_e + A\_n + \([\(2  k\)\/\(\(c\_p\) Δx\)]\)\_w + \([\(h/c\_p\)\/\(1 + hΔy/2  k\_f\)]\)\_s + 
        M\_p\/Δt\)

matrixT  =  ZeroMatrix[16];
Do[  
 Do[
      matrixT[[matrixGrid[i, j], 
                      matrixGrid[i, j]]]  =  Ap;  ArrayB[matrixGrid[i, j]]    =  Bp ;
      If[j + 1 ≤  4,matrixT[[matrixGrid[i, j], 
                          matrixGrid[i, j + 1]]] = -Ae,
             matrixT[[matrixGrid[i, j], 
                          matrixGrid[i, 
                            j]]]  =  
           matrixT[[matrixGrid[i, j], matrixGrid[i, j]]]  -  
                        Ae   +   ((2 k)/Δ)x; 
           ArrayB[
                        matrixGrid[i, j]]    =    ArrayB[
                          matrixGrid[i, 
                            j]]  +   ((2 k Tf)/Δ)x  ];

      If[  j - 1  ≥   1,matrixT[[matrixGrid[i, j], 
                          matrixGrid[i, j - 1]]] = -Aw,
          matrixT[[matrixGrid[i, j], 
                          matrixGrid[i, 
                            j]]]  =  
           matrixT[[matrixGrid[i, j], matrixGrid[i, j]]]  -  
                        Aw   +   ((2 k)/Δ)x; 
           ArrayB[
                        matrixGrid[i, j]]    =    ArrayB[
                          matrixGrid[i, 
                            j]]  +   ((2 k Tf)/Δ)x     ];

      If[i + 1  ≤    4,matrixT[[matrixGrid[i, j], 
                          matrixGrid[i + 1, j]]] = -As,
           matrixT[[matrixGrid[i, j], 
                          matrixGrid[i, 
                            j]]]  =  
           matrixT[[matrixGrid[i, j], matrixGrid[i, j]]]  -  
                        As   +   (h/(1 +  (h/k) Δ(x/2))); 
           ArrayB[
                        matrixGrid[i, j]]    =   ArrayB[
                          matrixGrid[i, 
                            j]]  +  (h/(1 +  (h/k) Δ(x/2))) T∞ ]; 

      If[i - 1  ≥   1,matrixT[[matrixGrid[i, j], 
                          matrixGrid[i - 1, j]]] = -An,
             matrixT[[matrixGrid[i, j], 
                          matrixGrid[i, 
                            j]]]  =  
           matrixT[[matrixGrid[i, j], matrixGrid[i, j]]]  -  
                        An   +   ((2 k/Cp)/Δ)y; 
           ArrayB[
                        matrixGrid[i, j]]    =   ArrayB[
                          matrixGrid[i, 
                            j]]  + ((2 k Tf/Cp)/Δ)y  ];     
,{j, 1, 4}]
                               ,{i, 1, 4}];

MatrixForm[matrixT]

v  =  
        Table[ArrayB[o], {o, 1, Nx Ny}];
MatrixForm[v]

In this section we use the FVM to solve two conduction problems.
The first one is based on the solved example 4.3 (Incropera). It consists on a two dimensional steady state conduction problem in a infinite slab. The slab has a rectangular cross section with 1 m × 1 m (see figure below). The thermal conductivity is 1 W/ m.K. The lower side is submitted to convection with air at 300 K and heat transfer coefficient of 10 W/ m2.K. All the other sides are kept at 500 K.

The solution procedure follows the same steps shown in the previous section. The Mathematica lines are shown on the next cell

<< LinearAlgebra`MatrixManipulation`

\(Nx\  = 20; \ \ \  (*\ 
      number\ of\ volumes\ in\ the\ x\ direction\ *) 
Ny\  = \ 20; \ \ \  (*\ 
      number\ of\ volumes\ in\ the\ y\ direction\ *) 
Δy = 1/Ny; \  (*\ 
      lenght\ of\ y\ side\ of\ volume\ *) 
Δx = 1/Nx;\ \  (*\ 
      lenght\ of\ x\ side\ of\ volume\ *) \)
Δz   =  1;
MatrixForm[Array[matrixGrid, {Ny, Nx}]]
MatrixForm[Array[ArrayB, {Nx*Ny}]]
k  =  0 ;
Do[  
 
        Do[
      k  =  k  + 1 ; 
          matrixGrid[i, j]  =  k; , {j, 1, 
            Nx}]
                               , {i, 1, Ny}];
ClearAll[k];

Cp = 1    
k = 1      
  
Mp = 0;     
Apo  = 0;     
  
Sp = 0;    
Sc = 0 ;   

h = 10;   
Tf =  500;   
T∞ = 300;   
  
Ae  =  (k/(Cp Δx)) Δy Δz ;
Aw  =  (k/(Cp Δx)) Δy Δz ;
An  =  (k/(Cp Δy)) Δx Δz ;
As   =  (k/(Cp Δy)) Δx Δz ;
Ap = Ae + Aw + An + As + (Mp/Δ)t-Sp Δx Δy Δz ;
Bp  =  
            Apo Tpo  +  
              Sc Δx Δy Δz ;

matrixT  =  ZeroMatrix[Nx*Ny];
Do[  
 Do[
      matrixT[[matrixGrid[i, j], 
                            matrixGrid[i, j]]]  =  Ap;           ArrayB[matrixGrid[i, j]]    =  Bp ;
      If[j + 1 ≤  Nx, 
                        matrixT[[matrixGrid[i, j], 
                              matrixGrid[i, 
                                j + 1]]]  =  -Ae, 
             matrixT[[matrixGrid[i, j], 
                              matrixGrid[i, 
                                
                                j]]]  =  
           matrixT[[matrixGrid[i, j], matrixGrid[i, j]]]  -  
                            Ae   +   ((2 k/Cp)/Δ)x; 
           ArrayB[
                            matrixGrid[i, j]]    =    ArrayB[
                              matrixGrid[i, 
                                j]]  +   ((2 k Tf/Cp)/Δ)x  ];

      If[  j - 1  ≥   1, 
                        matrixT[[matrixGrid[i, j], 
                              matrixGrid[i, 
                                j - 1]]]  =  -Aw, 
          matrixT[[matrixGrid[i, j], 
                              matrixGrid[i, 
                                j]]]  =  
           matrixT[[matrixGrid[i, j], matrixGrid[i, j]]]  -  
                            Aw   +   ((2 k/Cp)/Δ)x; 
           ArrayB[
                            matrixGrid[i, j]]    =    ArrayB[
                              matrixGrid[i, 
                                j]]  +   ((2 k Tf/Cp)/Δ)x     ];
      If[i + 1  ≤   Ny,matrixT[[matrixGrid[i, j], 
                                matrixGrid[i + 1, j]]]  =  -As,
           matrixT[[matrixGrid[i, j], 
                                  matrixGrid[i, j]]] = 
           matrixT[[matrixGrid[i, j], 
                                    matrixGrid[i, j]]] - As  +((h/
                                        Cp)/(1 +  (h/k) Δ(y/2)));
           ArrayB[matrixGrid[i, j]]   =  ArrayB[matrixGrid[i, j]] +  ((h/
                                        Cp)/(1 +  (h/k) Δ(y/2))) T∞   ]; 

      If[i - 1  ≥   1,matrixT[[matrixGrid[i, j], 
                              matrixGrid[i - 1, j]]]  =  -An,
             matrixT[[matrixGrid[i, j], 
                                matrixGrid[i, j]]] = 
           matrixT[[matrixGrid[i, j], 
                                  matrixGrid[i, j]]] - An  +  ((2 k/Cp)/Δ)y;
           ArrayB[matrixGrid[i, j]]   =  ArrayB[matrixGrid[i, j]] +((2 k Tf/Cp)/Δ)y   ];
    
,{j, 1, Nx}]
                               ,{i, 1, Ny}];

v  =  Table[ArrayB[o], {o, 1, Nx Ny}]

(*For solving the linear equation set we use the function LinearSolve*)

sol=N[LinearSolve[matrixT,v]]

s  =  Partition[sol, 
          Nx]//MatrixForm

ShowLegend[ListContourPlot[s,ColorFunction->(Hue[.75- .75  #]&),DisplayFunction -> Identity],{(Hue[ .75  #] &),12,ToString[IntegerPart[Max[s]]],ToString[IntegerPart[Min[s]]],LegendPosition -> {1.1, - .4}}]

The two dimensional unsteady conduction problem consists on slab of 1 m × 1 m of cross-sectional area. The thermal difusivity (= (k/(ρ (Subscript[c, p])))) is specified. Initially, all the slab is at 100 C. After that all the faces are submited to a prescribed temperature condition of 0 C (see Figure below). We are going to solve this problem using time steps of Δt = 0.1 s to the final time of 1 s.

The solution procedure follows the same steps shown in the previous sections. The Mathematica lines are shown on the next cell

<< LinearAlgebra`MatrixManipulation`

    << Graphics`Legend`

    Off[General::spell];  Off[General::spell1];

Nx  = 30;    
Ny  = 30;   

  
MatrixForm[Array[matrixGrid, {Ny, Nx}]];
k  =  0 ;
Do[  
 
          Do[
      k  =  k  + 1 ; 
            matrixGrid[i, j]  =  k; , {j, 1, 
              Nx}]
                               , {i, 1, Ny}];
ClearAll[k];

  

Δt = 0.1;   
it = 10;   

  
Δy = 1/Ny;
Δx = 1/Nx;
\(\(Δz = 1;\)\(
\)\)
\(Sp = 
        0; \  (*No\ thermal\ energy\ sources\  - \ see\ Eq . \ 18*) 

      Sc = 0;  (*No\ thermal\ energy\ sources\  - \ see\ Eq . \ 18*) 

 (*\ 
        Side\ temperatures\ in\ Celsius\ *) 

      Tf1 = 0; \  (*North\ face*) 
Tf2 = 
        0; \  (*East\ face*) 
Tf3 = 
        0; \  (*South\ face*) 
Tf4 = 
        0; \  (*West\ face*) 

 (*Coefficients\ for\ the\ fully\ implicit\ discretization, \ see\ Eq . \ 28\ to\ 31*) \ 
α = 1\/π\^2;  (*\ 
        Thermal\ difusivity\ *) 
 (*\ 
        Note : \ 
            For\ comparison\ purposes\ with\ the\ analytical\ solution\ we\ are\ introducing\ the\ thermal\ difusivity\ α\  = \ \(k\/ρc\_p\)\  . From\ this\ definition\ the\ coefficients\ corresponding\ to\ Eqs . \ 
              
              28\ to\ 31\ are\ given\ next\ *) 
Ae\  = \ \(α\/\(\(\ \)\(Δx\)\)\) Δy\ Δz\ ; \)
Aw  =   (α/( Δx)) Δy Δz ;
An  =  (α/( Δy)) Δx Δz ;
As   =   (α/( Δy)) Δx Δz ;
Ap  =  
          Ae  +  Aw  +  An  +  
            As  +  ((Δx Δy Δz)/Δ)t - 
            Sp Δx Δy Δz ;
Apo  = Δx Δy Δz/Δt;
 (*Initial temperature field*) 
Tpo = Table[100, {Nx*Ny}];
G = {};
Do[
matrixT  =  ZeroMatrix[Nx*Ny];
Array[ArrayB, {Nx*Ny}];
Do[  
 Do[
      matrixT[[matrixGrid[i, j], 
                                matrixGrid[i, j]]]  =  Ap;  ArrayB [ matrixGrid[i, j] ]     =  
                            Apo*Tpo[[   matrixGrid[i, j] ]  ] ;

      If[j + 1 ≤  Nx,matrixT[[matrixGrid[i, j], 
                                    matrixGrid[i, j + 1]]]  =  -Ae,
             matrixT[[matrixGrid[i, j], 
                                      matrixGrid[i, 
                                        j]]]  =  
           matrixT[[matrixGrid[i, j], matrixGrid[i, j]]]  -  
                                    Ae   + ((2*α)/( Δx))  ;
           ArrayB[matrixGrid[i, j]]   =   ArrayB[
                                      matrixGrid[i, 
                                        j]]  +   ((2*α*Tf2)/( Δx))  ];  


      If[  j - 1  ≥   1,matrixT[[matrixGrid[i, j], 
                                    matrixGrid[i, j - 1]]]  =  -Aw,
          matrixT[[matrixGrid[i, j], 
                                      matrixGrid[i, 
                                        j]]]  =  
           matrixT[[matrixGrid[i, j], matrixGrid[i, j]]]  -  
                                    Aw   +   ((2*α)/( Δx));
           ArrayB[
                                    matrixGrid[i, j]]    =    ArrayB[
                                      matrixGrid[i, 
                                        j]]  + ((2*α*Tf4)/( Δx))     ];
  


      If[i + 1  ≤   Ny,matrixT[[matrixGrid[i, j], 
                                    matrixGrid[i + 1, j]]]  =  -An,
           matrixT[[matrixGrid[i, j], 
                                      matrixGrid[i, 
                                        j]]]  =  
           matrixT[[matrixGrid[i, j], matrixGrid[i, j]]]  -  
                                    An   +   ((2*α)/( Δy)) ;
           ArrayB[
                                    matrixGrid[i, j]]    =   ArrayB[
                                      matrixGrid[i, 
                                        j]]  + ((2*α*Tf1)/( Δy))  ];
  
 

      If[i - 1  ≥   1,matrixT[[matrixGrid[i, j], 
                                    matrixGrid[i - 1, j]]]  =  -As,
             matrixT[[matrixGrid[i, j], 
                                      matrixGrid[i, 
                                        j]]]  =  
           matrixT[[matrixGrid[i, j], matrixGrid[i, j]]]  -  
                                    As   +  ((2*α)/( Δy));
           ArrayB[
                                    matrixGrid[i, j]]    =   ArrayB[
                                      matrixGrid[i, 
                                        j]]  +  ((2*α*Tf3)/( Δy))    ]; 
  
    
,{j, 1, Nx}]; 
 ,{i, 1, Ny}];
v  =  Table[ArrayB[o], {o, 1, Nx Ny}];
ClearAll[Tpo, ArrayB];
  
 Tpo = LinearSolve[matrixT, v];
   
G = Append[G, Partition[  N[Tpo] , Nx]];
,{n, 1, it}];

ArrayframeN  =  { };     
Do[ArrayframeN = Append[ArrayframeN,
 ShowLegend[ListDensityPlot[G[[vet]],ColorFunction->(Hue[.75- .75  #]&),PlotRange -> {0, 100},DisplayFunction -> Identity],{(Hue[ .75  #] &),10,ToString[IntegerPart[Max[G[[vet]]]]],ToString[IntegerPart[Min[G[[vet]]]]],LegendPosition -> {1.1, - .4}}]             ],{vet, 1, it}];
\(\(SelectionMove[EvaluationNotebook[], \ All, \ 
        GeneratedCell]\)\(;\)\(\ \)\(FrontEndTokenExecute["CellGroup"]\)\(;\)\(\ \)\(FrontEndTokenExecute["OpenCloseGroup"]\)\(;\)\(\ \)\)

(* 
      Next cell allows saving the  results as a animated gif *)

SetDirectory["C:\FolderName"];   
Export["NameGifFile.gif",  Arrayframe  N,  
          ConversionOptions -> {"Loop" -> True}] ;

ArrayframeN2  =  { };   
Do[ArrayframeN2 = Append[ArrayframeN2,
 ShowLegend[ListContourPlot[G[[vet]],ColorFunction->(Hue[.75- .75  #]&),PlotRange -> {0, 100},DisplayFunction -> Identity],{(Hue[ .75  #] &),10,ToString[IntegerPart[Max[G[[vet]]]]],ToString[IntegerPart[Min[G[[vet]]]]],LegendPosition -> {1.1, - .4}}]             ],{vet, 1, it}];
\(\(SelectionMove[EvaluationNotebook[], \ All, \ 
        GeneratedCell]\)\(;\)\(\ \)\(FrontEndTokenExecute["CellGroup"]\)\(;\)\(\ \)\(FrontEndTokenExecute["OpenCloseGroup"]\)\(;\)\(\ \)\)

mmax  =  10 ;
nmax  =  10 ;
T  =  (400/(π^2)) (∑_(n = 1)^nmax)(∑_(m = 1)^mmax)(((1 - (((-1))^m)) (1 - (((-1))^n)))/(m n)) Exp[-((m^2) + (n^2)) t] Sin[m π x] 
                Sin[n π y];

Arrayframe  =  { };   
Do[    Arrayframe = Append[Arrayframe,ShowLegend[ DensityPlot[T,{x, 0, 1},  {y, 0, 1},PlotPoints -> 30,PlotRange -> {0, 100},ColorFunction->(Hue[.75- .75  #]&),DisplayFunction -> Identity] ,{(Hue[ .75  #] &),10, ToString[x = 0.5; y = 0.5; 
                          IntegerPart[N[T]]],"0",LegendPosition -> {1.1, - .4}}]  ]             ,{t, 0.1, 1, 0.1}];
SelectionMove[EvaluationNotebook[],  All,  
          GeneratedCell]; FrontEndTokenExecute["CellGroup"]; FrontEndTokenExecute["OpenCloseGroup"];

ABELL, M. L., BRASELTON, J. P., Mathematica by example, 2nd edition, Academic Press, 1997
CHUNG, T. J., Computational Fluid Dynamics, Cambridge University Press, 2002
INCROPERA, F. P., DeWITT, D. P., Fundamentals of Heat and Mass Transfer, 5th edition, Wiley, 2002
MALISKA, C. R., Heat transfer and Computational Fluid Mechanics (in portuguese), LTC, Brazil, 2004.
PATANKAR, S. V., Numerical heat transfer and Fluid Flow, Hemisphere, 1980.