MATLAB Examples

Natural convection in a tall cavity

This code solves the test problem of a thermally driven flow in a rectangular enclosure with an aspect ration of 8:1, as described in Christon et al. (2002). The method relies on Taylor-Hood finite elements on a regular graded triangular mesh for the spatial discretization and the Crank-Nicolson scheme with uniform step-size for the temporal discretization. The solver is documented in the PhD thesis of Sebastian Ullmann "POD-Galerkin Modeling for Incompressible Flows with Stochastic Boundary Conditions". Please reference the thesis or one of the related papers if you use the solver or parts of it for your own scientific work.

Copyright (C) Sebastian Ullmann 2015

Contents

Clean up storage

close('all');
clear('classes');

filename = '/worklocal/matlab.mat';  % use your favorite storage directory
if exist(filename,'file')
  delete(filename);
end
matObj = matfile(filename);

Create finite element solver object

To converge towards a periodic solution and not to a steady solution, a highly resolving mesh is needed, e.g. 'ndx', 32 and 'ndy', 128.

solverOptions = struct();
solverOptions.Pr = 0.71;   % Prandtl number
solverOptions.Ra = 3.4e5;  % Rayleigh number
solverOptions.ConditionsAtPoints = @conditionsAtPoints; % extra file
solverOptions.ConditionsAtEdges  = @conditionsAtEdges;  % extra file
solverOptions.Mesh ...
  = @()meshRectangularUnionJack('ndx',  32, ...
                                'ndy', 128, ...
                                'xmin', -0.5, ...
                                'xmax',  0.5, ...
                                'ymin', -4.0, ...
                                'ymax',  4.0, ...
                                'xgrading', 8, ...
                                'ygrading', 4);

fem = SolverUnsteadyBoussinesq(solverOptions);

Plot finite element mesh

figure();
trimesh(fem.mesh.triangles(:,1:3), ...
        fem.mesh.points(:,1), ...
        fem.mesh.points(:,2), ...
        zeros(size(fem.mesh.points,1),1), ...
        'EdgeColor','k');
view(2);
axis('equal');

Create preconditioner object

The preconditioner accelerates the iterative linear algebraic solver.

preconditioner = PreconditionerILU(fem);

Specify simulation output

The output object defines what data is saved to disk in each time step and specifies the graphical output. Save computation time by disabling the graphical output via 'display','none' instead of 'display','solution'.

output = Output(fem,'display','solution');

if ~strcmp(output.display,'none')
  % Open graphical window.
  figure('Renderer','zbuffer','Position',[100 100 640 480]);
end

Solve the problem

User options can be defined via keyword/value pairs. Have a look at the class file for more options and default values. What you see in the figure are spatial plots of the four components of the solution (pressure, x-velocity, y-velocity and temperature) as well as the time evolution at some point in the domain. The whole computation may take a few hours.

ticTime = tic();
fem.solve(matObj, ...
          'Time', 0:0.025:1200, ...
          'LinearSolver','bicgstab', ...
          'LinearTolerance',1e-2, ...
          'NonlinearTolerance',1e-12, ...
          'PreconditionerUpdate',@preconditioner.update, ...
          'PreconditionerApply',@preconditioner.apply, ...
          'OutputFunction',@output.compute);
toc(ticTime);

Compute the benchmark outputs

Be sure to have a fine enough spatial mesh and fine enough time step in order to obtain a fully periodic solution after about 1000 time units.

if matObj.t(1,end)>=1000
  out = output.postprocess(matObj);
end