# 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

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 — Laplacian .

The boundary conditions are given by

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

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

where

### Matlab codes

Determine the most amplified flow structures for the 3D linearized Navier-Stokes equations in a pressure-driven channel flow with , , and .

```% 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,:))');

% 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
```