# Homework Help: CFD: 2-Dim Rayleigh-Benard problem

1. Apr 1, 2013

### SebastianF

1. The problem statement, all variables and given/known data

Hi guys,

we are supposed to program the velocity and temperature fields for the 2-dim Navier-Stokes, using the pressure correction method and Boussinesq. approx for the temperature field. The relevant set of equations and B.C. is displayed below:

2. Relevant equations

$$\frac{\partial u}{\partial x}+\frac{\partial v}{\partial x}=0$$

$$\frac{\partial u}{\partial t}+\frac{\partial P}{\partial x}=-\frac{\partial u^2}{\partial x}-\frac{\partial uv}{\partial y}+Pr\Big( \frac{\partial^2 u}{\partial x^2} +\frac{\partial^2 u}{\partial y^2}\Big)$$

$$\frac{\partial v}{\partial t}+\frac{\partial P}{\partial y}=-\frac{\partial v^2}{\partial y}-\frac{\partial uv}{\partial x}+Pr\Big( \frac{\partial^2 v}{\partial x^2} +\frac{\partial^2 v}{\partial y^2}\Big)+f$$

$$\frac{\partial T}{\partial t}=-\frac{\partial vT}{\partial y}-\frac{\partial uT}{\partial x}+Pr\Big( \frac{\partial^2 T}{\partial x^2} +\frac{\partial^2 T}{\partial y^2}\Big)$$

where f = Ra * Pr * T
with Pr=Prandtl-number and Ra=Rayleigh Number

Boundary Conditions are:

x=0 : u=v=0 dT/dx=0
x=L_x: u=v=0 dT/dx=0
y=0 : u=v=0 T=T_bottom=0
y=L_y: u=v=0 T=T_top=0

Initial conditions are:

u=v=0
T=T_bottom+(y/L_y)*(T_top-T_bottom)

3. The attempt at a solution

I attempted to implement the above in Matlab:

The code runs but does not result in the expected vortices, but gives a strange velocity source at (1,1), I'm unsure where my error lies. Anyone have an idea?

Main Code:

Code (Text):
% Notes for the Rayleigh-Bernard problem
% 1. Approach
% Careful: Omega from the exercise sheet is a temperature difference, not
% the actual temperature
% Try to imitate the behaviour of V and see if that works for omega

lid_driven_cavity=0;

if (lid_driven_cavity==1)
% Parameters for test case I: Lid-driven cavity
% The Richardson number is zero, i.e. passive scalar.

Pr = 0.71;     % Prandtl number
Re = 10.;      % Reynolds number
Ri = 0.;       % Richardson number

dt = 0.00275;  % time step
Tf = 20;       % final time
Lx = 1;        % width of box
Ly = 1;        % height of box
Nx = 30;       % number of cells in x
Ny = 30;       % number of cells in y
ig = 200;      % number of iterations between output

% Boundary and initial conditions:
Utop = 1.;
% IF TEMPERATURE: Tbottom = 1.; Ttop = 0.;
% IF TEMPERATURE: namp = 0.;
else
% Parameters for test case II: Rayleigh-Bénard convection
% The DNS will be stable for Ra=1705, and unstable for Ra=1715
% (Theoretical limit for pure sinusoidal waves
% with L=2.01h: Ra=1708)
% Note the alternative scaling for convection problems.

Pr = 0.71;     % Prandtl number
Ra = 2715;     % Rayleigh number

Re = 1./Pr;    % Reynolds number
Ri = Ra*Pr;    % Richardson number

dt = 0.0005;   % time step
Tf = 20;       % final time
Lx = 10.;      % width of box
Ly = 1;        % height of box
Nx = 200;      % number of cells in x
Ny = 20;       % number of cells in y
ig = 200;      % number of iterations between output

% Boundary and initial conditions:
Utop    = 0.;
% IF TEMPERATURE: Tbottom = 1.; Ttop = 0.;
namp    = 0.1;
Tbottom = 1.;
Ttop    = 0.;
end

%------------------------------------------

