Объем, центр масс, моменты инерции тела имея только mesh поверхности

в 20:30, , рубрики: CAD/CAM, python, масса, математика, Объем, тензор инерции, центр масс
Объем, центр масс, моменты инерции тела имея только mesh поверхности - 1

Для начала нужно обзавестись этим самым "mesh"-ем поверхности, или триангуляцией поверхности, полигональной сеткой, разбиением двумерного многообразия. В данном случае работа будет вестись именно с треугольной сеткой, но все ниже представленные формулы и код (если немного модифицировать), будет работать с сеткой состоящей из любых полигонов. Главное, чтобы они были малые, от этого зависит точность, чем меньше - тем лучше.

Код здесь написан языке Python, mesh импортируется из .stl файла. Откуда он будет импортироваться совершенно неважно, Python дальше будет работать с данными, которые имеют тип массивов numpy.

данные выглядят так
# Треугольники (список треугольников, в каждый из которых входят 3 точки)
[[[-5.0500002e+00  0.0000000e+00  0.0000000e+00]
  [-4.9604182e+00  2.5921434e-02 -3.2746387e-03]
  [-4.9604182e+00  2.5664667e-02 -4.8957970e-03]]

 [[-4.9604182e+00  2.6075900e-02 -1.6405566e-03]
  [-4.9604182e+00  2.5921434e-02 -3.2746387e-03]
  [-5.0500002e+00  0.0000000e+00  0.0000000e+00]]

 [[-4.9604182e+00  2.6127456e-02  0.0000000e+00]
  [-4.9604182e+00  2.6075900e-02 -1.6405566e-03]
  [-5.0500002e+00  0.0000000e+00  0.0000000e+00]]

 ...

 [[-4.8888855e+00  4.6034887e-02 -2.8962698e-03]
  [-4.9604182e+00  2.5921434e-02 -3.2746387e-03]
  [-4.9604182e+00  2.6075900e-02 -1.6405566e-03]]

 [[-4.8888855e+00  4.6034887e-02 -2.8962698e-03]
  [-4.9604182e+00  2.6075900e-02 -1.6405566e-03]
  [-4.8888855e+00  4.6125907e-02  0.0000000e+00]]

 [[-4.8888855e+00  4.6125907e-02  0.0000000e+00]
  [-4.9604182e+00  2.6075900e-02 -1.6405566e-03]
  [-4.9604182e+00  2.6127456e-02  0.0000000e+00]]]
# Нормали к треугольникам (вектора)
[[ 0.27986652 -0.94821906  0.15018432]
 [ 0.27986658 -0.95577824  0.09034719]
 [ 0.27986655 -0.95956516  0.03015531]
 ...
 [ 0.26912326 -0.95883155  0.09063575]
 [ 0.26912323 -0.96263045  0.03025224]
 [ 0.26912326 -0.96263045  0.03025234]]
И здесь это дело импортируется из STL файла, используя библиотеку numpy-stl
# pip install numpy-stl  # чтобы установить библиотеку
my_mesh = mesh.Mesh.from_file('название_файла.STL')
Normals = my_mesh.normals 
Triangles = my_mesh.vectors

# и сразу нормализуем нормали, на всякий пожарный
Nnorm = np.linalg.norm(Normals, axis=1) 
Nnorm = np.tile(Nnorm.reshape(len(Nnorm), 1), (1, 3))
Normals = Normals/Nnorm

Теперь у нас есть треугольники и единичные нормали к ним. Но для дальнейшей работы еще понадобится список площадей треугольников и их центров масс.

Получаем площади треугольников и центры масс треугольников
A, B, C = Triangles[:, 0], Triangles[:, 1], Triangles[:, 2] 
R1, R2 = C-A, B-A

# список центров масс
R = 1/3*(A+B+C)

# список площадей 
S = np.cross(R1, R2)
S = np.linalg.norm(S, axis=1)/2

Теперь у нас есть:

  • S - список площадей треугольников

  • R - список центров масс треугольников

  • Normals - список нормалей к треугольникам

  • Triangles - список треугольников (треугольник содержит 3 точки)

Теперь следует напомнить как вычисляются (чисто математически) массы, объемы, площади, центры масс, моменты инерции и что вообще это такое и зачем вообще нужно.

А те кто в теме и сомневаются, что все эти вещи можно вычислять для объемных тел используя только 2d сетку, скажу: это делать можно. И нужно, если лень разбивать 3d область пространства на маленькие элементы - тетраэдры там, призмы... Да, это делать возможно, и мы это сделаем. Но всё по порядку.

Масса - вычисляется как интеграл по объему (площади, линии):

