1D Linear Convection (MATLAB code) (Updated with CFL condition)

clear; clc;% clear the screen and memory

Nx=41; %number of space nodes

Nt=82; %number of time nodes

Lx=2; %length of space (m)

Lt=1; %length of time (s)

dx=Lx/(Nx-1); %grid spacing

dt=Lt/(Nt-1); %time step

c=2; %speed of wave (constant)

a=c*dt/dx;

u=ones(Nx,1); %initialization of solution matrix

x=zeros(Nx,1); %initialization of space

hold on

for i=1:Nx-1 %space loop

    x(i+1)=x(i)+dx;

end

for i=0.5/dx:1/dx %initial conditions

    u(i)=2;

end

for t=1:Nt %time loop

    un=u; %u(i,t)

    for i=2:Nx %solution loop, backward in space forward in time

        u(i)=un(i)-a*(un(i)-un(i-1)); %discretized equation, u(i,t+1)

    end

    plot(x,u,'k') %plotting

    pause(0.1)

end



New Code, with the CFL condition satisfied.

clear; clc;% clear the screen and memory

Nx=13; %number of space nodes

Lx=2; %length of space (m)

Lt=1; %length of time (s)

dx=Lx/(Nx-1); %grid spacing

c=2; %speed of wave (constant)

dt=dx*0.1/c; %time step, Courant number = 0.1

Nt=(Lt/dt)+1; %number of time nodes

a=c*dt/dx;

u=ones(Nx,1); %initialization of solution matrix

x=zeros(Nx,1); %initialization of space

t=zeros(fix(Nt),1); %initialization of space

hold on

for i=1:Nt-1

    t(i+1)=t(i)+dt;

end

for i=1:Nx-1 %space loop

    x(i+1)=x(i)+dx;

end

for i=0.5/dx:1/dx %initial conditions

    u(i)=2;

end

for j=0:dt:Lt %time loop

    un=u; %u(i,t)

    for i=2:Nx %solution loop, backward in space forward in time

        u(i)=un(i)-a*(un(i)-un(i-1)); %discretized equation, u(i,t+1)

    end

    plot(x,u,'k') %plotting

    pause(0.1)

end

Solve PDEs using MATLAB (Code)

This is a post about a code I wrote to solve  1-D Poisson's equations. The method employed is called the finite difference method. For 1-D Laplace equation, the RHS of the equation mentioned below goes to zero. Hence all the terms in the B(i) loop equal zero. Also, B(1) and B(Nx) terms equal to the boundary conditions only. I hope this code helps you in learning something.

The example equation is mentioned below.
d²y/dx² = x²+1 and y(0)=10, y(1)=20

The result of the code is a plot between the dependent and the independent variables, as shown in Fig. 1. The accuracy of the solutions lie around 1.5 % of the analytical solution.

%solution to 1D Poisson's equations with dirichlet boundary conditions, 2nd order accurate, uniform grid
clear; clc;
Nx=101; %Enter the number of grid points
clc;
A=zeros(Nx,Nx); %Stiffness matrix
B=zeros(Nx,1); %Dependent variable vector, also used for boundary conditions and other known terms
x=zeros(Nx,1); %Independent variable vector
x(Nx)=1; %Domain size
dx=x(Nx)/(Nx-1); %Grid spacing
B(1)=(((x(1).^2)+1)*(dx.^2))-10; %Enter the function the PDE equals to multiplied by the square of the grid spacing minus the first boundary condition
B(Nx)=(((x(Nx).^2)+1)*(dx.^2))-20; %Enter the function the PDE equals to multiplied by the square of the grid spacing minus the second boundary condition
for i=2:Nx-1 %Off boundary values for the dependent variable vector
    B(i)=((x(i).^2)+1)*(dx.^2); %Enter the function the PDE equals to multiplied by the the square of the grid spacing
end
for i=1:Nx-1 %Loop for independent variable
    x(i+1)=(i)*dx;
end
for i=1:Nx %Diagonal terms for the stiffness matrix
    A(i,i)=-2;
