ProtoCYPVis.m

% Grid Points Visualizer
% Gauss Quadrature Points Mapped in Space in Each Direction using a
% Symmetric Chebyshev Distribution (Initial Guess to Legendre)
% Written by Stephen Tang

clear all
close all
clc

% Read u+ y+ data
D1yp = xlsread('uPyP.csv','uPyP','A1:A63');
D1up = xlsread('uPyP.csv','uPyP','B1:B63');
D2yp = xlsread('LogLaw.csv','LogLaw','A1:A70');
D2up = xlsread('LogLaw.csv','LogLaw','B1:B70');

% Progression Ratio
n1 = 0.86; % Horizontal (Airfoil L.E to Outer)
n2 = 1; % Vertical Outer Domain
n3 = 1.00; % Vertical Inner Domain
n4 = 1;  % Horizontal Outer Domain
n5 = 1.01; % Horizontal Inner Domain

% Number of Layers [Grid Spread]
m1 = 90; % Horizontal (Airfoil L.E to Outer)
m2 = 20; % Vertical Outer Domain
m3 = 100; % Vertical Inner Domain
m4 = 150; % Horizontal Outer Domain
m5 = 200; % Horizontal Inner Domain


% Domain 10
xD10 = 0.0685;
yD10 = 10;
% First Cell Size y near Airfoil
C1yD10 = yD10 * ((n1)^(m1-1));
% Second Cell Size y near Airfoil
C2yD10 = yD10 * ((n1)^(m1-2));
% Sum of 1st & 2nd Cell 
CTyD10 = C1yD10 + C2yD10;

% Calculate Domain of Individual Cell Height
k = 30; % Number of h-cells
CHT(1) = 0;
for i = 1:k
    CyHI(i) = yD10 * ((n1)^(m1-(i)));    
end

for i = 1:k
    CHT(i+1) = CHT(i) + CyHI(i);
end

% Initial Guess to Legendre Gauss where Chebyshev Gauss is Used
% Gauss quadrature points in the mapped space are placed in each direction
% using a symmetric Chebyshev distribution

P = 2; % Polynomial Order
Pa = linspace(0,P,P+1);
%Pa = Pa + 0.5;
for i = 1 : length(Pa)
X(i) = 0.5 * (1-cos((((i-1)+0.5)/(P+1))*pi)); % X_i+1/2
%X2(i) = 0.5 * (1-cosd((((i-1)+0.5)/(P+1))*pi)); % X_i+1/2
end

% Add in Points for P-order
CHPT(1) = CHT(1);
for i = 1:k
    for j = 1 : length(Pa)
        X(j) = 0.5 * (1-cos((((j-1)+0.5)/(P+1))*pi)); % X_i+1/2
        XAd(j) = X(j).*CyHI(i) + CHPT(i+((i-1)*length(Pa)));
        CHPT(i+((i-1)*length(Pa))+j) = XAd(j);
    end
CHPT(i+1+(i*length(Pa))) = CHT(i+1);  
end