m=iiint mathit{rho dV}

здесь под интегралом (тройным) - плотность вещества. Объемная плотность вещества, имеющая размерность кг/м^3. Но, если в наличии имеется поверхностная плотность (кг/м^2), или линейная (?? кг/м), то интегралы соответственно будут двойные и простые.

Объяснять зачем нужна масса мало кому нужно, вернее об это знают не только лишь все. Масса на ускорение - это сила. Второй закон Ньютона, известный всем еще с седьмого класса.

Для большинства приложений в инженерии достаточно работать с плотностями, которые постоянны в пространстве и во времени =) Потому мы упростим себе жизнь и вынесем за все интегралы с которыми будем работать все плотности. И даже работать мы будем не с массами, а с площадями и объемами. А когда будем вычислить центры масс, плотности вообще сами по себе сократятся, и неважны.

Кстати, центр масс:

overrightarrow r_{mathit{цм}}=frac{iiint overrightarrow rmathit{dV}} V

Ну или момент первого порядка, или как там его называют в статистике.

В инженерном деле (конструкторском) центр масс важен, например, для описания движения твердых (недеформируемых) тел. Конструктора помещают туда связанную с телом координатную систему (как-то располагая её относительно тела, под каким-то углом в смысле) и с этих пор тело и связанная кс неразлучны друг с другом и неподвижны друг относительно друга.

Зачем помещают связанную систему в центр масс? Ну так принято, пойдем дальше.

Вот мы пришли к моментам инерции. Момент инерции - это то, что влияет на выбор конструктора, когда перед ним встает задача: как именно закрепить на теле связанную координатную систему.

Удачный выбор закрепления будет влиять вполне конкретно на уравнения движения этого самого тела. Момент инерции - это аналог массы во втором законе Ньютона для точки. Еще говорят "тензор инерции" - он участвует в описании вращательного движения тела, но это далеко не скалярная величина, это не просто число. Тензор инерции имеет валентность матрицы, то есть таблицы порядка 2 на 2 или 3 на 3, в зависимости от размерности пространства.

Ну а закрепляют кс исходя из симметрии тела (корабли, самолеты, дирижабли.. ) имеют зачастую симметрию относительно вертикальной плоскости. Вот поэтому одна из координатных плоскостей кс совпадает с ней. А дальше у нас остается только одна степень свободы как расположить - угол. И обычно одну из осей кс, находящихся в вертикальной плоскости, направляют горизонтально, параллельно земле. Тогда тензор инерции будет иметь почти идеальный диагональный вид, и, кстати, его можно было бы сделать диагональным в любом случае, для любого тела, но всё же кс располагают так как я описал. Ну чтобы диагональные элементы имели более ясный, практичный физический смысл. Что-то далеко я отошел от темы, вычисляется момент инерции так:

J_i=m_ileft[begin{matrix}y_i^2+z_i^2&-x_iy_i&-x_iz_i\…&x_i^2+z_i^2&-y_iz_i\…&…&x_i^2+y_i^2end{matrix}right]

это момент инерции для точки, для тела бы записалось как сумма по всем его точкам, или интеграл. А точнее, куча интегралов, в точности 6 различных штук, ведь тензор инерции - симметричная матрица (и еще её всегда можно диагонализировать - повернуть в пространстве координатную систему, относительно которой и вычисляется момент, так, чтобы матрица приняла диагональный вид - все элементы 0, которые не на диагонали).

А еще тензор инерции можно пересчитывать в другие координатные системы смещенные относительно центра масс по теореме Гюйгенса. То есть, есть формула, связывающая моменты инерции в разных координатных системах, и зная в одном можно найти в другом имея только лишь вектор смещения. Но важно! теорема Гюйгенса работает только когда одна из кс расположена в центре масс. Если есть две кс вне центра масс, в формуле будут появляться дополнительные члены (моменты первого порядка).

не нужно это читать

если статья наберет 1000 лайков, в следующих статьях докажу теорему Гюйгенса и повороты тензора инерции в элегантных тензорных обозначениях, которые полюбил Эйншейн..

Вообще и дизлайками не брезгую..

Итак, достаточно. Всё что нужно напомнилось (было напомнено), короче, вспомнено. Восполнено, ладно, клоунада как клоунада, она иногда надо.

А давайте посчитаем все эти вещи просто для поверхности, чтобы набить немного руку?

Масса поверхности получается как сумма площадей всех треугольников умноженная на поверхностную плотность:

m_{mathit{об}}=ρ_{mathit{пленки}} cdot sum _iS_i

def get_m_obolochka(S, rho):
    m = np.sum(S) * rho
    return m

