Three-dimensional channel flow

Three-dimensional channel flow geometry.

The linearized evolution model governing the dynamics of three-dimensional velocity fluctuations in channel flows is given by

$begin{array}{rcl} Delta , phi_{1, t} (k_x, k_z, y, t) & !! = !! & left( {displaystyle frac{1}{R}} Delta^{2} , - , mathrm{i} , k_x , U(y) , Delta , + , mathrm{i} , k_x , U''(y) right) phi_{1}(k_x, k_z, y, t) [0.3cm] && , - , mathrm{i} , k_x , D^{(1)} , d_{1}(k_x, k_z, y, t) , - , kappa^{2} , d_{2}(k_x, k_z, y, t) , - , mathrm{i} , k_z , D^{(1)} , d_{3}(k_x, k_z, y, t) [0.3cm] phi_{2,t} (k_x, k_z, y, t) & !! = !! & left( {displaystyle frac{1}{R}} Delta , - , mathrm{i} , k_x , U(y) right) phi_{2}(k_x, k_z, y, t) , - , mathrm{i} , k_z , U'(y) , phi_{1}(k_x, k_z, y, t) [0.3cm] && , + , mathrm{i} , k_z , d_{1}(k_x, k_z, y, t) , - , mathrm{i} , k_x , d_{3}(k_x, k_z, y, t), ;;;; y , in , left[ -1, , 1 right] end{array}$

where

,	—	wall-normal velocity and vorticity
, ,	—	streamwise, wall-normal, and spanwise forcing
	—	Reynolds number
,	—	streamwise and spanwise wavenumbers
	—	for shear-driven flow
	—	for pressure-driven flow
$Delta = D^{(2)} - kappa^2$	—	Laplacian

$Delta^2 = D^{(4)} - 2 , kappa^2 , D^{(2)} + kappa^4$ .

The boundary conditions are given by

$begin{array}{rcl} phi_{1}(k_x, pm 1, t) & !! = !! & D^{(1)} phi_{1}(k_x, pm 1, t) ; = ; 0 [0.15cm] phi_{2}(k_x, pm 1, t) & !! = !! & 0. end{array}$

The desired outputs are the streamwise, wall-normal, and spanwise velocity fluctuations,

$begin{array}{rcl} u (k_x, k_z, y, t) & !! = !! & {displaystyle frac{1}{kappa^{2}}} , left( mathrm{i} , k_x , D^{(1)} , phi_{1}(k_x, k_z, y, t) , - , mathrm{i} , k_z , phi_{2}(k_x, k_z, y, t) right) [0.25cm] v (k_x, k_z, y, t) & !! = !! & phi_{1}(k_x, k_z, y, t) [0.25cm] w (k_x, k_z, y, t) & !! = !! & {displaystyle frac{1}{kappa^{2}}} left( mathrm{i} , k_z , D^{(1)} , phi_{1}(k_x, k_z, y, t) , + , mathrm{i} , k_x , phi_{2}(k_x, k_z, y, t) right). end{array}$

The input-output differential equations representing the frequency response operator are given by

$begin{array}{rcl} {displaystyle left( frac{1}{R} , D^{(4)} , + , a_{11,2}(y) , D^{(2)} , + , a_{11,0}(y) right) , phi_1 (y) } & !! = !! & mathrm{i} , k_x , D^{(1)} , d_{1}(y) , + , kappa^2 , d_{2}(y) [0.25cm] && , + , mathrm{i} , k_z , D^{(1)} , d_{3}(y) [0.35cm] {displaystyle left( frac{1}{R} , D^{(2)} , - , a_{22,0} (y) right) phi_2 (y) , - , a_{21,0}(y) , phi_1 (y) } & !! = !! & - mathrm{i} , k_z , d_1 (y) , + , mathrm{i} , k_x , d_{3} (y) end{array}$

where

$begin{array}{rcl} a_{11,2}(y) & !! = !! & {displaystyle -left( frac{2 , kappa^2}{R} , + , mathrm{i} , k_x , U(y) , + , mathrm{i} , omega right) } [0.45cm] a_{11,0}(y) & !! = !! & {displaystyle frac{kappa^4}{R} , + , mathrm{i} , k_x , kappa^{2} , U(y) , + , mathrm{i} , k_x , U''(y) , + , mathrm{i} , omega , kappa^2 } [0.45cm] a_{21,0}(y) & !! = !! & - , mathrm{i} , k_x , U'(y). end{array}$

