Finding neutral axis for arbitrary 2d shape - aerofoil

16 views (last 30 days)
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!
  2 Comments
Dyuman Joshi
Dyuman Joshi on 1 Dec 2023
What is the definition of 'nuetral axis' in this context?
Selina
Selina on 1 Dec 2023
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.

Sign in to comment.

Accepted Answer

Mathieu NOE
Mathieu NOE on 1 Dec 2023
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--');
  8 Comments
Selina
Selina on 5 Dec 2023
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
Mathieu NOE on 8 Dec 2023
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])

Sign in to comment.

More Answers (0)

Categories

Find more on Programming in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!