Finding neutral axis for arbitrary 2d shape - aerofoil

13 次查看（过去 30 天）
Selina
Selina 2023-11-30
评论： Mathieu NOE 2023-12-8
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 个评论
Selina
Selina 2023-12-1
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.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Mathieu NOE
Mathieu NOE 2023-12-1
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 个评论
显示 6更早的评论
Selina
Selina 2023-12-5
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 2023-12-8
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])

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Gas Dynamics 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by