Tutorial 21 – Unsteady Free Convection

Introduction

The aim of this tutorial is to show an unsteady flow simulation due to natural convection in a simple geometry.

Simulation setup

The shape of the domain is rectangular with a heated spot at the bottom (T=293K) and cold one on the top (T=273K). The remaining horizontal pieces have no BC prescribed, hence they are insulated. The vertical side walls have the symmetry boundary condition imposed. The mesh is of moderate size, notice the extra boundary layer generated at the top and bottom walls. We want to resolve the heat flux precisely. The whole fluid is initalized with zero velocity and the constant temperature equal 283K. The properties of fluid are somewhat artificial in order to speed-up the calculations. The flow will be driven by the buoyancy forces. In reality they arize due to small densiry variations in an almost-incompressible flow. Since the simulations are run with the incompressibility constraint, the so-called Bousinesq approximation is employed, where the body force results from the temperature difference,

Here beta is the thermal expansion coefficient, which is a materil property, and g is the gravity constant.
Of course the larger the the temperature difference, the stronger the buoyancy currents.

Post-processing

In order to visualize the simulations results we create a movie out of the simulation snapshots.
Again, to save time and memory, the setting for full-HD resolution are commented out. Feel free to test that at full resolution.

Watch the corculation currents develop.

clear;
clc;

% Import and display mesh
[p,e,t] = importMeshGmsh('twoPlates2.msh');

% Re-scale mesh to cm
p=p/50;

% Add boundary layer
layers.h1 = 0.0001;
layers.q = 1.5;
layers.n = 5;
wallIds = [11, 12,13];
[p,e,t] = extrudeLayers2D(p,e,t,wallIds,layers);

displayMesh2D(p,t);

% Convert mesh to the second order
[p,e,t,nVnodes,nPnodes,indices] = convertMeshToSecondOrder(p,e,t);
% Define fluid
nu = 1e-6; % kinematic viscosity [m^2/s]
lambda = 0.001; % thermal conductivity [W/(m*K)]
rho = 1.0; % density [kg/m^3]
cp = 100; % heat capacity [J/(kg*K)]
k = lambda/(rho*cp); %thermal diffusivity
Tref = 283; % temperature [K]
beta = 1/Tref; %thermal expansion coefficient [1/K]
g = -10; % gravity constant

% Init solution
[u, convergence] = initSolution(p,t,[0 0],0);
T = initScalarSolution(p,283);
T_low=273;
T_high=293;

T_range=[T_low, T_high];
U_range=[0, 0.1];

% Define timestepping
dt = 0.01;
nstep = 10;

% Assemble mass matrix
M = assembleMassMatrix2D(p,t,1/dt);
MT = assembleMassMatrix2D(p,t,1/dt);
% Shoot movie
movieObj1 = VideoWriter('unsteadyNaturalConvection_T.avi');
movieObj2 = VideoWriter('unsteadyNaturalConvection_Umag.avi');
% movieObj1.Quality=100;
% movieObj2.Quality=100;
x0=0;
y0=0;
% width=1920;
% height=1080;
open(movieObj1);
open(movieObj2);
for step = 1:nstep
time = step*dt;

% Assemble Navier-Stokes matrix
[NS, F] = assembleNavierStokesMatrix2D(p,e,t,nu,u(indices.indu),u(indices.indv),'nosupg');
% force field
fb = [zeros(1,nVnodes); -g*beta*(T'-Tref)];
F = [assembleVectorSourceTerm2D(p,t,fb); zeros(nPnodes,1)];
% Add unsteady terms
[NS,F] = addUnsteadyTerms(NS,F,M,u,indices);

% Impose boundary conditions
[NS, F] = imposeCfdBoundaryCondition2D(p,e,t,NS,F,[11 12 13],'wall');
[NS, F] = imposeCfdBoundaryCondition2D(p,e,t,NS,F,[14],'slipAlongY');

u = NS\F;

% Solve unsteady heat convection problem
FT = MT*T;

D = assembleDiffusionMatrix2D(p,t,k);
D = D + assembleScalarConvectionMatrix2D(p,t,k,u(indices.indu),u(indices.indv),'supgDoublyAsymptotic');
KT = MT + D;

% Apply thermal boundary conditions - temperatures on the walls
[KT,FT] = imposeScalarBoundaryCondition2D(p,e,KT,FT,12,'value',T_low);
[KT,FT] = imposeScalarBoundaryCondition2D(p,e,KT,FT,11,'value',T_high);

% Solve system of equations
T = KT\FT;

%display solutions and save movie snapshots
figure(2)
displaySolution2D(p,t,T,'Temperature',T_range);
drawnow()
%set(gcf,'units','points','position',[x0,y0,width,height])
currFrame1 = getframe(gcf);
writeVideo(movieObj1,currFrame1);
figure(3)
vmag=sqrt(u(indices.indu).^2+u(indices.indv).^2);
displaySolution2D(p,t,vmag,'Velocity magnitude',U_range);
% set(gcf,'units','points','position',[x0,y0,width,height])
currFrame2 = getframe(gcf);
writeVideo(movieObj2,currFrame2);
% Report time step number
disp(step);
end

close(movieObj1);
close(movieObj2);