CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > General Forums > Main CFD Forum

Generating elliptic grid around NACA SC2-0714 airfoil

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree5Likes

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   June 12, 2015, 18:49
Default Generating elliptic grid around NACA SC2-0714 airfoil
  #1
New Member
 
Justin
Join Date: Jan 2015
Posts: 28
Rep Power: 11
DA6righthand is on a distinguished road
I'm attempting to generate an elliptical grid around the said airfoil. Please keep in mind I'm a newbie to CFD so if you could keep responses as simple as possible, I'd super appreciate it! I'm using MATLAB for generating the grid and other calculations. I have several questions though. From what I'm told, it's reasonable to have the grid extend out in all directions roughly fifteen times the chord length. Attached is a plot of the airfoil shifted on the x-axis in order to make room for the elliptic grid. The plot was generated via coordinates found online. Where I'm having trouble is actually generating the grid by only knowing the coordinates of the airfoil.

What I know is that I'll have a predetermined number of grid points in both the plane and plane. The plane will have a rectangular grid representing the airfoil where the Navier-Stokes equations will be solved (which I'm not concerned about yet). What I've done thus far is use the following equations for the grid transformation:




where , , and . I then will use finite differences and solve for x^{n+1}_{i,j} and y^{n+1}_{i,j} where n represents the iteration. Using central differences for \alpha, \beta, and \gamma, there are no x or y at grid point i,j...so I'll spare myself typing the finite difference representations. So, solving the first two equations in their finite difference equivalents at point i,j we get




I also know I'll be using a loop to determine each x and y at point i,j and terminate the loop once the above differences have converged on a point. I'm just not sure about what to do next. How do I generate the elliptical grid by only knowing the x and y values of the airfoil? I've Googled many different things, searched this forum itself, tried using the Joe Thompson book on grid generation on Mississippi State University's website, none of which have truly helped. Any advice is greatly appreciated. I'm not looking for an outright answer...but rather where to look or some hints. I feel like I'm close but being that this is my first endeavor with grid generation, it's pretty tough.
Attached Images
File Type: jpg NASASC0714AirfoilPlot.jpg (11.2 KB, 131 views)

Last edited by DA6righthand; June 27, 2015 at 11:23.
DA6righthand is offline   Reply With Quote

Old   June 15, 2015, 10:10
Default
  #2
Member
 
A. S.
Join Date: Apr 2009
Location: Raipur (INDIA)
Posts: 54
Rep Power: 17
apoorv is on a distinguished road
Hi
First you have make grid using Transfinite Interpolation, between airfoil and outer domain. Then you have to fix the boundary vertices and march with elliptic grid generation method to get better grid. Also recheck the formulation of finite difference you are using more specifically the last term. It should be i,j+1 - i,j-1 if I remember correctly. If you need further help drop me a message.

Regards
Apoorv
DA6righthand likes this.
apoorv is offline   Reply With Quote

Old   June 15, 2015, 11:20
Default
  #3
New Member
 
Justin
Join Date: Jan 2015
Posts: 28
Rep Power: 11
DA6righthand is on a distinguished road
Quote:
Originally Posted by apoorv View Post
Also recheck the formulation of finite difference you are using more specifically the last term. It should be i,j+1 - i,j-1 if I remember correctly...
Thank you! And you're correct, I just typed it in wrong because it was a lot of copy and pasting.

I'll look up transfinite interpolation and see if I get anywhere. I appreciate it!
DA6righthand is offline   Reply With Quote

Old   June 16, 2015, 01:01
Default
  #4
Senior Member
 
Join Date: Oct 2011
Posts: 242
Rep Power: 16
naffrancois is on a distinguished road
As said by apoorv, in order to run an elliptic solver you need to prescribe both inner and outer boundaries and provide an initial "guess" for the interior nodes. The need for an accurate initial guess will strongly depend on the robustness of your elliptic solver but usually transfinite interpolation is necessary.

On the other side, hyperbolic mesh generation does not need information coming from the interior (thanks to the hyperbolicity) thus the algorithm generates successive layers from the only information coming from the airfoil coordinates (note that in this case it is difficult to get control on outer boundary node distribution). Though ,this method naturally imposes orthogonality and is order of magnitude faster computationaly. Look for Sorensen technical paper if interested
apoorv, DA6righthand and hamayun like this.
naffrancois is online now   Reply With Quote