% Number of iterations
Nit = ceil(Tf/dt);
% Spatial grid: Location of corners -> position of the vel. nodes
x =  linspace(0,Lx,Nx+1);
y =  linspace(0,Ly,Ny+1);
% Grid spacing
hx = Lx/Nx;
hy = Ly/Ny;
% Boundary conditions:
uN = x*0+Utop;    vN = avg(x,2)*0;
uS = x*0;         vS = avg(x,2)*0;
uW = avg(y,2)*0;  vW = y*0;
uE = avg(y,2)*0;  vE = y*0;
% all row vectors for now
% actual computation of Boundary condtions at Ue, Ve

% Initial conditions
% Remember that i=1...x values end up in one colum U(i,j)=U(row,colum) and
% all           j=1...y values        in one row   V(i,j)=V(row,colum)
U = zeros(Nx-1,Ny);
V = zeros(Nx,Ny-1);
% linear profile for T with random noise
% IF TEMPERATURE: T = ... + namp*rand(Nx,Ny)
T     = ones(Nx,Ny)* Tbottom + (Ttop - Tbottom)*repmat(avg(y,2),Nx,1) + namp*rand(Nx,Ny);
% Omega = diff(T,1,1);

% T needs wider Boundary conditions vector: Nx, Ny
% Note: Although the spacing is incorrect this is ultimately irrelavant. xT
% will be used to save (B.C.) temperatures, not positions!
xT =  linspace(0,Lx,Nx);
yT =  linspace(0,Ly,Ny);
% temperature boundary conditions:
tN = xT*0+Ttop;   tS = xT*0+Tbottom;
% tW = T(1,:);      tE = T(end,:);

% Time series
tser = [];
Tser = [];

%-----------------------------------------

% Compute system matrices for pressure
% First set homogeneous Neumann condition all around
% Laplace operator on cell centres: Fxx + Fyy
Lp = kron(speye(Ny), DD(Nx,hx) ) + kron( DD(Ny,hy) ,speye(Nx));
% Set one Dirichlet value to fix pressure in that point
Lp(1,:) = 0 ; Lp(1,1) = 1 ;
% rhs = 0 inside the loop; further down
% For more speed, you could pre-compute the LU decomposition
% [LLp,ULp] = lu(Lp);
%-----------------------------------------

% Progress bar
fprintf(...
'[         |         |         |         |         ]\n')

%-----------------------------------------

% Main loop over iterations

for k = 1:Nit

% include all boundary points for u and v (linear extrapolation
% for ghost cells) into extended array (Ue,Ve)
Ue     = [uW; U; uE;];
% rightmost column (translates to bottommost row in actual matrix)
% choosen, so that the average value implemented in Ua will yield uN
Ue     = [(2*uS')-Ue(:,1) Ue (2*uN')-Ue(:,end)];
Ve     = [vS' V  vN'];
Ve     = [(2*vW)-Ve(1,:); Ve; (2*vE)-Ve(end,:);];

% boundary conditions for the Temperature
Te     = [(2*tS')-T(:,1) T (2*tN')-T(:,end)];
% Add row top & bottom vectors.
Te     = [Te(1,:);  Te;  Te(end,:);];

% averaged (Ua,Va) of u and v on corners
% Dimensionaltity of Ua,Va should be (nx+1)*(ny+1) for both
% (nx+1)*(ny+2) => (nx+1)*(ny+1)
% averaging in y-direction (Up/down)
Ua = avg(Ue,2);
% (nx+2)*(ny+1) => (nx+1)*(ny+1)
% averaging in x-direction (left/right)
Va = avg(Ve,1);
% averaging projects on cell corners

% Example RayleighBenardConvection used average value, so stay parallel
% and see what happens; Note that alone averaging will return:
% (nx+2)*(ny+2) => (nx+1)*(ny+2) => projection on V
Ta = avg(Te,1);
% (nx+1)*(ny+2) => (nx+1)*(ny+1) => projection on U
Ta = avg(Ta,2);
% Note: taking the difference with the lower boundary value
% (y-direction) should happen at some point -> Omega!

% construct individual parts of nonlinear terms
dUVdx = diff(Ua.*Va, 1, 1)/hx;
dUVdy = diff(Ua.*Va, 1, 2)/hy;
% diff in x-direction (up/down) in T-Matrix
% leaving v-distro
dUOdx = diff(Ua.*Ta, 1, 1)/hx;
% diff in y-direction (left/right) in T-Matrix
% leaving u-distro
dVOdy = diff(Va.*Ta, 1, 2)/hy;