Matlab codes

Determine the most amplified flow structures for the 3D linearized Navier-Stokes equations in a pressure-driven channel flow with R = 2000 , k_x = k_z = 1 , and omega = -0.385 .

% System parameters:
R = 2000;  % Reynolds number

kxval = [1 1 -1 -1];        % streamwise wave-number
kzval = [1 -1 1 -1];        % spanwise wave-number
omval  = 0.385*[-1 -1 1 1]; % temporal frequency

dom = domain(-1,1);     % domain of your function
y = chebfun('y',dom);
fone = chebfun(1,dom);      % fone(y) = 1
fzero = chebfun(0,dom);     % fzero(y) = 0

U = diag(1 - y.^2); % in pressure-driven flow
Uy = diag(-2*y);
Uyy = diag(-2*fone);

N = 100;    % number of collocation points for plotting
yd = chebpts(N);

% Boundary conditions
Wa0{1} = [1, 0, 0, 0; 0, 1, 0, 0];
Wa0{2} = [1, 0];
Wb0{1} = [1, 0, 0, 0; 0, 1, 0, 0];
Wb0{2} = [1, 0];

% Looping over kx, kz, and om
uvec = zeros(N,length(kxval)); wxvec = zeros(N,length(kxval));

for n = 1:length(kxval),

    kx = kxval(n); kz = kzval(n); om = omval(n);
    kx2 = kx*kx; kz2 = kz*kz;
    k2 = kx2 + kz2; k4 = k2*k2;

    % coefficients of the operator A0
    % a110*phi1 + 0*phi1' + a112*phi1'' + 0*phi1''' + (1/R)*phi1''''
    a112 = -(1/R)*2*k2*fone - 1i*kx*U*fone - 1i*om*fone;
    a110 = (1/R)*k4*fone + 1i*kx*k2*U*fone + 1i*kx*Uyy*fone ...
        + 1i*om*k2*fone;

    % a220*phi2 + 0*phi2' + (1/R)*phi2'' + a210*phi1
    a220 = -(1/R)*k2*fone - 1i*kx*U*fone - 1i*om*fone;
    a210 = -1i*kz*Uy*fone;

    A0{1,1} = [a110, fzero, a112, fzero, (1/R)*fone];
    A0{1,2} = [fzero, fzero, fzero];
    A0{2,1} = [a210, fzero, fzero, fzero, fzero];
    A0{2,2} = [a220, fzero, (1/R)*fone];

    % coefficients of the operator B0
    % 0*d1 + i*kx*d1' + k2*d2 + 0*d2' + 0*d3 + i*kz*d3'
    B11 = [fzero, 1i*kx*fone]; B12 = [k2*fone, fzero];
    B13 = [fzero, 1i*kz*fone];

    % -i*kz*d1 + 0*d1' + 0*d2 + 0*d2' + i*kx*d3 + 0*d3'
    B21 = [-1i*kz*fone, fzero]; B22 = [fzero, fzero];
    B23 = [1i*kx*fone, fzero];

    B0 = {B11, B12, B13; B21, B22, B23};

    % coefficients of the operator B0
    % u = 0*phi1 + (1/k2)*i*kx*phi1' - (1/k2)*i*kz*phi2 + 0*phi2'
    C11 = [fzero,(1/k2)*1i*kx*fone]; C12 = [-(1/k2)*1i*kz*fone, fzero];
    % v = 1*phi1 + 0*phi1' + 0*phi2 + 0*phi2'
    C21 = [fone, fzero]; C22 = [fzero, fzero];
    % v = 0*phi1 + (1/k2)*i*kz*phi1' + (1/k2)*i*kx*phi2 + 0*phi2'
    C31 = [fzero,(1/k2)*1i*kz*fone]; C32 = [(1/k2)*1i*kx*fone, fzero];

    C0 = {C11, C12; C21, C22; C31, C32};

    % solving for the left principle singular pair
    [Sfun, Sval] = svdfr(A0,B0,C0,Wa0,Wb0,1,1);

    % velocities
    ui = Sfun{1}; % streamwise velocity
    vi = Sfun{2}; % wall-normal velocity
    wi = Sfun{3}; % spanwise velocity

    % streamwise vorticity, wx = w' - i*kz*v
    wx = diff(wi(:,1)) - 1i*kz*vi(:,1);
    % discretized values for plotting
    uvec(:,n) = ui(yd,1); wxvec(:,n) = wx(yd,1);

end

Determine the streamwise velocity and vorticity in the physical space.

kx = abs(kxval(1)); kz = abs(kzval(1));
zval = linspace(-7.8, 7.8, 100);    % spanwise coordinate
xval = linspace(0, 12.7, 100);      % streamwise coordinate

Up = zeros(length(zval),length(xval),N);    % physical value of u
Wxp = zeros(length(zval),length(xval),N);   % physical value of wx

for indx = 1:length(xval)

    x = xval(indx);

    for indz = 1:length(zval)

        z = zval(indz);

        for n = 1:4

            kx = kxval(n); kz = kzval(n);

            Up(indz,indx,:) =  squeeze(Up(indz,indx,:)) + ...
                uvec(:,n)*exp(1i*kx*x + 1i*kz*z);

            Wxp(indz,indx,:) =  squeeze(Wxp(indz,indx,:)) + ...
                wxvec(:,n)*exp(1i*kx*x + 1i*kz*z);

        end
    end
end

Up = real(Up); Wxp = real(Wxp); % only real part exist

Plot most amplified structures of streamwise velocity fluctuations (3D isosurface plots).

p = patch(isosurface(xval,zval,yd,Up,0.65));
isonormals(xval,zval,yd,Up,p)
set(p,'FaceColor','red','EdgeColor','none');
daspect('auto')
view(3);
axis([xval(1) xval(end) zval(1) zval(end) yd(1) yd(end)]);
camlight
lighting phong

hold on
p = patch(isosurface(xval,zval,yd,Up,-0.65));
isonormals(xval,zval,yd,Up,p)
set(p,'FaceColor','green','EdgeColor','none');
daspect('auto')
view(3);
axis([xval(1) xval(end) zval(1) zval(end) yd(1) yd(end)]);
camlight
lighting phong

xlab = xlabel('x', 'interpreter', 'tex');
ylab = ylabel('z', 'interpreter', 'tex');
zlab = zlabel('y', 'interpreter', 'tex');
set(xlab, 'FontName', 'cmmi10', 'FontSize', 25);
set(ylab, 'FontName', 'cmmi10', 'FontSize', 25);
set(zlab, 'FontName', 'cmmi10', 'FontSize', 25);
h = get(gcf,'CurrentAxes');
set(h,'FontName','cmr10','FontSize',15,'xscale','lin','yscale','lin');
hold off

Plot most amplified structures of streamwise velocity (2D color plots) and streamwise vorticity (contour plots) fluctuations.

xslice = 1;
pcolor(zval,yd,squeeze(Up(:,xslice,:))');
shading interp;

% normalizing the streamwise vorticity for plotting
Wxpn = Wxp/max(max(max(Wxp)));

hold on
contour(zval,yd,squeeze(Wxpn(:,xslice,:))',...
    linspace(0.1*max(max(max(Wxpn))),0.9*max(max(max(Wxpn))),4),'k--','LineWidth',1.1)
contour(zval,yd,squeeze(Wxpn(:,xslice,:))',...
    linspace(-0.9*max(max(max(Wxpn))),-0.1*max(max(max(Wxpn))),4),'k','LineWidth',1.1)
cb = colorbar('vert');
xlab = xlabel('z', 'interpreter', 'tex');
ylab = ylabel('y', 'interpreter', 'tex');
set(xlab, 'FontName', 'cmmi10', 'FontSize', 20);
set(ylab, 'FontName', 'cmmi10', 'FontSize', 20);
h = get(gcf,'CurrentAxes');
set(h,'FontName','cmr10','FontSize',15,'xscale','lin','yscale','lin');
set(cb, 'FontName', 'cmr10', 'FontSize', 15);
hold off