Finding neutral axis for arbitrary 2d shape - aerofoil

I am interested in determining the neutral axis of an aerofoil. I have the shape separated into lower and upper side (attached as txt files). When compiled it should look similar to this. I would like to determine the neutral axis of the profile. I have found this finding the axis for least moment of inertia of an object in 2D binary image - MATLAB Answers - MATLAB Central (mathworks.com) which uses a picture but sadly I do not have the necessary toolbox and my data is in x,y coordinates.

Is there a way to calculate it from this? I would appreciate any help!

[Dyuman Joshi]

What is the definition of 'nuetral axis' in this context?

[Selina]

It is the principle axis of the product of inertia I believe. So basically when load is applied to the object, where stress is zero.

Answer 1

[Mathieu NOE]

hello @Selina

the equations are still valid , see below how we can use them :

result

code :

u = readmatrix('upper.txt');
l = readmatrix('lower.txt');

x = [u(:,2);l(:,2)];
y = [u(:,3);l(:,3)];
n = numel(x, 1);% number of points

% centroids 
xc = mean(x);
yc = mean(y);

% compute moment of inertia
Ixx = sum(x.^2) / n;
Iyy = sum(y.^2) / n;
Ixy = sum(x.*y) / n;

% compute (ellipse) semi-axis lengths
common = sqrt( (Ixx - Iyy)^2 + 4 * Ixy^2);
ra = sqrt(2) * sqrt(Ixx + Iyy + common);
rb = sqrt(2) * sqrt(Ixx + Iyy - common);

% axes angle
theta = atan2(2 * Ixy, Ixx - Iyy) / 2;

% draw the neutral line 
slope = tan(theta);
xl = [min(x) max(x)];
yl = yc + slope*(xl-xc);
plot(x,y,'r',xc,yc,'dk',xl,yl,'b--');

[Selina]

Thank you for the help! I think the calculation to find the centroid might not be accurate. I believe for a composite part (unique shape, rather than elementary parts such as rectangles and triangles) it is not simply the mean of the x and y components. Statics: Centroids using Composite Parts (engineeringstatics.org)

But maybe I am wrong?

Edit: I managed to upload my coordinates as a binary figure and try the graphic version. It seems to provide a different centre. So I believe the centroid is calculated different.

[Mathieu NOE]

that's true indeed !

I was a bit too fast on this item

see correction below

the correct values are

xc = 0.4245

yc = 0.0182

instead of

xc = 0.5000

yc = 0.0126

u = readmatrix('upper.txt');
l = readmatrix('lower.txt');

x = [u(:,2);l(:,2)];
y = [u(:,3);l(:,3)];
n = numel(x, 1);% number of points

% centroids 
polyin = polyshape(x,y);
[xc,yc] = centroid(polyin);

% compute moment of inertia
Ixx = sum(x.^2) / n;
Iyy = sum(y.^2) / n;
Ixy = sum(x.*y) / n;

% compute (ellipse) semi-axis lengths
common = sqrt( (Ixx - Iyy)^2 + 4 * Ixy^2);
ra = sqrt(2) * sqrt(Ixx + Iyy + common);
rb = sqrt(2) * sqrt(Ixx + Iyy - common);

% axes angle
theta = atan2(2 * Ixy, Ixx - Iyy) / 2;

% draw the neutral line 
slope = tan(theta);
xl = [min(x) max(x)];
yl = yc + slope*(xl-xc);
plot(x,y,'r',xc,yc,'dk',xl,yl,'b--');

[Selina]

Oh yes that seems to align with the values I got from the graphical approach!

When I tried to use your code, it gave me the following error (see below). Did you modify the original input data in some way?

Edit: I got rid of this error message by removing the duplicate point at (0,0)

Warning: Polyshape has duplicate vertices, intersections, or other inconsistencies that may produce

inaccurate or unexpected results. Input data has been modified to create a well-defined polyshape.

> In polyshape/checkAndSimplify (line 523)

In polyshape (line 175)

In neutral_axis_2 (line 10)

Insufficient number of outputs from right hand side of equal sign to satisfy assignment.

Error in neutral_axis_2 (line 11)

[xc,yc] = centroid(polyin);

