Natural Convection Simulation in QuickerSim CFD Toolbox for MATLAB

We are going to present in this code example, how a natural convection problem can be solved in QuickerSim CFD Toolbox for MATLAB. Our benchmark problem is defined by unstationary Navier-Stokes equations for the fluid flow, transient convection-diffusion equation for convective heat transfer and the momentum source term is supplemented with additional force due to buyoancy modelled by Boussinesq approximation.

The movie below shows the result of simulation:

The problem can be solved in the full version of QuickerSim CFD Toolbox (also student licenses are accessible) with the code given below:

% QuickerSim CFD Toolbox: Solve unsteady natural convection problem in 2-D

clear;
clc;

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

% Convert mesh to second order
[p,e,t,nVnodes,nPnodes,indices] = convertMeshToSecondOrder(p,e,t);

% Define fluid
nu = 0.005; % kinematic viscosity
lambda = 0.6; % thermal conductivity
rho = 1.225;
cp = 100;
k = lambda/(rho*cp);
Tref = 283;
beta = 1/Tref;
g = -9.81;

% Init solution or solve a stationary problem
[u, convergence] = initSolution(p,t,[0 0],0);
T = initScalarSolution(p,283);

% Define timestepping settings
dt = 0.5;
nstep = 40;

% Assemble mass matrix for flow problem and temperature problem
M = assembleMassMatrix2D(p,t,1/dt);
MT = M;

% Shoot movie
movieObj = initMovie('unsteadyNaturalConvection.avi',12);

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');
    
    % Introduce buyoancy as momentum source term according to Boussinesq
    % approximation
    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,[7 8 9 10],'wall');
    
    u = NS\F;
    
    % Solve unsteady heat convection problem (left wall with T = 0, right
    % with T = 20, top and bottom walls insulated)
    FT = MT*T;
    
    D = assembleDiffusionMatrix2D(p,t,k);
    D = D + assembleScalarConvectionMatrix2D(p,t,k,u(indices.indu),u(indices.indv),'nosupg');
    KT = MT + D;

    % Apply boundary conditions
    [KT,FT] = imposeScalarBoundaryCondition2D(p,e,KT,FT,10,'value',273);
    [KT,FT] = imposeScalarBoundaryCondition2D(p,e,KT,FT,8,'value',293);

    % Solve system of equations
    T = KT\FT;
    
    % Preview flow field
    figure(3)
    subplot(1,2,1)
    displaySolution2D(p,t,u(indices.indv),'y-velocity');
    subplot(1,2,2)
    displaySolution2D(p,t,T-273,'Temperature [deg C]');
    drawnow;
    writeVideo(movieObj,getframe(gcf));
    
    % Report time step number
    disp(step);
end

close(movieObj)

Tags

Search