Determine Camber and Thickness of a Airfoil - Given the xy coordinates

Dear All,

I have got airfoil coordinates, in the form of a closed polygon (xy points given below).

I need a function to determine the maximum possbile inscribed circle about each point (vertex) in the polygon.

ps: inscribed circle at a vertex point is defined as the maximum possible circle drawn inside the polygon and also tangent to the vertex.

The airfoil below given contains 33 points, and i require 33 inscribed circles...

Otherwise need to determine the camber and thickness distribution of the airfoil as referred in the comments below.

Any help is highly appreciated.

Thanks in advance,

K Vijay Anand

%Coordinates of Airfoil

xy = [1.00000 0.94908 0.89816 0.79638 0.69491 0.59377 0.49322 0.39328 0.29389 0.19489 0.14562 0.09650 0.07203 0.04764 0.02338 0.01143 0.00000 0.01330 0.02596 0.05095 0.07573 0.10048 0.14992 0.19932 0.29821 0.39752 0.49722 0.59742 0.69804 0.79883 0.89957 0.94979 1.00000;

0.00000 0.02835 0.05669 0.11138 0.15658 0.19180 0.20853 0.20678 0.18805 0.15733 0.13473 0.10764 0.09134 0.07255 0.04976 0.03287 0.00000 -0.02457 -0.02966 -0.02934 -0.02254 -0.01473 0.00237 0.02098 0.05519 0.07642 0.08566 0.07942 0.06019 0.03596 0.01324 0.00637 0.00000];

plot(xy(1,:),xy(2,:),'o-','LineWidth',2); grid on; axis equal;

[VBBV]

Could you elaborate what "maximum possible inscribed circle" means ?

[DGM]

I'll have to give it some thought, but I don't know that you're going to be able to solve the problem for all vertices. Only vertices where the polygon is not strictly convex will be able to be tangent to an inscribed circle (of nonzero radius). Consider that the inscribed circle of a square is tangent to none of the square's vertices.

EDIT:

Have you taken a look at this discussion?

https://stackoverflow.com/questions/3783207/determining-a-mean-camber-line

You may also be able to get some insight by seeing how this calculates the camber

https://www.mathworks.com/matlabcentral/fileexchange/12889-airfoil-analyzer

[Vijay Anand]

Thanks a lot Mr.DGM for your quick response !!!

I ve already gone through the second reference you have provided(Airfoil Analyzer), they have followed a simple method of taking average of top and bottom surface points to determine the mean camber line ( which is approximate and not good for high cambered airfoils) . . .

The idea being to get back the thickness and camber distribution for any arbitrary airfoil (with high camber) ...

The first reference also I have seen, but could not implement effectively the different approcahes suggested over there . . .

Inscribed Circle over square is a good take awway for me and Thanks again for your sincere inputs and will wait for others comments also ...

-Vijay

[Vijay Anand]

Dear Mr. @VBBV,

I am looking forward to create similar circles...

[DGM]

Ah. I didn't notice that they were solving it that way.

Answer 1

[DGM]

This isn't exactly a great way, but I'm more used to abusing image processing tools than polyshape() and such. The mention of Voronoi diagrams made me think that a distance map would be a decent place to start.

profileline = [1.00000 0.94908 0.89816 0.79638 0.69491 0.59377 0.49322 0.39328 0.29389 0.19489 0.14562 0.09650 0.07203 0.04764 0.02338 0.01143 0.00000 0.01330 0.02596 0.05095 0.07573 0.10048 0.14992 0.19932 0.29821 0.39752 0.49722 0.59742 0.69804 0.79883 0.89957 0.94979 1.00000;
    0.00000 0.02835 0.05669 0.11138 0.15658 0.19180 0.20853 0.20678 0.18805 0.15733 0.13473 0.10764 0.09134 0.07255 0.04976 0.03287 0.00000 -0.02457 -0.02966 -0.02934 -0.02254 -0.01473 0.00237 0.02098 0.05519 0.07642 0.08566 0.07942 0.06019 0.03596 0.01324 0.00637 0.00000];

stepsize = 0.0001; % effective resolution
xrange = [-0.1 1.1];
yrange = [-0.1 0.25];

% use a dummy image display to generate an image of the profile
x = xrange(1):stepsize:xrange(2);
y = yrange(1):stepsize:yrange(2);
h = image(xrange,yrange,ones(numel(y),numel(x)));
L = images.roi.Polygon(gca);
L.Position = profileline.';
mask = ~createMask(L);

% find camberline coordinates from the ridgeline of the distance map
% this is only an incomplete solution due to the leading edge angle
dmap = bwdist(mask);
[~,idx] = max(dmap,[],1);
camberline1 = [x; y(idx)];