Old   June 16, 2015, 17:25
Default
  #5
New Member
 
Justin
Join Date: Jan 2015
Posts: 28
Rep Power: 11
DA6righthand is on a distinguished road
I've done a bit of reading on transfinite interpolation and from what I understand is that I need a set of curves in order to set up the interpolation equation. Suppose I enclose the airfoil with a basic circle. Then I'd have curves enclosing the airfoil on the top and bottom...but how to I account for the airfoil in the center of the circle? Every internet resource I've looked seems to indicate that transfinite interpolation needs four curves enclosing the domain: top, bottom, left, and right.

I read the Thompson grid generation book sections on Lagrange and transfinite interpolation and his notation is a little confusing/hard to read but I think I'm understanding it a bit better.
DA6righthand is offline   Reply With Quote

Old   June 17, 2015, 00:53
Default
  #6
Senior Member
 
Join Date: Oct 2011
Posts: 242
Rep Power: 16
naffrancois is on a distinguished road
You are right you need a set of four curves, but this is also what you need for your elliptic solver. The definition of these four boundaries will define the type of grid you are constructing, have a look to C-grid, O-grid and H-grid
naffrancois is online now   Reply With Quote

Old   June 24, 2015, 12:12
Default
  #7
New Member
 
Justin
Join Date: Jan 2015
Posts: 28
Rep Power: 11
DA6righthand is on a distinguished road
I figured I'd let you know how I generated my grid. I searched about transfinite interpolation and couldn't figure it out. So here's what I did. I have the coordinates of the airfoil and the circle enclosing it. So between each coordinate of the airfoil and circle, starting from the trailing edge of the airfoil and 0*pi, I created a predetermined amount of linearly spaced x-values and used linear interpolation to find each new y coordinate. It worked just fine, as did generating the elliptic grid.

The only issue I'm having now is the amount of time it takes to generate the elliptic grid. Right now, with 42025 i,j points it takes 4 hour and 50 minutes to form the elliptic grid. Currently I'm using a while loop with two for loops within it. The for loops solve the finite difference equations in my first post for the elliptic grid. The while loop keeps iterating the for loops until the error between the x or y coordinate matrices reach a certain level (set at 0.00001). The way I'm evaluating the error is the Frobenius norm. For example, the error for the x matrix is



The error decreases very slow. Is there a more efficient way to evaluate error? Should I use a different norm to evaluate error? I'm looking for anything to speed up generating the grid. If anyone wants to see my code just let me know...I'll post it.
DA6righthand is offline   Reply With Quote

Old   June 24, 2015, 15:27
Default
  #8
Senior Member
 
Join Date: Oct 2011
Posts: 242
Rep Power: 16
naffrancois is on a distinguished road
I'll help you figure out how to implement transfinite interpolation, I did it some time ago, just need to find the code.

About tolerance, I do not see a critical problem, maybe I would have taken the Linfinity norm or would have scaled your actual norm by the total number of nodes. Anyway, if your solver is slow it is likely due to the way you solve the discrete equations. Could you post your code ?

Cheers
naffrancois is online now   Reply With Quote

Old   June 25, 2015, 13:59
Default
  #9
Senior Member
 
Join Date: Oct 2011
Posts: 242
Rep Power: 16
naffrancois is on a distinguished road
Hi, I found the routine I used to compute TFI as an initial guess for the elliptic solver.

As for the elliptic solver itself, here is the sequence I used:

- TFI
- elliptic solver without forcing terms
- elliptic solver + forcing terms

I solved the discrete equations using Multigrid + gauss-seidel (I am not saying it is the best option though, but computational time was acceptable, but still magnitude slower than the hyperbolic solver however)
Attached Files
File Type: f tfi.f (2.5 KB, 135 views)
hamayun likes this.
naffrancois is online now   Reply With Quote

Old   June 27, 2015, 11:12
Default
  #10
New Member
 
