1. Not finding help here? Sign up for a free 30min tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

CFD: 2-Dim Rayleigh-Benard problem

  1. Apr 1, 2013 #1
    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

    [tex]
    \frac{\partial u}{\partial x}+\frac{\partial v}{\partial x}=0
    [/tex]

    [tex]
    \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)
    [/tex]

    [tex]
    \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
    [/tex]

    [tex]
    \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)
    [/tex]

    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)
       % add Sidevectors
       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  
       % Add column Side vectors
       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  
       %Add column Side vectors
       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  
         % Add column Side vectors
         % 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.
     

    Attached Files:

    • avg.m
      avg.m
      File size:
      513 bytes
      Views:
      81
    • DD.m
      DD.m
      File size:
      1.3 KB
      Views:
      83
    • SG2212_2_03.m
      SG2212_2_03.m
      File size:
      10.1 KB
      Views:
      94
  2. jcsd
  3. Apr 1, 2013 #2
    Pdf-plots
     

    Attached Files:

  4. Aug 22, 2016 #3
    Hi,I guess in governing equation instead of Pr it should be (1/Pr)
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?
Draft saved Draft deleted



Similar Discussions: CFD: 2-Dim Rayleigh-Benard problem
  1. 2 Counting problems (Replies: 1)

Loading...