%     for j = 1 : length(Pa)
%         X(j) = 0.5 * (1-cos((((j-1)+0.5)/(P+1))*pi)); % X_i+1/2
%         XAd(j) = X(j).*CyHI(1+1) + CHPT(1+length(Pa));
%         CHPT(1+(j) = XAd(j);
%     end
% CHPT(1+length(Pa)) = CHT(1+1);  

% Calculate y+

% Input Chord Length of Airfoil
Lbl = 1; % L_boundary layer

% Input Freestream Velocity
Ma = 0.15; % Mach Number

% Method 1: Calculate Based on Length of Chord 

% Assume Sea Level for Conversion
Uf = Ma * 343; % Speed of Sound = 343 m/s
rho = 1.225; % Density of Air [kg/m^3]
mu = 1.81 * (10^(-5)); % Viscosity of Air [kg/(m*s)]

% Method 2: Calculate Based on Reynolds Number
ReA = [1000, 3000, 30000, 100000, 1000000, 3000000, 6000000, 9000000]; % Reynolds Number Array

% Estimate Skin Friction [Flat Plate]

% Schlichting Skin-Friction
CfSchA = (2.*(log10(ReA))-0.65).^(-2.3); 

% Compute Wall Shear Stress:
TauWA = CfSchA.*0.5.*rho.*(Uf.^2);

% Compute Friction Velocity
UfricA = sqrt(TauWA./rho);

% Initiate Count Var
nLSA = 0;
nBLA = 0;
nLLA = 0;
nWRA = 0;

% Array of Dimensionless Distance y+
ypA = [1, 5, 10, 15, 30, 350, 10000];
for i = 1:length(ypA)
    if (0 <= ypA(i)) && (ypA(i) <= 5)
    nLSA = nLSA + 1; % Number of Elements in Laminar Sublayer
    elseif (5 < ypA(i)) && (ypA(i) <= 30)
    nBLA = nBLA + 1; % Number of Elements in Buffer Layer
    elseif (30 < ypA(i)) && (ypA(i) <= 350)
    nLLA = nLLA + 1; % Number of Elements in Log Layer
    elseif (ypA(i) > 350)
    nWRA = nWRA + 1; % Number of Elements in Wake Region
    end
end

% Calculate Wall Distance
for h = 1: length(UfricA)
    for j = 1:length(ypA)
yA(j,h) =  (ypA(j).*mu)./(rho.*UfricA(h));
    end
end
% Calculate Dimensionless Velocity u+
% up = U/Ufric 

% Laminar Sublayer
for i = 1:nLSA
upA(i) = ypA(i);
end

% Buffer Layer & Low Law Region
for i = nLSA:(nLSA + nBLA + nLLA)
upA(i) = 2.5*log(ypA(i)) + 5.45; 
end

% Turbulent Zone
for i = (nLSA + nBLA + nLLA):(nLSA + nBLA + nLLA + nWRA)
upA(i) = 5.6*log10(ypA(i)) + 4.9;
end

% Calculate y+
for h = 1:length(UfricA)
CHPTP(:,h) = (CHPT .* rho .* UfricA(h))./mu ;
end

% % Calculate y+
% for h = 1:length(UfricA)
% CHPTP(:,h) = (CHPT .* rho .* UfricA(h))./mu ;
% end

% Sort y+ into respective layers
% Initiate Count Var
nLSA2(1:length(UfricA)) = 0;
nBLA2(1:length(UfricA)) = 0;
nLLA2(1:length(UfricA)) = 0;
nWRA2(1:length(UfricA)) = 0;

for j = 1:length(UfricA)
for i = 1:length(CHPTP(:,j))
    if (0 <= CHPTP(i,j)) && (CHPTP(i,j) <= 5)
    nLSA2(j) = nLSA2(j) + 1; % Number of Elements in Laminar Sublayer
    elseif (5 < CHPTP(i,j)) && (CHPTP(i,j) <= 30)
    nBLA2(j) = nBLA2(j) + 1; % Number of Elements in Buffer Layer
    elseif (30 < CHPTP(i,j)) && (CHPTP(i,j) <= 350)
    nLLA2(j) = nLLA2(j) + 1; % Number of Elements in Log Layer
    elseif (CHPTP(i,j) > 350)
    nWRA2(j) = nWRA2(j) + 1; % Number of Elements in Wake Region
    end
end
end

for j = 1:length(UfricA)
% Laminar Sublayer
for i = 1:nLSA2(j)
uPTP(i,j) = CHPTP(i,j);
end

% Buffer Layer & Low Law Region
for i = nLSA2:(nLSA2 + nBLA2 + nLLA2)
uPTP(i,j) = 2.5*log(CHPTP(i,j)) + 5.45; 
end

% Turbulent Zone
for i = (nLSA2 + nBLA2 + nLLA2):(nLSA2 + nBLA2 + nLLA2 + nWRA2)
uPTP(i,j) = 5.6*log10(CHPTP(i,j)) + 4.9;
end

end
% Calculate u+

%figure(1)
%plot(X.*C1yD10,0,'*-');

figure(2)
semilogx(D1yp,D1up,'b-.')
hold on
semilogx(D2yp,D2up,'r-.')
semilogx(CHPTP(:,1),uPTP(:,1),'k-x')
title('Law of the Wall Plot')
xlabel('y+')
ylabel('u+')
legend('u+ = y+','u+ = 2.5ln(y+) + 5.45','Re = 1000 & p = 2')
hold off

figure(3)
semilogx(D1yp,D1up,'b-.')
hold on
semilogx(D2yp,D2up,'r-.')
semilogx(CHPTP(:,2),uPTP(:,2),'k-x')
title('Law of the Wall Plot')
xlabel('y+')
ylabel('u+')
legend('u+ = y+','u+ = 2.5ln(y+) + 5.45','Re = 3000 & p = 2')
hold off

figure(4)
semilogx(D1yp,D1up,'b-.')
hold on
semilogx(D2yp,D2up,'r-.')
semilogx(CHPTP(:,3),uPTP(:,3),'k-x')
title('Law of the Wall Plot')
xlabel('y+')
ylabel('u+')
legend('u+ = y+','u+ = 2.5ln(y+) + 5.45','Re = 30000 & p = 2')
hold off

figure(5)
semilogx(D1yp,D1up,'b-.')
hold on
semilogx(D2yp,D2up,'r-.')
semilogx(CHPTP(:,4),uPTP(:,4),'k-x')
title('Law of the Wall Plot')
xlabel('y+')
ylabel('u+')
legend('u+ = y+','u+ = 2.5ln(y+) + 5.45','Re = 100000 & p = 2')
hold off

figure(6)
semilogx(D1yp,D1up,'b-.')
hold on
semilogx(D2yp,D2up,'r-.')
semilogx(CHPTP(:,5),uPTP(:,5),'k-x')
title('Law of the Wall Plot')
xlabel('y+')
ylabel('u+')
legend('u+ = y+','u+ = 2.5ln(y+) + 5.45','Re = 1000000 & p = 2')
hold off

figure(7)
semilogx(D1yp,D1up,'b-.')
hold on
semilogx(D2yp,D2up,'r-.')
semilogx(CHPTP(:,6),uPTP(:,6),'k-x')
title('Law of the Wall Plot')
xlabel('y+')
ylabel('u+')
legend('u+ = y+','u+ = 2.5ln(y+) + 5.45','Re = 3000000 & p = 2')
hold off