Justin
Join Date: Jan 2015
Posts: 28
Rep Power: 11
DA6righthand is on a distinguished road
naffrancois, thanks for your replies. Here is what I did in MATLAB (the only way I can fluently code). Between each point on the airfoil and enclosing circle I designated a set amount of x values linearly spaced apart. Then I used simple linear interpolation to find the y values. This generated an algebraic grid but when I run my elliptic solver, the elliptic grid is not reasonable. In other words, the elliptic grid looks good at the leading edge (sufficiently dense grid points) but the trailing edge of the airfoil is not dense enough. If you're able to run MATLAB files, attach is my code. If you cannot run MATLAB, here is my code below.

I'm working on another code also. The limitation to my attached code is that it is limited to a set number of grid points. The .txt file containing the airfoil coordinates has 205 x-values and 205 y-values. The way I wrote my code is the enclosing circle will have 205 x-value and 205 y-values. Conveniently, the airfoil coordinates begin at the trailing edge and progress counterclockwise around the airfoil. I generated the circle in the same manner beginning at 0\pi and ending 2\pi. That way when generating the algebraic grid, it begins at the trailing edge of the airfoil and the point at 0\pi and increments counterclockwise.

I'm just not entirely sure how to write a code that allows for an alternate number of grid points. I though about using several MATLAB builtin functions such as polyfit and polyval to generate a polynomial that fits the airfoil. This would allow any x and y in the domain of the airfoil to be evaluated. Although, doing so would introduce some amount of error. I experimented with the polyfit functions and was able to produce polynomials that very closely fit the original airfoil. In order to do this, I had to make two polynomials to represent the airfoil: one to plot the top half of the airfoil, the other to represent the bottom of the airfoil. MATLAB tells me that the polynomials are "poorly conditioned" though because I had to use polynomials of high degree (degree 19 for the top of airfoil, 20 for the bottom of airfoil). Other than what I described, I had no idea how to write a code that allows for an alternate number of grid points. My current code generates 205 by 205 x and y matrices. That yields 42025 grid points. I imagine I'll want more grid points when performing calculations for compressible supersonic flow.

Code:

function AirfoilTrans

B=load('NASASC20714AirfoilCoordinatesSelig.txt');

x1=B(1:end,1);
y1=B(1:end,2);

%center airfoil at the origin
for i=1:length(x1)
x1(i) = x1(i) - 0.5;
end
R=5;
L=length(y1);

theta=linspace(0,2*pi,L);
x2=R*cos(theta);
y2=R*sin(theta);

% figure(1)
% plot(x1,y1,x2,y2)

x=zeros(L,L);
y=x;
z=y;

for i=1:L
x(i,1)=x1(i);
x(i,end)=x2(i);
y(i,1)=y1(i);
y(i,end)=y2(i);
end

%Create algebraic grid by linearly spaced elements between known x
%coordinates on airfoil and circle then linearly interpolate to find values
%of y
for i=1:L
deltaX=linspace(x(i,1),x(i,end),L);
for j=2:L-1
x(i,j)=deltaX(j);

y(i,j)=y(i,1)+(y(i,end)-y(i,1))*((x(i,j)-x(i,1))/(x(i,end)-x(i,1)));
end
end

figure(1)
surf(x,y,z)

errX=1;
errY=1;
err = 0.00001;

xold=x;
yold=y;

xi=linspace(0,1,L);
dxi=xi(2)-xi(1);

eta=linspace(1,abs(R/max(y1)),L);
deta=eta(2)-eta(1);

iter=0;

while errX > err || errY > err

for i=2:L-1
for j=2:L-1
alpha1= (1/(2*deta))*(xold(i,j+1)-x(i,j-1));
alpha2= (1/(2*deta))*(yold(i,j+1)-y(i,j-1));

gamma1 = (1/(2*dxi))*(xold(i+1,j)-x(i-1,j));
gamma2 = (1/(2*dxi))*(yold(i+1,j)-y(i-1,j));

Alpha = alpha1^2 + alpha2^2;
beta = alpha1*gamma1 + alpha2*gamma2;
gamma = gamma1^2 + gamma2^2;

factor = 1/(4*(Alpha*deta^2 + gamma*dxi^2));

