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

Для начала нужно обзавестись этим самым «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 r\mathit{dV}} V$

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

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

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

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

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

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

$J_i=m_i\left[\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^2\end{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 _i\overrightarrow r_iS_i=\frac{\sum _i\overrightarrow 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 n\right)\mathit{dS}=\iiint div \left(\overrightarrow F\right)\mathit{dV}$

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

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

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

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

$\overrightarrow F=\frac 1 3 \left[x,y,z\right]$

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

$div \left(\overrightarrow F\right)=\frac{∂F_x}{∂x}+\frac{∂F_y}{∂y}+\frac{∂F_z}{∂z}=1$

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

$\iiint \mathit{dV}=V=\iint \left(\overrightarrow F,\overrightarrow n\right)\mathit{dS}≈\sum _i\left(\overrightarrow F_i,\overrightarrow n_i\right)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 r\mathit{dV}} V$

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

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

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

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

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

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

$x_{\mathit{цм}}=\frac 1 V\iiint \mathit{xdV}=\frac 1 V\iiint div \left(\overrightarrow F^x\right)\mathit{dV}=\frac 1 V\iint \left(\overrightarrow F^x,\overrightarrow n\right)\mathit{dS}≈\frac 1 V\sum _i\left(\overrightarrow F^x_i,\overrightarrow n_i\right)S_i$

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

$\overrightarrow a\ast \overrightarrow b=\left[a_xb_x,a_yb_y,a_zb_z\right]$ $\overrightarrow r_{\mathit{цм}}=\frac 1 V\sum _i\left(\overrightarrow F_i\ast \overrightarrow n_i\right)S_i$

где

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

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

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^2\right)\mathit{dV}=ρ\iiint div \left(\overrightarrow F\right)\mathit{dV}=ρ\iint \left(\overrightarrow F,\overrightarrow n\right)\mathit{dS}=ρ\sum _i\left(\overrightarrow F_i,\overrightarrow n_i\right)S_i$

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

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

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

$div \left(\overrightarrow F\right)=-xy$

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

$\overrightarrow F=\left[0,0,-xyz\right]$

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

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% совпадет с реальностью…

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

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