dU2dx = diff(avg(Ue(:,2:end-1).^2,1), 1, 1)/hx;
dV2dy = diff(avg(Ve(2:end-1,:).^2,2), 1, 2)/hy;

% treat viscosity explicitly
viscu = diff( Ue(:,2:end-1),2,1 )/hx^2 + diff( Ue(2:end-1,:),2,2 )/hy^2;
viscv = diff( Ve(:,2:end-1),2,1 )/hx^2 + diff( Ve(2:end-1,:),2,2 )/hy^2;
%       (nx+2)*(ny+2)=>(nx)*(ny)       + (nx+2)*(ny+2)=>(nx)*(ny)
viscT = diff( Te(:,2:end-1),2,1 )/hx^2 + diff( Te(2:end-1,:),2,2 )/hy^2;

% buoyancy term
% IF TEMPERATURE: fy = ...
% (nx)*(ny)
% dUOdx(:,2:end-1) or dUOdx(:,1:end-1) or dUOdx(:,2:end)
T  = T + dt*viscT - dt*(dUOdx(:,1:end-1)+dVOdy(1:end-1,:));
% T  = T + viscT - (dUOdx(:,2:end)+dVOdy(2:end,:));
% driving force is the temperature DIFFERENCE:
fy = Ri*diff(T,1,2);

% compose final nonlinear term + explicit viscous terms
U = U + dt/Re*viscu - dt*( dU2dx + dUVdy(2:end-1,:));
% V = V + dt/Re*viscv - dt*( dV2dy + dUVdx(:,2:end-1)); % + % IF TEMPERATURE: dt*fy;
V = V + dt/Re*viscv - dt*( dV2dy + dUVdx(:,2:end-1)) + dt*fy;

% pressure correction, Dirichlet P=0 at (1,1)
rhs = (diff([uW; U; uE],1,1)/hx + diff([vS' V  vN'],1,2)/hy)/dt;
rhs = reshape(rhs,Nx*Ny,1);
rhs(1) = 0;
P = Lp\rhs;
% alternatively, you can use the pre-computed LU decompositon
% P = ULp\(LLp\rhs);
P = reshape(P,Nx,Ny);

% apply pressure correction
U = U - dt*diff(P,1,1)/hx;
V = V - dt*diff(P,1,2)/hy;

% Temperature equation: Recomputation with pressure corrected velocities!
% is enforcement of B.C. necessary at this step?
% boundary conditions for the Temperature
Te = [(2*tS')-T(:,1) T (2*tN')-T(:,end)];
% Add row top & bottom vectors.
Te = [Te(1,:);  Te;  Te(end,:);];
% IF TEMPERATURE: Tu = ...
Ta = avg(Te,1);
Ta = avg(Ta,2);
dUOdx = diff(Ua.*Ta, 1, 1)/hx;
% IF TEMPERATURE: Tv = ...
dVOdy = diff(Va.*Ta, 1, 2)/hy;
% IF TEMPERATURE: H = ...
viscT = diff( Te(:,2:end-1),2,1 )/hx^2 + diff( Te(2:end-1,:),2,2 )/hy^2;
% IF TEMPERATURE: T = T + dt*H;
T  = T + dt*viscT - dt*(dUOdx(:,2:end)+dVOdy(2:end,:));

%-----------------------------------------

% progress bar
if floor(51*k/Nit)>floor(51*(k-1)/Nit), fprintf('.'), end

% plot solution if needed
if k==1|floor(k/ig)==k/ig

% compute divergence on cell centres
if (1==1)
div = diff([uW;U;uE])/hx + diff([vS' V vN'],1,2)/hy;