x(i,j) = factor*(2*Alpha*deta^2*(xold(i+1,j)+x(i-1,j))-beta*dxi*deta*(xold(i+1,j+1)-xold(i+1,j-1)-xold(i-1,j+1)+x(i-1,j-1)) + 2*gamma*dxi^2*(xold(i,j+1)+x(i,j-1)));

y(i,j) = factor*(2*Alpha*deta^2*(yold(i+1,j)+y(i-1,j))-beta*dxi*deta*(yold(i+1,j+1)-yold(i+1,j-1)-yold(i-1,j+1) +y(i-1,j-1)) + 2*gamma*dxi^2*(yold(i,j+1)+y(i,j-1)));
end
end
iter=iter+1;
if mod(iter,10)==0
iter
errX = norm(xold-x,1)/norm(xold,1)
errY = norm(yold-y,1)/norm(yold,1)
end

xold=x;
yold=y;
end

figure(2)
surf(x,y,z)
Attached Files
File Type: txt NASASC20714AirfoilCoordinatesSelig.txt (4.4 KB, 139 views)
DA6righthand is offline   Reply With Quote

Old   June 30, 2015, 01:50
Default
  #11
Senior Member
 
Join Date: Oct 2011
Posts: 242
Rep Power: 16
naffrancois is on a distinguished road
Hello, sorry for the late reply, I am pretty busy at the moment and will give you a more complete answer in the coming days.

About node distribution, you are right, you should have two distinct node distributions, one describing the geometry (the input coordinates describing the airfoil) and your mesh boundaries. Have a look to cubic splines interpolation. It does not try to fit a single polynomial to the whole set of airfoil coordinates, instead it builds piecewise cubic polynomials between each nodes.

You could also try simple linear interpolation between each airfoil coordinate nodes at first, but keeping in mind that the error will be high
naffrancois is online now   Reply With Quote

Old   June 30, 2015, 11:27
Default
  #12
New Member
 
Justin
Join Date: Jan 2015
Posts: 28
Rep Power: 11
DA6righthand is on a distinguished road
Linear interpolation is kind of what I did. I placed a certain amount of x values between each coordinate on the airfoil and circle, then used basic linear interpolation to determine the y coordinates of the algebraic grid.

So you think a spline would work better the fitting a high degree polynomial to the airfoil? I'll try this later today and see what happens. Thanks!
DA6righthand is offline   Reply With Quote

Old   June 30, 2015, 18:00
Default
  #13
Senior Member
 
Join Date: Oct 2011
Posts: 242
Rep Power: 16
naffrancois is on a distinguished road
I think we misunderstood each other, if not excuse me if you already understood the following. About linear interpolation or cubic splines, I was talking about how you could get any number of nodes on j=1 line (on your airfoil then).

For now, the code you gave me (which runs well by the way), allows for exactly the same number of nodes on j=1 (and as a consequence on j=jmax) than the number of coordinates you give in input. In order to allow for a more flexible mesh generation, the first step in to interpolate the input coordinates of your airfoil to a new set of nodes (on j=1). I believe that you cannot fit these coordinates by a single polynomial, instead either you can have a look at cubic splines, or can try to simply linearly interpolate between each airfoil coordinates (but you may have surprises when running later your CFD code on it).

Then, you have an independent j=1 node distribution and you can go for the generation of the interior mesh:

- algebraic initial guess (linear interpolation as you did is ok for your configuration because the convexity of your airfoil is not very high, or more robust guesses like TFI)
- elliptic smoothing

Then, you will have a look at forcing terms to enforce grid stretching and orthogonality.
naffrancois is online now   Reply With Quote

Old   June 30, 2015, 18:42
Default
  #14
Senior Member
 
Join Date: Oct 2011
Posts: 242
Rep Power: 16
naffrancois is on a distinguished road
About speed of the elliptic solver, one thing you should look at is the initial guess (algebraic grid). If you zoom close to the trailing edge, you will see crossing lines, in my opinion this may affect the rate of convergence badly
naffrancois is online now   Reply With Quote

Old   June 30, 2015, 18:47
Default
  #15
New Member
 
Justin
Join Date: Jan 2015
Posts: 28
Rep Power: 11
DA6righthand is on a distinguished road
About the trailing edge of the airfoil: I was told it might help to slightly alter the very last few coordinates of the airfoil and this could help.

