Бенчмарк аналитикой SCAD++, Lira и ammonit3d. Тест на точность с одним конечным элементом

Предыстория

В предыдущей статье «Облако своими руками для расчета пространственных стержней методом конечных элементов на Node js, React js и Three js» представлен краткий обзор облачного SPA приложения ammonit3d по моделированию пространственных стержневых систем (ферм, балок, рамных и связевых конструкций, опор ЛЭП) методом конечных элементов с численно-аналитическим решением для каждого конечного элемента, в основе которого математическая модель Эйлера-Бернулли — механическая модель упругой балки или стержня длиной L с заданной изгибной жёсткостью EJ на которую действуют сосредоточенная сила F или момент M, а также распределённая сила q(x) или момент m(x) по длине стержня, продольная ось которого x1, вертикальная x2 проходит через начало стержня и x3 направлена на нас из точки пересечения x1, x2. Полагая сечения плоскими до и после изгиба при одноосном деформированном состоянии в рамках краевой задачи линейной теории упругости, уравнение упругой оси стержня можно представить в виде обыкновенного неоднородного дифференциального уравнения 4-го порядка:

$EJ_3 \dfrac{d^4v}{dx^4} = q(x) - \dfrac{dm}{dx}, \quad q(x) = kx + q_0, \quad m(x) = rx + m_0,\quad x\in[0,L],$

общее решение которого состоит из частного решения правой части (первые два слагаемых) и решения однородного дифференциального уравнения:

$v(x)=\dfrac{kx^5}{120EJ_3}+\dfrac{(q_0-r)x^4}{24EJ_3}+a_4x^3+a_3x^2+a_2x+a_1.$

Если записать граничные условия

$\begin{equation}\begin{cases}\begin{array} \, v(0)=v^1,\quad v,_1(0)=\varphi^1,\\v(L)=v^2, \quad v,_1(L)= \varphi^2,\end{array}\end{cases}\end{equation}$

то неизвестные коэффициенты можно выразить через свёртку (сумму произведений) прогибов и поворотов на концах стержня с эрмитовыми полиномами N:

$\begin{array} \, v(x) = \dfrac{k}{EJ_3}\left(\dfrac{x^5}{120}-\dfrac{x^3L^2}{40}+\dfrac{x^2L^2}{60}\right)+\dfrac{q_0-r}{EJ_3}\left(\dfrac{x^4}{24}-\dfrac{x^3L}{12}+\dfrac{x^2L^2}{24}\right)+N_p^i\alpha_i^p,\\N_p^1\alpha_1^p=\left(\dfrac{2x^3}{L^3}-\dfrac{3x^2}{L^2}+1\right)\,v^1 + \left(-\dfrac{2x^3}{L^3}+\dfrac{3x^2}{L^2}\right)\,v^2,\quad i=1,2, \quad p = 1,2,\\ N_p^2\alpha_2^p=\left(\dfrac{x^3}{L^2}-\dfrac{2x^2}{L}+x\right)\,\varphi^1 + \left(\dfrac{x^3}{L^2}-\dfrac{x^2}{L}\right)\,\varphi^2,\quad \alpha_1^p=v^p, \quad \alpha_2^p=\varphi^p. \end{array}$

Сравнение результатов SCAD++, Lira и ammonit3d на одном конечном элементе с аналитикой

Рис.1 Модель балки длиной L=1 м, с жёстким защемлением концов и сечением 25х25мм из стали E=2.06e+11, v=0.3.

Рассмотрим модель балки длиной L=1 м с жёстким защемлением концов. Сечение балки — квадрат 25×25 мм из стали с модулем Юнга E = 2.06e+11 Па и коэффициентом Пуассона 0.3. По всей длине на балку сверху действует статическая равномерно распределённая нагрузка 3 кН/м, см. Рис 1. В этом случае из общего решения обыкновенного неоднородного дифференциального уравнения 4-го порядка (см. выше) можно выписать аналитические выражения для функции прогиба v(x) и изгибающего момента M(x):

$v(x) = \dfrac{q_0}{EJ}\left(\dfrac{x^4}{24}-\dfrac{x^3L}{12}+\dfrac{x^2L^2}{24}\right), \quad M(x) = -EJ\dfrac{d^2v}{dx^2} = -q_0\left(\dfrac{x^2}{2}-\dfrac{xL}{2}+\dfrac{L^2}{12}\right).$

Используя эти соотношения можно построить графики функций прогиба v(x) и момента M(x) по длине балки, см. Рис. 2.

Рис. 2. Зависимости функции прогиба v(x) и момента M(x) по длине балки.

Из аналитических выражений и графиков, представленных на Рис. 2, видно, что наибольший прогиб v(0.5) = -1.165e-3 м, достигается в середине балки. Наибольший изгибающий момент достигается на концах балки М(0) = М(1) = -250Нм. В середине балки момент равен М(0.5) = 125Нм. На Рис. 3, 4 представлены графики прогиба и момента по длине балки, найденный в программе SCAD++, Lira и ammonit3d, соответственно.

Рис. 3 Графики прогиба v(x) и момента M(x) по длине балки, найденный в программе SCAD++

Рис. 4 Графики момента M(x) и прогиба v(x) по длине балки, найденный в программе Lira

Рис. 5 Графики момента M(x) и прогиба v(x) по длине балки, найденный в ammonit3d, модель: analytic

Из таблицы 1 видно, что параметры, найденные аналитически, совпадают с численным решением модели из одного конечного элемента для SCAD++ и Lira, так и для ammonit3d. Однако SCAD++, Lira в рамках одного конечного элемента не отображает эпюры прогибов на конечно-элементной расчётной схеме(?), а лишь в отдельном окне эпюр. При этом ammonit3d обеспечивает отображение и прогибов и моментов на пространственной схеме.

Что сравниваем	$v(0.0), 10^{-3} м$	$M(0.5), \,Н\cdotм$	$M(0.0), \,Н\cdotм$
Аналитическая формула	-1.165000	125	250
SCAD++	-1.164999	125	250
Lira	-1.165000	125	250
ammonit3d	-1.165000	125	250

Вывод

Безусловно, коммерческие пакеты SCAD++ и Lira обеспечивают в рамках модели Эйлера-Бернулли аналитическую точность на одном конечном элементе, как и ammonit3d. Однако облачное приложение демонстрирует современную графику и кроссплатформенность, т. е. возможность создавать модели и выполнять расчёты на компьютерах, планшетах, смартфонах и других гаджетах, имеющих доступ к интернету. И главная особенность – моделями 3Д расчётных схем можно делиться в соцсетях, по почте, размещать ссылки публичных моделей на различных сайтах или делиться приватными моделями в своём сообществе.

ссылка на оригинал статьи https://habr.com/ru/articles/1024876/