end
for i=1:Nx-1 %Off diagonal terms for the stiffness matrix
    A(i,i+1)=1;
    A(i+1,i)=1;
end
X=A\B;
hold on
set(gca,'FontSize',28)
xlabel('Independent Variable','FontSize',38.5)
ylabel('Dependent Variable','FontSize',38.5)
title('Dependent VS Independent Variable (d^2y/dx^2 = x^2+1)','FontSize',38.5);
plot(x,X)

Fig. 1. Plot from the code

Solve ODEs using MATLAB (Code)

I have decided to teach myself some CFD coding using my own intuition and some free online tutorials. This is a post about code I wrote to solve a first order ordinary differential equation. The method employed is called the finite difference method. I will start by coding equations using this method and will slowly move to finite volume method. I hope this code helps you in learning something.

The equation in question is mentioned below.
dy/dx = x+4 and y(0)=3

A uniform grid was employed, as I have not yet figured out how to make non uniform grid. When I figure it out, I will update it here, which will be soon.

The code, mentioned below, starts with a variable called Nx, which represents the total number of grid points in which the domain is divided in to. A rule of thumb, the more points in the domain, more accurate the solution becomes at the cost of computation resources, mainly CPU time and RAM and storage.

x and y are the vectors in which the known and the unknown values are stored. x has the known values, like the length of a pipe through which a fluid flow or a rod through which heat flows etc. y stores the unknown values like temperature of the rod or velocity of the fluid etc. The values are stored at nodes.

y(1) means first value in the vector y. which also happens to be our boundary condition. y is the value of the function at x=0. x=0 is the first value in in the x vector, written as x(1).

x(Nx-1) represents the domain size. Minimum domain size is usually zero e.g. starting point of a rod or a pipe and maximum can be any real number.

To calculate the grid spacing, represented by dx, the domain size was divided by the number of elements. In this example, we have four nodes, two at corners and two in-between, forming a total of three elements.

First order finite difference method was used to discretize the given equation. And a for loop was used for the solution. Using 4 grid points, the solution has an accuracy of 2.22% against the analytical solution of the equation which can be verified by hand calculations. Thank you for reading.

%solution to dy/dx=x+4 y(0)=3
%uniform grid
clear; clc;
Nx=4; %number of grid points
x=zeros(Nx,1); %x vector, known, goes from some known value to maximum value of the domain
y=zeros(Nx,1); %y vector, one value known due to boundary condition rest unknown
y(1)=3; %boundary condition
x(Nx-1)=1; %domain size
dx=x(Nx-1)/(Nx-1); %grid spacing
for i=1:Nx-1
    x(i+1)=(i)*dx;
    y(i+1)=(dx.*(x(i)+4))+y(i); %discretized differential equation
end

Wing with Leading Edge Tubercles: Update 01

Tubercles (bumps) added to the leading edge of the aeronautic wings are a class of vortex generators. As shown in the Figure 1.

Figure 1, Modified wing geometry

Figure 2, Baseline wing geometry

The wing's cross section is NACA 0012. The dimensions used for this analysis were as used by [Chen, J. H., S. S. Li, and V. T. Nguyen. "The effect of leading edge protuberances on the performance of small aspect ratio foils." 15th International Symposium on Flow Visualization. 2012]. The flow conditions were reproduced from the research paper mentioned. In the wing with the leading edge tubercles, the tubercle wavelength was kept at two-thirds of the wing span. The tubercle amplitude was kept at half of the wing chord.

Computational fluid dynamics analysis was conducted to analyze the effect of the leading edge tubercles on wing performance. For validation of the numerical methodology, the lift and drag produced by the wing was compared to the experiments conducted for the mentioned research paper at 25 degree angle of attack. The forces calculated using computational fluid dynamics were within 7.3 % and 1.4 % of the experimental values, respectively.

The addition of tubercles at the leading edge of a wing creates a non-uniform pressure distribution on the wing’s surface, as shown in Figure 3. This non uniform pressure distribution and the changed wing geometry is responsible for the creation of a counter rotating chord-wise vortex pair behind each tubercle trough, as shown in Figure 4. To read more about non uniform pressure distribution, refer to Non Uniform Pressure Distribution.