If you ran my code, you'll see the elliptic grid is fairly condense at the leading edge while the trailing edge is not. Do forcing functions help control the grid density or is the problem arising from the algebraic grid?

I really appreciate all of your help, by the way. Thanks!
DA6righthand is offline   Reply With Quote

Old   July 1, 2015, 02:13
Default
  #16
Senior Member
 
Join Date: Oct 2011
Posts: 242
Rep Power: 16
naffrancois is on a distinguished road
I am not refering to the coordinates of the trailing edge, but to the algebraic grid shown in figure(1) of your matlab code, if you zoom near the trailing edge, right below you will see that lines are crossing each other. This might be problematic and slow down the elliptic solver.

About grid density, you are right, forcing functions will help a lot in clustering nodes where you want or impose orthogonality to the airfoil.
naffrancois is online now   Reply With Quote

Old   July 1, 2015, 16:14
Default
  #17
Senior Member
 
Join Date: Oct 2011
Posts: 242
Rep Power: 16
naffrancois is on a distinguished road
Here are some pictures of results of applying different approaches I did some time ago:
- TFI
- Laplace
- Laplace + forcing functions
- Hyperbolic

As you see, the TFI and Laplace without forcing terms give a poor quality grid. The reason is that they do not have any mechanism to enforce orthogonality and much control on the interior grid spacing. With forcing functions you can force the grid to some extend

I think I took as a reference the source codes coming with the book of Blazek:

- Computational fluids dynamics: principles and applications

grid_TFI.jpg

grid_Laplace.jpg

grid_Poisson.jpg

grid_Hyperboic.jpg
naffrancois is online now   Reply With Quote

Old   July 21, 2015, 19:07
Default
  #18
New Member
 
Justin
Join Date: Jan 2015
Posts: 28
Rep Power: 11
DA6righthand is on a distinguished road
Quote:
Originally Posted by naffrancois View Post
...As you see, the TFI and Laplace without forcing terms give a poor quality grid. The reason is that they do not have any mechanism to enforce orthogonality and much control on the interior grid spacing. With forcing functions you can force the grid to some extent...
Sorry for the long lapse in replies. I've been busy also with getting ready to implement the SIMPLE procedure for airflow around the same airfoil. I started looking into forcing functions and found something potentially useful in the CFD book by T. J. Chung. Similar to the equations I have in the original post, I'll need to solve for the x,y grid points using

\alpha x_{\xi\xi}-2\beta x_{\eta\xi}+\gamma x_{\eta\eta}=-J^2\left(Px_\xi+Qx_\eta\right)

\alpha y_{\xi\xi}-2\beta y_{\eta\xi}+\gamma y_{\eta\eta}=-J^2\left(Py_\xi+Qy_\eta\right)

where \alpha, \beta, and \gamma are as defined in the original post, J=\frac{\partial(x,y)}{\partial(\xi,\eta)}=x_\xi y_\eta-x_\eta y_\xi is the Jacobian, and P and Q are given by
P(\xi,\eta)=-\sum_{i=1}^{n}a_i\left | \xi-\xi_i\right |e^{-c_i\left | \xi-\xi_i \right |}-\sum_{i=1}^{m}b_i\left | \xi-\xi_i \right |e^{-d_i\sqrt{(\xi-\xi_i)^2+(\eta-\eta_i)^2}}

Q(\xi,\eta)=-\sum_{i=1}^{n}a_i\left | \eta-\eta_i\right |e^{-c_i\left | \eta-\eta_i \right |}-\sum_{i=1}^{m}b_i\left | \eta-\eta_i \right |e^{-d_i\sqrt{(\xi-\xi_i)^2+(\eta-\eta_i)^2}}

In P and Q, n is the number of \xi lines and m is the number of \eta lines. The amplification factors are a_i and b_i. According to the Chung book,
- a_i>0 lines \xi are attracted to lines \xi_i (similarly for \eta coordinates)
- b_i>0 lines \xi are attracted to points (\xi_i,\eta_i) (similarly for \eta coordinates)
- and the decay factors c_i,d_i modulate the amplification factors
So should I choose distinct amplification factors for each point? Or could I just make it a constant? Basically my goal now is to condense the \eta lines closer to the airfoil and condense \xi lines closer around the leading and trailing edges. But, being a newbie to CFD, I'm not sure how to manipulate these forcing functions. Although it seems the only way to enhance my current elliptic grid is by using these forcing functions.

