Finding neutral axis for arbitrary 2d shape - aerofoil
5 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Selina
am 30 Nov. 2023
Kommentiert: Mathieu NOE
am 8 Dez. 2023
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.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1557084/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1557089/image.png)
Is there a way to calculate it from this? I would appreciate any help!
2 Kommentare
Akzeptierte Antwort
Mathieu NOE
am 1 Dez. 2023
hello @Selina
the equations are still valid , see below how we can use them :
result
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1557659/image.png)
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 Kommentare
Mathieu NOE
am 8 Dez. 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])
Weitere Antworten (0)
Siehe auch
Kategorien
Mehr zu Numerical Integration and Differential Equations finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!