Edit2: This error persits and when running the code step by step it says:

Insufficient number of outputs from right hand side of equal sign to satisfy assignment.

Error in neutral_axis_2 (line 11)

[xc,yc] = centroid(polyin);

[Mathieu NOE]

try again with this code (I needed clearvars first to avoid the error message, haven't yet figured out where the root cause is)

I haven't touched to your data files

I also corrected a small bug that has no impact on final result anyway

n = numel(x);% number of points

instead of

n = numel(x,1);% incorrect usage of numel

full code

clc
clearvars

u = readmatrix('upper.txt');
l = readmatrix('lower.txt');
x = [u(:,2);l(:,2)];
y = [u(:,3);l(:,3)];
n = numel(x);% number of points
% centroids 
polyin = polyshape(x,y);
[xc,yc] = centroid(polyin);
% compute moment of inertia
Ixx = sum(x.^2) / n;
Iyy = sum(y.^2) / n;
Ixy = sum(x.*y) / n;
% compute (ellipse) semi-axis lengths
common = sqrt( (Ixx - Iyy)^2 + 4 * Ixy^2);
ra = sqrt(2) * sqrt(Ixx + Iyy + common);
rb = sqrt(2) * sqrt(Ixx + Iyy - common);
% axes angle
theta = atan2(2 * Ixy, Ixx - Iyy) / 2;
% draw the neutral line 
slope = tan(theta);
xl = [min(x) max(x)];
yl = yc + slope*(xl-xc);
plot(x,y,'r',xc,yc,'dk',xl,yl,'b--');

[Selina]

Oh amazing the clearvars seemed to have worked. Yeah, I was certainly confused because one of the inputs for centroid is polyshape and it told me, I did not give it the correct input...

[Mathieu NOE]

glad it works now !

[Selina]

One follow up comment: when comparing the section properties to ones you obtain through a CAD software, the neutral axis actually does not align (the matlab code predicts a slope of 3.034° while the CAD provides 3.61°). When running it via python (based on Python module for section properties - All this (leancrew.com)), it actually provides the correct slope. All three methods show the correct centroid. I assume the inertia calculations might be varying between the methods. I have not yet found a way to modify the code but wanted to add this in case anyone does not have the CAD or python code to verify.

[Mathieu NOE]

I wonder if there is a sign mistake in the Python code or in our matlab code

numerically speaking I have the same value as the matlab code I provided - and that's normal as both codes rely on the same equations

theta_deg = 0.9299

theta_deg2 = -0.9299

and this is different from what you announce above (and from the picture it seems to me you are running the code on another set of data)

u = readmatrix('upper.txt');
l = readmatrix('lower.txt');
x = [u(:,2);l(:,2)];
y = [u(:,3);l(:,3)];
n = numel(x);% number of points

% centroids 
polyin = polyshape(x,y);
[xc,yc] = centroid(polyin);

% compute moment of inertia
Ixx = sum(x.^2) / n;
Iyy = sum(y.^2) / n;
Ixy = sum(x.*y) / n;

% % compute (ellipse) semi-axis lengths
common = sqrt( (Ixx - Iyy)^2 + 4 * Ixy^2);
ra = sqrt(2) * sqrt(Ixx + Iyy + common);
rb = sqrt(2) * sqrt(Ixx + Iyy - common);

% axes angle
theta = atan2(2 * Ixy, Ixx - Iyy) / 2;
theta_deg = 180/pi*theta

% section below from link : 
% https://leancrew.com/all-this/2018/01/python-module-for-section-properties/
%'Principal moments of inertia (I1 I2) and orientation.'
I_avg = (Ixx + Iyy)/2;
I_diff = (Ixx - Iyy)/2;
I1 = I_avg + sqrt(I_diff^2 + Ixy^2);
I2 = I_avg - sqrt(I_diff^2 + Ixy^2);
theta2 = atan2(-Ixy, I_diff)/2;
theta_deg2 = 180/pi*theta2

% draw the neutral line 
slope = tan(theta);
xl = [min(x) max(x)];
yl = yc + slope*(xl-xc);
plot(x,y,'r*-',xc,yc,'dk',xl,yl,'b--');
ylim([-0.1 0.15])