Tutorial 8 – Unsteady Heat Transfer

Posted on July 17, 2017 - Tutorials

This tutorial is intended for the full version of the toolbox.

Introduction

This tutorial is continuation of Tutorial 7 on Steady State Heat Conduction in a Two-Component Material and it is adviced that you first analyse in details the previous tutorial, since we assume here that certain concepts and portions of code will be clear for you.

Problem Definition

We will simulate transient heat conduction in a two-component material composed of two distinct phases which differ in heat conductivity. Geometry with ids of Physical Lines and ids of Physical Surfaces are shown in figure below.

The problem is described by the following mathematical model

    $$\rho c_p \frac{\partial T}{\partial t} = \nabla \cdot (\lambda \nabla T)$$

We assume the following boundary conditions:

  • Line 12: Constant heat flux of 10 000 W/m^2
  • Line 13: Insulated
  • Line 14: Heat convection with htc = 20 W/(m2K) and ambient temperature of 293 K.

The initial condition is a temperature of 293 K everywhere.

Code

In this section we are going to show the Matlab script which solves the problem in QuickerSim CFD Toolbox. We will mainly describe the code in comments, since majority of steps should be clear after completing Tutorial 7.

clc;
clear;


% Import and display mesh
[p,e,t] = importMeshGmsh('simpleRib.msh');
displayMesh2D(p,t);
<img class="size-full wp-image-2396 aligncenter" src="http://quickersim.com/cfdtoolbox/wp-content/uploads/2017/07/QuickerSim-Heat-Transfer-Rib-mesh.png" alt="" width="560" height="420" />


% Define timestepping parameters
nstep = 120; % number of steps
dt = 240; % time step size


% Define material parameters
rho = 2000; % density
cp = 1200; % heat capacity
lambda1 = 50; % thermal conductivity of the surrounding material
lambda2 = 400; % thermal conductivity of the rib material


% Assign thermal conductivities to appropriate subdomains
nelements = size(t,2);
lambda = lambda2*ones(nelements,1);
lambda(t(end,:)==19) = lambda1;


% Init solution with constant value
[T,convergence] = initScalarSolution(p,293);


% Precompute mass matrix and diffusion matrix
M = assembleMassMatrix2D(p,t,rho*cp/dt);
D = assembleDiffusionMatrix2D(p,t,lambda);

We will shoot a movie of the temperature field evolution. To accompish that, we will use $\tt{initMovie}$ function from our toolbox. Additionally, we exploit Euler method for timestepping. This resuts in following iterative process for timestepping (k in upper index denotes time step index):

    $$\rho c_p \frac{T^{k+1} - T^{k}}{\Delta t} = \nabla \cdot (\lambda \nabla T^{k+1})$$

The discretized for with mass matrix M and diffusion matrix D reads:

    $$M \cdot (T^{k+1}-T^k) = -D\cdot T^{k+1}$$

which can be rearranged as

    $$(M + D) \cdot T^{k+1} = M \cdot T^{k}$$

% We realize this algorithm in the timestepping loop below
% Init movie
movieObj = initMovie('ribHeatCondution.avi',16);


% Begin timestepping
for step = 1:nstep
    time = step*dt;

    % Compute transient part of the right hand side
    F = M*T;

    % Assemble problem matrix
    K = M + D;

    % Impose boundary conditions
    [K, F] = imposeScalarBoundaryCondition2D(p,e,K,F,12,'flux',10000);
    [K, F] = imposeScalarBoundaryCondition2D(p,e,K,F,14,'robin',20*293,20);

    % Solve for temperature field at the new time level
    T = K\F;

    % Report step number (uncomment code line below)
    % disp(step);

    % Display current temperature field and shoot movie frame
    figure(1);
    displaySolution2D(p,t,T,'Temperature [deg C]');
    drawnow;
    writeVideo(movieObj,getframe(gcf));
end
<img class="size-full wp-image-2397 aligncenter" src="http://quickersim.com/cfdtoolbox/wp-content/uploads/2017/07/QuickerSim-Heat-Transfer-Rib-Temperature-Map.png" alt="" width="560" height="420" />



% Close movie file
close(movieObj);


 

Tags