Центр масс поверхности (оболочки):

overrightarrow r_{mathit{цм}mathit{оболочки}}=frac 1{m_{mathit{об}}}iint overrightarrow rρ_{mathit{пленки}}mathit{dS}≈frac{ρ_{mathit{пленки}}}{m_{mathit{об}}}sum _ioverrightarrow r_iS_i=frac{sum _ioverrightarrow r_iS_i}{sum _iS_i}

def get_cm_obolochka(S, R):
    S_sum = np.sum(S)
    cm = np.einsum('ij,i', R, S) / S_sum
    return cm

Момент инерции оболочки:

J_{mathit{об}}=ρ_{mathit{пленки}}sum _iJ_i cdot S_i

где матрица J_iвычисляется по приведенной где-то выше формуле, а еще здесь она не умноженная на i-тую массу.

def get_J_for_ob(rho, R, S):
	x, y, z = R[:, 0], R[:, 1], R[:, 2]
	x2, y2, z2 = x**2, y**2, z**2

	Jx, Jy, Jz = y2+z2, x2+z2, x2+y2
	Jxy, Jyz, Jxz = -x*y, -y*z, -x*z

	J = np.array([	[Jx, Jxy, Jxz],
									[Jxy, Jy, Jyz], 
									[Jxz, Jyz, Jz]
							])
	J = J.T # это нужно, потому что если подумать, то нам нужен список J[i] для всех треугольников
	# короче делаем выше правильный shape =)
  J = np.einsum('ijk,i->jk', J, S)*rho
	return J

Теперь, когда размялись, можно приступать к основному, тому, что написано в названии статьи. Вычислим все эти характеристики для объемного тела, которое ограничено mesh-ем, используя только mesh.

Для этого понадобится, внимание... теорема Остроградского-Гаусса:

iint left(overrightarrow F,overrightarrow nright)mathit{dS}=iiint div left(overrightarrow Fright)mathit{dV}

На мой взгляд это самая красивая теорема математики... на самом деле она вроде как есть частный случай операций над дифференциальными формами, но я не слишком математик, поверхностно знаком ;) Да... там и теорема Гаусса-Бонне и все в этом духе, это топология, детка. Это всё очень красиво, и может быть я этим когда-нибудь займусь, когда меня освободят из рабства.

Формула позволяет переходить от вычислений по объему, к вычислениям по поверхности, ограничивающей заданный объем. Это можно интерпретировать как проекцию. Мы первое что видим, когда рождаемся - это проекцию трехмерной мамы на нашу двухмерную сетчатку глаза....

Да и куда ни глянь, всё есть проекцией, всё является лишь тенью более сложной "настоящей" вещи, вещи из "реального" мира.. а она в свою очередь - еще одна проекция из следующей итерации, следующего под-уровня реальности, под-Матрицы Вачовской.

Пренебрегая философией, и подбирая векторное поле следующим образом:

overrightarrow F=frac 1 3  left[x,y,zright]

мы видим, что дивергенция становится равна:

div left(overrightarrow Fright)=frac{∂F_x}{∂x}+frac{∂F_y}{∂y}+frac{∂F_z}{∂z}=1

Это означает, что теперь зная лишь mesh поверхности, мы можем вычислить объем:

iiint mathit{dV}=V=iint left(overrightarrow F,overrightarrow nright)mathit{dS}≈sum _ileft(overrightarrow F_i,overrightarrow n_iright)S_i

Можно сказать, что здесь:

overrightarrow F_i=frac 1 3 overrightarrow R_i

(Центры масс треугольников, напоминаю)

Код:

def Volume(S, R, Normals):
	V = np.einsum('ij,ij,i', R, Normals, S)*1/3 # einsum как работает можно найти в документации
	return V

Но можно обнаглеть, и пойти дальше объема:

overrightarrow r_{mathit{цм}}=frac{iiint overrightarrow rmathit{dV}} V

Используя ту же теорему и задав вектор F так (для вычисления координаты x):

overrightarrow F^x=frac 1 2 left[x^2,0,0right]

(чего-то индекс x не туда залез, я его хотел чутка ниже поместить)

Его дивергенция:

div left(overrightarrow F^xright)=frac 1 2left(frac{∂x^2}{∂x}+0+0right)=x

Координата x центра масс тела:

x_{mathit{цм}}=frac 1 Viiint mathit{xdV}=frac 1 Viiint div left(overrightarrow F^xright)mathit{dV}=frac 1 Viint left(overrightarrow F^x,overrightarrow nright)mathit{dS}≈frac 1 Vsum _ileft(overrightarrow F^x_i,overrightarrow n_iright)S_i

