================================================================================
 final_LSA - Linear Stability Analysis of Channel Flow over a Porous Layer
================================================================================
 
OVERVIEW
--------
This MATLAB code performs a temporal linear stability analysis (LSA) of a
pressure-driven flow in a channel that is split in two by a porous layer.
It does two things:
 
  1. Computes the laminar base-flow velocity profile U(y) across the full
     wall-to-wall domain, including flow through the porous layer, for a
     prescribed bulk Reynolds number.
 
  2. Solves a modified Orr-Sommerfeld eigenvalue problem for that profile to
     find the least stable (or most unstable) 2D disturbance modes at a
     given streamwise wavenumber, and reconstructs the velocity field of a
     chosen mode.
 
The porous layer is modelled with a Darcy-Forchheimer drag: a viscous
(Darcy) term set by a permeability tensor K and an inertial (Forchheimer)
term set by a tensor C. Both tensors are 2x2 (x-y components), so
anisotropic media can be represented. All quantities are dimensional (SI).
 
 
GEOMETRY
--------
  y = +(2*delta + t/2)  ================== top wall (no-slip)
                           fluid channel,  height 2*delta
  y = +t/2              ------------------
                           porous layer,   thickness t
  y = -t/2              ------------------
                           fluid channel,  height 2*delta
  y = -(2*delta + t/2)  ================== bottom wall (no-slip)
 
  delta = mat.delta     (half-height of each fluid channel)
  t     = mat.thickness (porous layer thickness)
 
The bottom and top channels can be driven by different pressure gradients
(mat.dpdx is a two-element vector: [bottom top]).
 
 
FILES
-----
modes.m
    Driver script. Sets the physical and numerical parameters, calls the
    base-flow solver and the eigenvalue solver, then post-processes one
    eigenmode. This is the file to run.
 
ctcf.m
    Base-flow solver.
      [u, du, ddu, y, dy, solid] = ctcf(mat)
    Solves the steady 1D momentum equation
      mu*U'' - mu*(K^-1)_22*U - rho*(C^-1)_22*|U|*U = dp/dx
    (drag terms active only inside the porous layer) on a uniform grid
    using a 4th-order central finite-difference stencil and U = 0 at both
    walls. The nonlinear Forchheimer term is handled by relaxed Picard
    iteration (inner loop). An outer loop rescales dp/dx until the bulk
    Reynolds number of the fluid regions matches mat.target_reb
    (relative tolerance 1e-4).
    Outputs:
      u, du, ddu - base velocity and its first/second derivatives
      y, dy      - grid (row vector) and uniform spacing
      solid      - logical mask, true inside the porous layer
 
modified_OS_FDM.m
    Eigenvalue solver.
      [c, v] = modified_OS_FDM(u, du, ddu, dy, solid, alp, mat)
    Discretises the Orr-Sommerfeld equation for the wall-normal velocity
    disturbance v, with extra Darcy and Forchheimer terms inside the
    porous layer, using 4th-order finite-difference stencils for the 1st,
    2nd and 4th derivatives. The two outermost grid points at each wall
    are removed (v = 0, which also approximates v' = 0 at the walls), and
    the resulting generalised eigenproblem  A*v = c*B*v  is solved with
    eigs for the 20 eigenvalues with the largest imaginary part.
    Outputs:
      c - eigenvalues (see "Interpreting the eigenvalues" below)
      v - eigenvectors, one per column, length N-4 (wall points removed)
    A commented-out 2nd-order version of the discretisation is kept in the
    file for reference.
 
 
REQUIREMENTS
------------
  - MATLAB R2018b or newer (uses max(...,'all') and eigs name-value
    options such as 'SubspaceDimension' and 'Tolerance').
  - No additional toolboxes.
 
 
HOW TO RUN
----------
  1. Put all three .m files in the same folder (or on the MATLAB path).
  2. Open MATLAB in that folder and run:
         >> modes
  3. Results are left in the workspace (see "Outputs of modes.m").
 
Run time is short for the default N = 500 grid points.
 
 
INPUT PARAMETERS (the "mat" struct, set in modes.m)
---------------------------------------------------
  Field            Default           Meaning
  ---------------  ----------------  -----------------------------------------
  mat.delta        0.01 m            Half-height of each fluid channel
  mat.mu           1.7893e-5 Pa s    Dynamic viscosity (air)
  mat.rho          1.225 kg/m^3      Density (air)
  mat.nu           mu/rho            Kinematic viscosity
  mat.N            500               Number of grid points, wall to wall
  mat.thickness    0.004 m           Porous layer thickness
  mat.K            1e-7 * eye(2)     Permeability tensor [m^2]
  mat.C            1e-3 * eye(2)     Forchheimer tensor (inertial drag term
                                     uses its inverse)
  mat.re_tau       [1e3 1e3]         Friction Reynolds numbers used only to
                                     build the initial dp/dx guess
  mat.dpdx         from re_tau       Initial pressure gradient [bottom top];
                                     rescaled by ctcf to hit target_reb
  mat.target_reb   4.5e3             Target bulk Reynolds number,
                                     Re_b = U_b*delta/nu, with U_b averaged
                                     over both fluid channels
 
  Other settings in modes.m:
  alp              2*pi/(6*delta)    Streamwise wavenumber [1/m]
                                     (wavelength = 6*delta)
  mode_number      1                 Which eigenmode to post-process
 
Setting K or C to a singular matrix (e.g. all zeros) switches off the
corresponding drag term.
 
 
INTERPRETING THE EIGENVALUES
----------------------------
Despite the variable name c, the eigenvalue returned by modified_OS_FDM is
the complex angular frequency omega = alpha * c_phase, for disturbances of
the form
 
    v'(x, y, t) = v(y) * exp( i*(alpha*x - omega*t) )
 
This is how modes.m uses it:
  - imag(c)        temporal growth rate [1/s]; > 0 means unstable
  - real(c)/alp    phase (convection) speed of the mode [m/s]
 
 
OUTPUTS OF modes.m (workspace variables)
----------------------------------------
  u, du, ddu, y        Base-flow profile and grid
  solid                Porous-layer mask
  c, v                 Eigenvalues and eigenvectors
  growth               Growth rate of the selected mode
  vp, up               Wall-normal and streamwise disturbance profiles of
                       the selected mode (up from continuity:
                       u' = -v'_y / (i*alpha))
  u_b                  Bulk velocity
  re_conv, re_ratio    Reynolds number based on the mode's phase speed, and
                       its ratio to the bulk Reynolds number
  t_final              Time at which the mode has grown ~10x (times 1.13)
  u_prime_local,       2D disturbance velocity fields over one wavelength
  v_prime_local        at t_final, normalised by the peak of v
 
modes.m does not currently produce any figures; the "Spatial plot"
section only prepares the fields above. They can be plotted with, e.g.:
    quiver(x_grid, y_grid, u_prime_local, v_prime_local)