Finding neutral axis for arbitrary 2d shape - aerofoil

조회 수: 15 (최근 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
Dyuman Joshi
Dyuman Joshi 2023년 12월 1일
What is the definition of 'nuetral axis' in this context?
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
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 CenterFile Exchange에서 Programming에 대해 자세히 알아보기

Community Treasure Hunt

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

Start Hunting!

Translated by