figure(1);clf; hold on;
contourf(avg(x,2),avg(y,2),div');colorbar
axis equal; axis([0 Lx 0 Ly]);
title(sprintf('divergence at t=%g',k*dt))
drawnow
end

% compute velocity on cell corners
% renew Ua, Va matrixes with freshly computed values of U and V
Ue = [uW; U; uE;];
Ue = [(2*uS')-Ue(:,1) Ue (2*uN')-Ue(:,end)];
Ve = [vS' V  vN'];
Ve = [(2*vW)-Ve(1,:); Ve; (2*vE)-Ve(end,:);];
Ua = avg(Ue,2);
Va = avg(Ve,1);
% Ua = ...
% Va = ...
Len = sqrt(Ua.^2+Va.^2+eps);

figure(2);clf;hold on;
%contourf(avg(x,2),avg(y,2),P');colorbar
contourf(x,y,sqrt(Ua.^2+Va.^2)',20,'k-');colorbar
quiver(x,y,(Ua./Len)',(Va./Len)',.4,'k-')
axis equal; axis([0 Lx 0 Ly]);
title(sprintf('u at t=%g',k*dt))
drawnow

% IF TEMPERATURE: % compute temperature on cell corners
% IF TEMPERATURE: Ta = ...
% boundary conditions for the Temperature
% Te = [(2*tS')-T(:,1) T (2*tN')-T(:,end)];
% Add row top & bottom vectors.
% Te = [Te(1,:);  Te;  Te(end,:);];
Ta = avg(Te,1);
Ta = avg(Ta,2);

figure(3); clf; hold on;
contourf(x,y,Ta',20,'k-');colorbar
quiver(x,y,(Ua./Len)',(Va./Len)',.4,'k-')
axis equal; axis([0 Lx 0 Ly]);
title(sprintf('T at t=%g',k*dt))
drawnow

% Time history
if (1==1)
figure(4); hold on;
tser = [tser k*dt];
Tser = [Tser Ue(ceil((Nx+1)/2),ceil((Ny+1)/2))];
plot(tser,abs(Tser))
title(sprintf('Probe signal at x=%g, y=%g',...
x(ceil((Nx+1)/2)),y(ceil((Ny+1)/2))))
set(gca,'yscale','log')
xlabel('time t');ylabel('u(t)')
end
end
end
fprintf('\n')
Function: avg
Code (Text):
function B = avg(A,idim)
% AVG(A,idim)
%
% Averaging function to go from cell centres (pressure nodes)
% to cell corners (velocity nodes) and vice versa.
% avg acts on index idim; default is idim=1.
%
% This function belongs to SG2212.m

% nargin returns number of function arguments -> idim not entered
if nargin<2, idim = 1; end

if (idim==1)
B = [ A(2:end,:) + A(1:end-1,:) ]/2;
elseif (idim==2)
B = [ A(:,2:end) + A(:,1:end-1) ]/2;
else
error('avg(A,idim): idim must be 1 or 2')
end
Function: DD
Code (Text):
function A = DD(n,h)
% DD(n,h)
%
% One-dimensional finite-difference derivative matrix
% of size n times n for second derivative:
% h^2 * f''(x_j) = -f(x_j-1) + 2*f(x_j) - f(x_j+1)
%
% Homogeneous Neumann boundary conditions on the boundaries
% are imposed, i.e.
% f(x_0) = f(x_1)
% if the wall lies between x_0 and x_1. This gives then
% h^2 * f''(x_j) = + f(x_0) - 2*f(x_1) + f(x_2)
%                = + f(x_1) - 2*f(x_1) + f(x_2)
%                =              f(x_1) + f(x_2)
%
% For n=5 and h=1 the following result is obtained:
%
% A =
%
%    -1     1     0     0     0
%     1    -2     1     0     0
%     0     1    -2     1     0
%     0     0     1    -2     1
%     0     0     0     1    -1
%
% This function belongs to SG2212.m

% produce vector, filled with ones, for further use in A:
% Note that while diag will contain n numbers and therefore will
% be larger than required for the n-1 length sidediagonals, spdiags
% will simply cut-off the superflous information.
diag = ones(n,1);

A = spdiags([diag -2*diag diag], -1:1, n, n)/h^2;
A(1,1) = -1/h^2;
A(n,n) = -1/h^2;
Attachment contains the matlab m-files of the code fragments above and PDF-pictures of the Temp and velocity fields.

File size:
513 bytes
Views:
85
File size:
1.3 KB
Views:
86
File size:
10.1 KB
Views:
97
2. Apr 1, 2013

### SebastianF

Pdf-plots

#### Attached Files:

File size:
248.3 KB
Views:
140
• ###### VelField.pdf
File size:
249.5 KB
Views:
123
3. Aug 22, 2016

### Md Haidar

Hi,I guess in governing equation instead of Pr it should be (1/Pr)