% use the initial estimate to find where to bisect the profile
% so that the leading edge can be estimated 
[~,breakpoint] = max(camberline1(2,:));
breakpoint = round(breakpoint/2); % avoid shallow slopes
[~,idx] = max(dmap(:,1:breakpoint),[],2);
camberline2 = [x(idx); y];

% plot estimates for demonstration
subplot(2,1,1)
plot(profileline(1,:),profileline(2,:)); hold on; grid on
plot(camberline1(1,:),camberline1(2,:),'r:');
plot(camberline2(1,:),camberline2(2,:),'b:');
xlim(xrange)
ylim(yrange)

% clean up and merge estimates
camberline1(:,camberline1(2,:)==yrange(1)) = NaN;
camberline2(:,camberline2(1,:)==xrange(1)) = NaN;
camberline2(:,camberline2(1,:)==x(breakpoint)) = NaN;
camberline = [camberline2 camberline1(:,breakpoint+1:end)];

% plot merged estimate
subplot(2,1,2)
plot(profileline(1,:),profileline(2,:)); hold on; grid on
plot(camberline(1,:),camberline(2,:),'b--');
xlim(xrange)
ylim(yrange)

As mentioned, the leading edge area is problematic. I doubt either estimate is correct in that region, but the second looked marginally more reasonable, and I figured I'd at least demonstrate that a piecewise solution may be an option. It would take me more time to wrap my head around a more robust approach.

[Vijay Anand]

Dear Mr. @DGM

Thanks a lot your insight. I will try to understand the method you have developed, which uses image processing logic.

The robust method has to offer a solution which is inversible also. ie. once we find the camber and thickness distribution(measured normal to the camber line) from the profile line, the inverse function should able to reproduce the profile line from camber and thickness distributions.

From Airfoil Geometry Theory, the formula for finding profile lines { (xupper,yupper) and (xlower,ylower) } from thickness (t) and camber distribution(xmean,ymean),

The reference is taken from https://www.pdas.com/refs/aiaa5235.pdf

Thanks a lot again

K Vijay Anand

[DGM]

Well, again excepting the issue of the LE ambiguity, the above example should be invertible if modified to accomodate that intent. If the lines used to detect the ridgeline are adjusted, we can also get the value of the distance map at those locations. Since dmap is the euclidean distance map, each value of mx1 and mx2 are the radii of the largest inscribed circles at the given points along the respective paths.

[mx1,idx] = max(dmap,[],1);
% ...
[mx2,idx] = max(dmap(:,1:breakpoint),[],2);

I don't think that would give you the proper values for t, but it should at least be a means to reconstruct the profile line. That said, I'm not exactly suggesting that it's appropriate or practical. Just saying it's possible.

[Vijay Anand]

Dear Mr. @DGM and others,

Still thinking, is there a way to find this two way function or is not theorectically possible.

Say, we have a arbitrary shape ( closed polygon) with xy coordinates defined . is it possible to convert the polygon into a symetric polygon wrt x-axis (or some other axis) and a offset value ( like camber) ....

And need to ensure, the polygon is obtainable from the symmetry shape and camber line.

Thanks in advance,

K Vijay Anand

[DGM]

I imagine it's possible, but I'm really at a loss for a good approach, especially for a generalized solution. You might consider editing the title of your question to see if it can attract new answers.

[Vijay Anand]

Dear Mr.@DGM and others,

A small break through...

As i thought earlier, these circles need not be inscribed circles...

The definition of camber line says, the thickness should be normal to the camber line (measured at a point)... and also thickness is equal on top and bottom of the camber line.

Please refer to the figures 1,2,3

Figure 1
The dotted vertical lines are thickness lines plotted verticaly in y-axis. 
slanted lines are thickness lines plotted normal to local camber point...
The circles are drawn with thickness/2 as radius and center at camber point...

Figure 2

Figure 3

From the figures near leading edge (figure 2), it is evident that the circles are bigger . . .

Kindly give your valuable feed back on this...

Still, I am struck with the algorithm to determine the camber line and thickness distribution.

ps: these figures i have generated from known camber and thickness distributions equations (NACA 4312 Aifoil) and coordinate points

Thanks in advance,

K Vijay Anand

[nyhuma]

hey @Vijay Anand

can you provide the code for this analysis? i have implemented something similiar to this, but i am only looking for radii with r => "distance of voronoi_site to polygon". this way i am lacking voronoi sites at the leading and trailing edges of airfoils or in parts of high curvature. this makes the analysis of turbine airfoils pretty unreliable.

Answer 2

[Bruno Luong]

If you can tolerate some finite precision, you migh discretize the interior as b&w image, look for the skeleton

https://www.mathworks.com/help/images/ref/bwskel.html

For each point of the skeleton (center of the circle) the radius is the minimum of the distance from the boundary (not very difficult to compute).