Problem Formulation, Solutions, Results Presentation, And Interpretation

It consists of the translation of the problem statement from plain words into mathematical language. The problem formulation defines the mathematical model, stemming from an interpretation of the real-world problem and from the setting of a number of assumptions. Typically, a formulation is given in terms of algebraic equations; of algebraic functions with unknown coefficients; or in terms of difference or differential equations and initial and/or boundary conditions expressed in terms of the primary variable.

• Explanations and Solutions
Solutions refer to the problem formulated. Thus, a correct solution may be found for either a correctly or an incorrectly formulated problem. Typically, solutions consist of coefficients’ values; closed-forms that satisfy difference or differential equations; or values of discretised fields for numerical problems. Solutions may also entail quantities derived from them (i.e. velocity at a given location derived from obtained position field; or heat flux at a given location once temperature field has been obtained). Detailed explanations and justifications of what has been done and why throughout the problem-solving procedure –from formulation to results presentation– will also be assessed.
• Results Presentation
This does not refer to the professional look of the report, which is assessed as Report Presentation (see above). Instead, results presentation refers to the relevance, clarity and completeness of the results presented (e.g. equations, data, tables, figures, etc.) to describe and help interpret the solution(s) obtained for the problem at hand.

• Results Interpretation (Discussion):
This refers to the student exhibiting clear understanding of the meaning of the results obtained and what they entail, instead of merely enumerating them. The problem-solving procedure does not end when finding a result but when the latter has been interpreted. When doing this, it is quite common to find that results do not make sense and therefore there must be some mistake somewhere. The importance of this component is evidenced by the high weight assigned to it in the marking rubric below. Different grades will be awarded according to the depth of this analysis.Do not confuse depth with length, and remember: “good things, when short, are twice as good”.

Initial and boundary conditions

Initial and boundary conditions

A.1 Initial condition: Initial conditions are required when one is solving a numerical problem. Initial conditions are defined at the start of a problem when the time is zero and the problem is marching type either in space or time.

If one is solving a problem of heat transfer in a three-dimensional cavity, initial condition for the temperature can be written as,

T_o = 0 at time t_o = 0 for all values of x, y and z.

Where T_o is the initial temperature and t_o is the initial time.

Similarly if one is solving a problem of fluid flow in a three-dimensional cavity, initial conditions can be written as,

U_o = 0 at time t_o = 0 for all values of x, y and z.

V_o = 0 at time t_o = 0 for all values of x, y and z.

W_o = 0 at time t_o = 0 for all values of x, y and z.

Where U_o, V_o and W_o are the initial velocities and t_o is the initial time.

Note: One can also take the values of T_o,_, U_o, V_o and W_o to be equal to a constant value, this value will be treated as a garbage value for initialization of the problem.

Boundary conditions: Boundary conditions are defined at the boundary walls of the physical domain of the problem.

Boundary conditions for temperature can be of different types like, Dirichlet boundary wall, adiabatic boundary wall, constant flux wall and convective wall.

Dirichlet boundary condition: This type of boundary condition is applied when the value of temperature is known. This type of boundary condition is also known as constant temperature boundary condition. Atmospheric temperature can be considered as a Dirichlet type boundary condition. Mathematically Dirichlet boundary condition can be expressed as,

T _wall = Constant

Adiabatic boundary condition: This type of boundary condition is applied when the boundary wall is insulated or no heat is getting transferred through that wall. This type of boundary condition is also known as insulated type boundary condition. Plastic on the electric wire is an example of insulated type boundary condition. Mathematically adiabatic boundary condition can be expressed as,

Where, T is temperature and n is the normal direction.

Boundary conditions for velocity can be like, inlet boundary condition, outlet boundary condition, wall boundary condition, symmetric & axisymmetric boundary condition, cyclic boundary condition and periodic boundary conditions.

 Crank - Nicholson (Finite Difference)

A.2 Crank - Nicholson is a type of finite difference scheme used to discretize the partial differential equation. This scheme is a type of implicit scheme and is difficult to solve and also requires more computational time per time step.