Figure 3, Non-uniform pressure distribution. L-R; Modified Wing, Baseline Wing

Figure 4, Counter-rotating chord-wise vortex formation. L-R; modified wing, baseline wing. Notice the absence of vortex structures, represented by circles in the modified wing, in the baseline wing.

These vortices re-energize the boundary layer between them by carrying high momentum flow close to the wing surface, as shown by the dynamic pressure distribution in Figure 5, which leads to a delay in stall for the wing with leading edge tubercles, as shown in Figure 6. It can be seen clearly from both Figure 5-6 that the flow is attached behind the tubercle crest in the wing with leading edge tubercles while the baseline wing had stalled completely. To read more, refer to Counter Rotating Chord-wise Vortex Formation and Reduction in the Blade Tip Vortices, Delayed Stall and Reduction in the Span-wise Flow.

Figure 5, Dynamic pressure distribution around the wings. Top; modified wing. L-R; tubercle crest, tubercle trough. Bottom; baseline wing, same span-wise locations where the tubercle crest and trough is for the modified wing.

Figure 6, Streamlines around the wings. L-R baseline wing, modified wing.

Another effect of adding the leading edge tubercles to the wing is reduced span-wise flow and wing-tip vortices. The reason for this is that the counter rotating chord-wise vortices generated as a result of adding tubercles to the wing geometry act as a barrier to any span-wise flow, as shown in Figure 7-8.

Figure 7, Velocity flow trajectories. L-R modified wing, baseline wing. Notice the reduced size and strength of the tip vortex for the modified wing in comparison with the baseline wing.

Figure 8, Streamlines on the surface of the wings. L-R; modified wing, baseline wing. Reduced span-wise flow, (T-B of the screen), is clearly visible in the modified wing.

Thank you for reading. I hope this post added to your knowledge about tubercles. Next up will be the results of the counter-rotating configuration for the NREL Phase VI wind turbine. For verification and validation of the said wind turbine numerical simulations, refer to Verification and Validation.

Some of the reviewed literature:

[1] Watts, P., and Fish, F. E., “The Influence of Passive, Leading Edge Tubercles on Wing Performance,” Proceedings of the Unmanned Untethered Submersible Technology (UUST01), 2001.