My question is when I write the finite difference form with the forcing functions, will I have to sum to n and m for each iteration (i.e. for each i,j)? I'm a little confused on how to incorporate that into the differenced equations.

Last edited by DA6righthand; July 21, 2015 at 19:22. Reason: clicked submit reply rather than preview reply...I stoopid
DA6righthand is offline   Reply With Quote

Old   July 24, 2015, 17:01
Default
  #19
Senior Member
 
Join Date: Oct 2011
Posts: 242
Rep Power: 16
naffrancois is on a distinguished road
Hello, I did not use these forcing functions, I think I used the one based on orthogonality properties (it is mentionned in the book of Chung). Nevertheless, in order to test them, you should probably focus first on clustering at a single location. The first term clusters lines and the second term clusters towards points. You should try first to cluster only to eta=1 line rather than setting different clustering locations, as you will see you rapidly get an unstable solver.

Looking at the expressions for P and Q, it seems that you can precompute these quantities before the solving step, no need to recompute them as they do not depend on physical points.

Would not be a bad idea to test it on a simpler geometry like a square first

good luck and let me know !
naffrancois is online now   Reply With Quote

Old   November 2, 2015, 10:50
Default
  #20
Member
 
mechiebud
Join Date: Jan 2015
Posts: 49
Rep Power: 11
mechiebud is on a distinguished road
Quote:
Originally Posted by DA6righthand View Post
I'm attempting to generate an elliptical grid around the said airfoil. Please keep in mind I'm a newbie to CFD so if you could keep responses as simple as possible, I'd super appreciate it! I'm using MATLAB for generating the grid and other calculations. I have several questions though. From what I'm told, it's reasonable to have the grid extend out in all directions roughly fifteen times the chord length. Attached is a plot of the airfoil shifted on the x-axis in order to make room for the elliptic grid. The plot was generated via coordinates found online. Where I'm having trouble is actually generating the grid by only knowing the coordinates of the airfoil.

What I know is that I'll have a predetermined number of grid points in both the plane and plane. The plane will have a rectangular grid representing the airfoil where the Navier-Stokes equations will be solved (which I'm not concerned about yet). What I've done thus far is use the following equations for the grid transformation:




where , , and . I then will use finite differences and solve for x^{n+1}_{i,j} and y^{n+1}_{i,j} where n represents the iteration. Using central differences for \alpha, \beta, and \gamma, there are no x or y at grid point i,j...so I'll spare myself typing the finite difference representations. So, solving the first two equations in their finite difference equivalents at point i,j we get




I also know I'll be using a loop to determine each x and y at point i,j and terminate the loop once the above differences have converged on a point. I'm just not sure about what to do next. How do I generate the elliptical grid by only knowing the x and y values of the airfoil? I've Googled many different things, searched this forum itself, tried using the Joe Thompson book on grid generation on Mississippi State University's website, none of which have truly helped. Any advice is greatly appreciated. I'm not looking for an outright answer...but rather where to look or some hints. I feel like I'm close but being that this is my first endeavor with grid generation, it's pretty tough.
Hi I too am generating elliptic grid over an airfoil. I have generated the geometry, got the initial points using TFI but my code for the elliptic grid is not working as I am getting infinite values. Could you please guide me in this. Did you encounter any such error?
mechiebud is offline   Reply With Quote

Reply

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
SU2 AOA optimization 454514566@qq.com SU2 9 March 7, 2022 16:17
elliptic grid generation (orthogonal) vasu Main CFD Forum 8 October 28, 2015 15:20
[ICEM] O-grid mesh for NACA airfoil SR-71 ANSYS Meshing & Geometry 0 February 6, 2015 07:45
2D FFD Optimization RLangtry SU2 2 August 5, 2014 09:48
Structured Airfoil Grid Generator A.S. Main CFD Forum 1 September 30, 2008 18:29


All times are GMT -4. The time now is 03:29.