Crank-Nicholson (Finite Difference)

Disadvantage of this scheme is that, it is difficult to solve when compared with explicit scheme as more than one unknown appears while discretizing the partial differential equation. Implicit scheme is advantageous over explicit scheme, as it is free from the value of time step. This means a larger time step can be considered but increasing the value of time step also increases the truncation error due to the finite difference equation, so one has to optimize the value of accurate and appropriate solution.

Crank - Nicholson finite difference method is generally utilized for solving the heat conduction problem or diffusion problem. Formulation of the Crank-Nicholson method results in a tri-diagonal matrix problem.

1 - D transient problem: Crank - Nicholson method is unconditionally stable for 1 – D diffusion problem. For 1 – D problem Crank - Nicholson method has an accuracy of the order of [(?x) ² and (?t) ² ]

2 - D transient problem: Crank - Nicholson method is unconditionally stable for 2 – D diffusion problem. For 2 – D problem Crank - Nicholson method has an accuracy of the order of [(?x) ², (?y) ² and (?t) ² ]

Strong and weak form in finite elements

A.3 Both are used to solve the partial differential equation using finite element method. In strong form one solves the actual partial differential equation without any manipulation while in weak form one solves the partial differential equation in the integral form with the help of weigh functions.

Weak form does not mean that solution obtained will be inaccurate, weak form only help in reducing the complexity during the solution of differential equation by converting it into the integral form which is much easier to solve.

Damped harmonic oscillator with gravity

A.4 Equation of motion for the above problem can be written as,

The above equation is homogeneous and represents the free vibration. Above equation is linear differential partial differential equation and its solution can be written as,

Where, ‘e’ is base of natural algorithm, ‘t’ is time and s is a constant to calculate.

After differentiating the above equation and putting it back into the equation of motion will two values of ‘s’ as,

Given values:

m = 4 kg

k = 515+100 = 615 N / m

c = 30 kg / s

g = 9.81 m / s²

Above equation can be written as,

For x_o = 1 and v_o = 0, above equation can also be written as,

Strong and weak form in finite elements

Above equation can also be written The above equation is in the form of a homogeneous linear differential equation

The auxiliary equation can be written as,

f (m) = 0

m² + (k/m)=0

Putting the value of ‘k’ and ‘m’

k = 615 and m =4

m² + 153.75 = 0,

m₁ = 12.4i and m₂ = -12.4i

The roots are a pair of complex number, so the complementary function can be written as,

As the function is homogeneous its particular integral will be 0.

Now the complete solution can be written as

 Heated toroidal capacitor (conduction + convection)

A.5 Problem definition

As the toroidal capacitor is kept at the surface of the heat sink and heat is getting transferred to the heat sink and atmosphere by conduction and convection respectively. Conduction is assumed between the surfaces of heat sink and toroidal due to the fact that there is perfect contact at the interface as given in the problem while convection is assumed between the toroidal capacitor and ambient / atmosphere.

As the problem is symmetric in the radial direction we can consider the problem to be two dimensional with one boundary (heat sink) as constant temperature condition, two vertical boundaries as insulated boundary conditions and top wall as convective boundary condition.

Boundary conditions

Bottom wall: Constant temperature or Dirichlet type boundary condition

T = 40ºC

Top wall: Convective type boundary condition

Left wall: Adiabatic or insulated type boundary condition

Right wall: Adiabatic or insulated type boundary condition

After considering the above assumptions above problem can physically be represented

Above problem is a two-dimensional heat transfer problem with top & bottom boundaries at constant temperature and left & right boundaries as adiabatic or insulated wall. Above problem is also termed as Rayleigh-Benard convection problem.

Considering the problem to be independent of time we can write the governing equation of the problem asUsing the central differencing scheme we can write the above partial differential equation into algebraic form as,