[2] Fernandes, Irelyn, Yogesh Sapkota, Tania Mammen, Atif Rasheed, Calvin Rebello, and Young H. Kim, "Theoretical and Experimental Investigation of Leading Edge Tubercles on the Wing Performance," Proceedings of the Aviation Technology, Integration, and Operations Conference, Los Angeles, CA, 2013.
doi.org/10.2514/6.2013-4300
[3] Frank E. Fish, Paul W. Weber, Mark M. Murray, Laurens E. Howle, “The Tubercles on Humpback Whales' Flippers: Application of Bio-Inspired Technology,” Integrative and Comparative Biology, Vol. 51, No. 1, 2011, pp. 203–213.
doi.org/10.1093/icb/icr016
[4] Pedro, H. T. C., and Kobayashi, M. H., “Numerical Study of Stall Delay on Humpback Whale Flippers,” 46th AIAA Aerospace Sciences Meeting and Exhibit, AIAA Paper 2008-0584, Reno, NV, 2008.
doi.org/10.2514/6.2008-584
[5] Fish, F. E., and Lauder, G. V., “Passive and Active Flow Control by Swimming Fishes and Mammals,” Annual Review of Fluid Mechanics, Vol. 38, 2006, pp. 193–224.
doi.org/10.1146/annurev.fluid.38.050304.092201
[6] Weber, P. W., Howle, L. E., Murray, M. M., & Miklosovic, D. S., “Computational Evaluation of the Performance of Lifting Surfaces with Leading-Edge Protuberances,” Journal of Aircraft, Vol. 48, No. 2, 2011, pp. 591-600.
doi.org/10.2514/1.C031163
[7] Miklosovic, D. S., Howle, L. E., Murray, M. M., & Frank E. Fish, “Leading-edge tubercles delay stall on humpback whale (Megaptera novaeangliae) flippers,” Physics of Fluids, Vol. 16, No. 5, 2004.
dx.doi.org/10.1063/1.1688341
[8] N. Rostamzadeh, K. L. Hansen, R. M. Kelso & B. B. Dally, “The formation mechanism and impact of streamwise vortices on NACA 0021 airfoil's performance with undulating leading edge modification,” Physics of Fluids, Vol. 26, No. 10, 2014.
dx.doi.org/10.1063/1.4896748
[9] K. L. Hansen, R. M. Kelso & B. B. Dally, “Performance Variations of Leading-Edge Tubercles for Distinct Airfoil Profiles,” AIAA Journal, Vol. 49, No. 1, 2011.
doi.org/10.2514/1.J050631
[10] Ernst A. van Nierop, Silas Alben, and Michael P. Brenner, “How Bumps on Whale Flippers Delay Stall: An Aerodynamic Model,” Physical Review Letters, Vol. 100, No. 5, 2008.
doi.org/10.1103/PhysRevLett.100.054502
[11] K.L. Hansen, R.M. Kelso, B.B. Dally, E.R. Hassan, “Analysis of the Streamwise Vortices Generated Between Leading Edge Tubercles,” Proceedings of 6th Australian Conference on Laser Diagnostics in Fluid Mechanics and Combustion, Canberra, 2011
[12] Shi, Weichao, Mehmet Atlar, Rosemary Norman, Batuhan Aktas, and Serkan Turkmen, "Numerical optimization and experimental validation for a tidal turbine blade with leading-edge tubercles," Renewable Energy, Vol. 96, Part A, 2016, pp. 42-55.
doi.org/10.1016/j.renene.2016.04.064
[13] Ri-Kui Zhang and Jie-Zhi Wu, “Aerodynamic characteristics of wind turbine blades with a sinusoidal leading edge,” Wind Energy, Vol. 15, No. 3, 2012, pp 407-424.
doi.org/10.1002/we.479
[14] Ibrahim, I. H., and T. H. New, "A numerical study on the effects of leading-edge modifications upon propeller flow characteristics," Proceedings of Ninth International Symposium on Turbulence and Shear Flow Phenomena. Melbourne, 2015.
[15] Moore, K., and A. Ning, "Aerodynamic Performance Characterization of Leading Edge Protrusions," 54th AIAA Aerospace Sciences Meeting, AIAA Paper 2016-1786, San Diego, CA, 2016.
doi.org/10.2514/6.2016-1786
[16] Ning, Zhe, and Hui Hu, "An Experimental Study on the Aerodynamics and Aeroacoustic Characteristics of Small Propellers of UAV," 54th AIAA Aerospace Sciences Meeting, AIAA Paper 2016-1785, San Diego, CA, 2016.
doi.org/10.2514/6.2016-1785

Update 01

Following are my publications relating to the subject of this post.

Butt, F.R., and Talha, T., "A Numerical Investigation of the Effect of Leading-Edge Tubercles on Propeller Performance," Journal of Aircraft. Vol. 56, No. 2 or No. 3, 2019, pp. XX. (Issue/page number(s) to assigned soon. Active DOI: https://arc.aiaa.org/doi/10.2514/1.C034845)

Butt, F.R., and Talha, T., "A Parametric Study of the Effect of the Leading-Edge Tubercles Geometry on the Performance of Aeronautic Propeller using Computational Fluid Dynamics (CFD)," Proceedings of the World Congress on Engineering, Vol. 2, Newswood Limited, Hong Kong, 2018, pp. 586-595, (active link: http://www.iaeng.org/publication/WCE2018/WCE2018_pp586-595.pdf).

Butt, F.R., and Talha, T., "Optimization of the Geometry and the Span-wise Positioning of the Leading-Edge Tubercles on a Helical Vertical-Axis Marine Turbine Blade ," AIAA Science and Technology Forum and Exposition 2019, Turbomachinery and Energy Systems, accepted for publication.