Аналогичным образом вычисляются остальные две координаты. Компактно это можно записать если ввести операцию покоординатного умножения (использующегося в numpy):

overrightarrow aast overrightarrow b=left[a_xb_x,a_yb_y,a_zb_zright]overrightarrow r_{mathit{цм}}=frac 1 Vsum _ileft(overrightarrow F_iast overrightarrow n_iright)S_i

где

overrightarrow F=frac 1 2left[x^2,y^2,z^2right]

Воплощение в коде:

def get_cm_V(S, R, Normals, V_t):
	F = R**2/2
	cm = np.einsum('ij,ij,i->j', F, Normals, S)
	return cm/V_t

Для объемного тела подсчитать тензор инерции можно используя только триангуляцию оболочки, аналогично как делалось для подсчета объема и центра масс:

J_{mathit{xx}}=ρiiint left(y^2+z^2right)mathit{dV}=ρiiint div left(overrightarrow Fright)mathit{dV}=ρiint left(overrightarrow F,overrightarrow nright)mathit{dS}=ρsum _ileft(overrightarrow F_i,overrightarrow n_iright)S_i

Где вектор F можно, например, задать так:

overrightarrow F=left[xleft(y^2+z^2right),0,0right]

Для компоненты J_{xy} должно выполняться:

div left(overrightarrow Fright)=-xy

Соответственно можно задать F:

overrightarrow F=left[0,0,-xyzright]

Аналогично и с остальными компонентами.

def get_J_for_V(S, R, Normals, rho):
	x, y, z = R[:, 0], R[:, 1], R[:, 2]
	x2, y2, z2 = x**2, y**2, z**2

	Fx, Fy, Fz = x*(y2+z2), y*(x2+z2), z*(x2+y2)
	Fxyz = -x*y*z

	Fx, Fy, Fz = Fx[:, None], Fy[:, None], Fz[:, None]
	Fxyz = Fxyz[:, None]

	N = len(x)
	zero = np.zeros((N, 1))

	Fxx, Fyy, Fzz = np.hstack((Fx, zero, zero)), np.hstack((zero, Fy, zero)), np.hstack((zero, zero, Fz))
	Fxy, Fxz, Fyz = np.hstack((zero, zero, Fxyz)), np.hstack((zero, Fxyz, zero)), np.hstack((Fxyz, zero, zero))

	Jx, Jy, Jz = np.einsum('ij,ij,i', Fxx, Normals, S), np.einsum('ij,ij,i', Fyy, Normals, S), np.einsum('ij,ij,i', Fzz, Normals, S)
	Jxy, Jyz, Jxz = np.einsum('ij,ij,i', Fxy, Normals, S), np.einsum('ij,ij,i', Fyz, Normals, S), np.einsum('ij,ij,i', Fxz, Normals, S)

	J = np.array([	[Jx, Jxy, Jxz],
					[Jxy, Jy, Jyz], 
					[Jxz, Jyz, Jz]
				])*rho

	return J
# не так элегантно, но работает вроде.. 

На этом всё, спасибо за внимание.

еще

Если кому не нравится мой стиль повествования и раздражает или обижает, не обижайтесь, пожалуйста, я тоже часто грущу на самом деле... грущу чаще, чем смеюсь и радуюсь, это нормально для таких как мы... Просто, когда сажусь писать на Хабр, чет пробивает на шутки какие-то дурацкие и клоунство. Ну как написал, так написал. Сейчас выложу.
И еще. Кому интересно, я занимаюсь дирижаблестроением, пока что маленьким, но вижу как в будущем из этого зерна прорастает сперва зелень, потом колос, потом полное зерно в колосе. И кому это интересно, пусть связывается со мной, может будем сеять вместе, и вместе пожинать плоды.. Можно тут в личных сообщениях, также я есть на фейсбуке, в телеграме.

А может кто-то хочет помочь материально, тоже пусть свяжется. Так дело пойдет намного быстрее, потому что сейчас я не живу, а выживаю, такая правда.

И еще, скептикам, которые не верят, что маленькие дирижабли способны противостоять ветру:

да, это всего лишь симуляция, но я уверен, что она на 80-90 % совпадет с реальностью...

позже сделаю ветер порывистым, и случайным образом генерирующимся, как по величине скорости, так и по направлению, чтобы было более похоже на реальность...

Подписывайтесь на мой ютуб-канал, это мой видео-дневник дирижаблестроителя

Автор: Сергей

Источник


* - обязательные к заполнению поля


https://ajax.googleapis.com/ajax/libs/jquery/3.4.1/jquery.min.js