Above equation can be written as,Geometry shown in the above figure can be converted into different grid points as shown below. In this geometry X and Y are the two-mutual perpendicular direction while ‘i’ and ‘j’ are the vector notations. ?x and ?y are the step sizes in horizontal and vertical direction respectively.

The boundary conditions in the vector form can be written as,

Bottom boundary (Heat sink): Constant Temperature at 40ºC,

T (i, 1) = 40; where I is varying from 0 to 38mm

Damped harmonic oscillator with gravity

Left vertical wall: (Insulated boundary); As we are at the left side of the physical domain we will have to use forward differencing scheme for calculation of the temperature. That means at the walls, temperature will be calculated with the help of temperature on the adjacent node in the positive x-direction which are inside the domain.

Above equation will give

Right vertical wall: (Insulated boundary); As we are at the right side of the physical domain we will have to use backward differencing scheme for calculation of the temperature. That means at the walls, temperature will be calculated with the help of temperature on the adjacent node in the negative x-direction which are inside the domain.

Above equation will give

Top wall: (Convective boundary); As we are at the top wall and from this wall heat is losing to the surrounding by convection-conduction. This can be written as,

As we are at the top wall of the physical domain we will have to use backward differencing scheme for calculation of the temperature. That means at the walls, temperature will be calculated with the help of temperature on the adjacent node in the negative y-direction which are inside the domain. One can write the above equation as,Above equation can be written in simple form as,

Result:

Maximum temperature inside the toroidal capacitor is found to be equal to 58ºC.

Heat flux transferred to the heat sink is 0.1721 W/m²

Heat flux transferred to the ambient is 86.415 W/m²

Figure below shows the temperature

Discussion:

From the figure it can be observed that temperature is increasing in the vertical direction from the bottom wall. It can also be observed that the distribution of the temperature is symmetric, this is due to the fact that left and right vertical walls are insulated and the problem is symmetric along the vertical axis.

MATLAB code:

%%Defining all the required variables

T_hs = 40; %Temperature of the heat sink

T_amb = 60; %Temperature of the ambient

Q_gen = ?515; %Internal heat generation inside the toroidal capacitor

k = 0.6;   %Thermal conductivity of the isotropic material

hc = 43;   %Heat transfer coefficient of convection

w = 0.038; %Width of the toroidal capacitor

nx = 39;   %Number of division in the horizontal direction  ( Width )

h = 0.085; %Height of the toroidal capacitor

ny = 86;   %Number of division in the vertical direction  ( Height )

dx = w / (nx-1 ); %Step size in the horizontal direction  ( Width )

dy = h / (ny-1 ); %Step size in the vertical direction  ( Height )

A = 2 / dx^2;

B = 2 / dy^2;

C = A+B;

Nu= k / (hc*dy);

%% Initial values of temperature

Temp (1: nx, 1: ny) = 40; %Assuming that the toroidal capacitor is at uniform temperature before the start of the problem

%%Boundary conditions

Temp (1: nx, 1) = 40; %Defining the temperature at the bottom wall (heat sink)

Temp (1, 2: ny-1) = Temp (2, 1: ny); %Defining the insulated boundary condition at the left wall

Temp (nx, 2: ny-1) = Temp (nx-1, 2: ny-1); %%Defining the insulated boundary condition at the right wall

Temp (1: nx, ny) = (T_amb+Nu*Temp (i, j-1))/ (1+Nu); %Defining the convective boundary condition at the top wall

%%Calculating the updated temperature

for i = 2: nx-1

    for j = 2: ny-1

        Temp_n (i, j ) =  ( ( (Temp (i+1, j )+Temp (i-1, j ) )/dx^2 )...

            + ((Temp (i, j+1) + Temp (i, j-1)) / dy^2)...

            + (Q_gen / k))/ C;

    end

end

Q.6 Transient heated rod

A.6 Problem definition

The physical domain of the above problem is rod and temperature is changing only in the horizontal direction. Considering the radius of the rod is very small one can neglect the variation of the temperature in the radial direction. This assumption helps in achieving, insulated enclosing boundaries in the radial direction.

So after assuming this, above problem can be considered as one-dimensional transient problem with constant heat generation inside the rod. Mathematically it can be represented as,

Mathematical formulation:

Above equation can be discretized with the help of FTCS (forward in time and central in space) method. Scheme used here for time discretization is explicit scheme while scheme used for spatial discretization is central difference scheme. Explicit scheme is adopted here to solve the problem as it is easy to formulate when compared with the implicit scheme. Central difference scheme is adopted for the discretization as it is one order accurate when compared with forward and backward discretization differencing schemes. After discretization above equation can be written as

Initial condition:

T = 0 ºC (x = 1 to 6 m)

 = 515 W / m³

k = 10 W / m - ºC

α = 1.11 × 10^-4 m² / s

A = 0.75 m²

Boundary conditions:

T _left = 1 ºC (x = 0)

T _right = 2 ºC (x = 8 m)

Physically the problem can be treated as shown in figure below,

Solutions and results

To analyse the results temperature distribution is plotted along the horizontal direction. Temperatures at the left and right boundaries of the rod are 1ºC and 2ºC respectively as shown in figure.

Above figure shows the variation of the temperature inside the rod. It can be seen from the figure that temperature is increasing from the left side to the centre of the rod, and then it is decreasing continuously. This increment and decrement in the temperature is due to the internal generation of the heat inside the rod. The maximum temperature developed inside the rod after reaching the time = 25000 seconds is around 4.2061 ºC.

From this figure it can be concluded that the internal heat generation inside the cavity will increases the temperature of the rod, this increment in temperature is due to the fact other boundaries of the rod (except left and right) are considered to be insulated which will not let pass the heat from its boundaries and will accumulate all the heat generated inside the rod.

MATLAB Code:

%Defining all the required variables

dx = 1; %Step size in horizontal direction

dt = 1; %Step size in time

nx = 9; %Number of grid points in horizontal direction

ny = 3; %Number of grid point in vertical direction

dy = 1; % Step size in vertical direction

%Initial condition for temperature inside the physical domain

Temp (1:9, 1:3 ) = 0;

q = 15; %Amount of heat flux generated inside the physical domain

k = 10; %Thermal conductivity of the material of the rod

alpha = 1.11*10^  (-4 ); %Thermal diffusivity of the material of the rod

Area = 0.75; %Area of the rod

%%Boundary condition for temperature on different walls

Temp  (1,1:3 ) = 1; %Left wall

Temp  (9,1:3 ) = 2; %Right wall

Temp_n  (1,1:3 ) = 1; %Left wall

Temp_n  (9,1:3 ) = 2; %Right wall

stability = alpha*dt /  (dx^2+dy^2 ); %Calculation of the stability condition 

Temp (2:8, 3) = Temp (2:8, 2 ); %Top wall   (insulated or adiabatic )

Temp (2:8, 1) = Temp (2:8, 2 ); %Bottom wall   (insulated or adiabatic )

%%Giving the end time up to which the code is to be run

t_final = 25000;

time = dt;

while time< = 25000 %Condition for time, after reaching to a value of 25000 code will stop automatically

    time

    %Calculating the temperature of the inside node points of the rod

    for i = 2:nx-1

        j = 2;

        Temp_n  (i, j ) = Temp  (i, j ) + alpha*dt*  (Temp  (i+1, j ) + Temp  (i-1, j )- 2 * Temp  (i, j ) ) *  (1/dx^2 )...

            +alpha * dt *  (Temp  (i, j+1 ) + Temp  (i, j-1 )- 2 *Temp  (i,j ) ) *  (1/dy^2 )...

            +alpha * dt * q / k;

    end

    %Updating the values of temperature by new values calculated on all the

    %node points

    for i = 2:nx-1

        Temp (i, 2 ) = Temp_n  (i,2 );

        Temp (i, 3 ) = Temp_n  (i,2 );

        Temp (i, 1 ) = Temp_n  (i,2 );

    end

    time = time+dt;.