Прямое интегрирование уравнений движения

Для выполнения расчета с использованием прямого интегрирования уравнений движения следует выполнить следующие действия:

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

\[ \left\{ {{\begin{array}{*{20}c} {{\rm {\bf M\ddot{{u}}}}_{d} +{\rm {\bf C\dot{{u}}}}_{d} +{\rm {\bf K}}_{t} {\rm {\bf u}}_{d} ={\rm {\bf p}}_{d} \left( t \right)-{\rm {\bf K}}_{21}^{t} {\rm {\bf \bar{{u}}}}\left( t \right)-{\rm {\bf C}}_{21} {\rm {\bf \dot{{\bar{{u}}}}}}\left( t \right)} \\ {{\rm {\bf u}}_{d} \left( 0 \right)={\rm {\bf \dot{{u}}}}_{d} \left( 0 \right)=0} \\ \end{array} }} \right. \]

\({\bf u}_d\) – вектор узловых перемещений, обусловленных действием динамической нагрузки,

\({\bf M}\) – матрица масс,

\({\bf C}\) – матрица демпфирования,

\({\bf K}_t\) – матрица тангенциальной жесткости, определяемая по завершению статического расчета и остающаяся постоянной на всем этапе решения динамической задачи,

\({\bf p}_d(t)\) – вектор внешних сил динамического нагружения,

\( \bf \dot{\bar{u}}_{d} \) – зависящие от времени заданные перемещения,

\( \bf \bar{u}_{d} \) – зависящие от времени скорости, определяемые как производные по времени от заданных перемещений,

\( \bf{K}^{t}_{21} \) – подматрица тангенциальной жесткости \({\bf K}_t\), охватывающая степени свободы, по которым заданы вынужденные перемещения,

\( \bf{C}_{21} \) – подматрица диссипации C, охватывающая степени свободы, по которым заданы вынужденные перемещения.

Обоснование приведенного выше уравнения для линейной задачи дано в [2].

После решения динамической задачи полное решение определяется как суперпозиция статического решения и динамического:

\[ {\rm {\bf u}}\left( t \right)={\rm {\bf u}}_{stat} +{\rm {\bf u}}_{d} \left( t \right),\quad f\left( t \right)=f_{stat} +f_{d} \left( t \right) \]

где \( \bf u \) и \( \bf f \) – соответственно полный вектор перемещений и полный силовой фактор (внутреннее усилие, напряжение, деформация в любой точке), \( {\bf u}_{stat} \) и \( {\bf f}_{stat} \) – вектор перемещений и силовой фактор, обусловленные действием статической нагрузки, \( f_{d} \) – силовой фактор, обусловленный динамической нагрузкой, \( t \) – произвольный момент времени.

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

\[ \begin{array}{l} \left\{ {{\begin{array}{*{20}c} {{\rm {\bf M\ddot{{u}}}}+{\rm {\bf C\dot{{u}}}}+N\left( {{\rm {\bf u}}} \right)+{\rm {\bf \Gamma u}}={\rm {\bf p}}_{stat} +{\rm {\bf p}}_{d} \left( t \right)+{\rm {\bf \Gamma \bar{{u}}}}\left( t \right),} \\ {{\rm {\bf u}}\left( 0 \right)={\rm {\bf u}}_{0} ,\quad {\rm {\bf \dot{{u}}}}\left( 0 \right)=0,\quad } \\ \end{array} }} \right. \\ {\rm {\bf p}}_{d} \left( t \right)=\sum\limits_{p=1}^{N_{p} } {\mu_{p} {\rm {\bf k}}_{p}^{d} f_{p} \left( {t-t_{0}^{p} } \right)} ,\quad {\rm {\bf \bar{{u}}}}\left( t \right)=\sum\limits_{s=1}^{M_{s} } {\mu_{s} {\rm {\bf k}}_{s}^{d} f_{s} \left( {t-t_{0}^{s} } \right)} . \\ \end{array} \]

Здесь \( {\bf u} = {\bf u}_{t} \) – полный вектор перемещений, обусловленный действием как статической \(p_{stat}\), так и динамической \({\bf p}_d(t)\) нагрузки, \( {\bf u}_{0}\) – перемещения, полученные в результате решения нелинейной задачи статики \( N({\bf u}_0) = {\bf p}_{stat} \), \( N({\bf u}) \) – нелинейный оператор, возвращающий вектор внутренних сил системы. Для физически нелинейных задач тип этого оператора задается принятой теорией пластичности – деформационной или пластического течения. В случае деформационной теории пластичности оператор \( N({\bf u}) \) определяется диаграммой σ – ε для однородного материала, а для железобетона – диаграммами σ – ε для бетона и арматуры. В случае теории пластического течения оператор \( N({\bf u}) \) определяется видом заданной поверхности текучести – фон Мизеса для однородного материала и арматуры в случае железобетона и Друкера-Прагера или Гениева для бетона. Диаграммы σ – ε для однородного материала, бетона и арматуры определяют, как в процессе развития пластических деформаций поверхность текучести будет расширяться (изотропное упрочнение) или сжиматься (разупрочнение – типично для бетона при образовании и раскрытии трещин) и перемещаться в пространстве главных напряжений.

Статическая нагрузка \( {\bf p}_{stat} \) и неоднородные начальные условия \( {\bf u}_{0}\) передаются в динамическую часть задачи автоматически, поэтому пользователь не должен задавать статическую нагрузку в исходных данных для динамической задачи.

Динамическая нагрузка \({\bf p}_d(t)\) представляется в виде разложения по независящим от времени векторам пространственной конфигурации нагрузки \( \bf{k}_{d}(t) \), умноженным на функции времени \( f_{p}(t – t_0^p) \), где \( t_0^p) \) – время запаздывания. Здесь \( {\bf N}_p \) – количество ссылочных загружений, а \( p \) – текущий номер статического загружения при расчете на силовые воздействия, задающего пространственную конфигурацию нагрузки \( {\bf k}^d_p \). Коэффициенты \( μ_p \) играют роль масштабных множителей.

При наличии заданных перемещений используется метод параметров штрафа, в результате чего возникает диагональная матрица Γ. При отсутствии заданных смещений \( {\bf \Gamma } = 0 \). Заданные смещения представляются аналогично динамической нагрузке – в виде разложения по независящим от времени векторам пространственной конфигурации нагрузки \( {\bf k}^d_s \), умноженным на функции времени \( f_s(t – t_0^s) \), где \( t_0^s \) – время запаздывания, \( M_s \) – количество ссылочных загружений, а \( s \) – текущий номер статического загружения при расчете на заданные перемещения, задающего пространственную конфигурацию нагрузки \( {\bf k}^d_s \).

Для активации режима решения нелинейных уравнений движения следует использовать маркер Нелинейный анализ.

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

Кнопка позволяет вызвать диалоговое окно и осуществить тонкую настройку параметров нелинейного расчета:

Опция Игнорировать силы инерции и диссипации позволяет на основе алгоритма, разработанного для нелинейной динамики, решать задачи статики

\[ \left\{ {{\begin{array}{*{20}c} {N\left( {{\rm {\bf u}}} \right)+{\rm {\bf \Gamma u}}={\rm {\bf p}}_{stat} +{\rm {\bf p}}\left( t \right)+{\rm {\bf \Gamma \bar{{u}}}}\left( t \right)} \\ {{\rm {\bf u}}\left( 0 \right)={\rm {\bf u}}_{0} } \\ \end{array} }} \right. \]

При \( {\bf p}_{stat} = 0 \) в диалоге Управление шаговым процессом, вызываемом из узла дерева Моделирование нелинейных нагрузок, следует задать только один шаг с коэффициентом загружения равным нулю для инициализации физически нелинейных конечных элементов. При этом  \( {\bf u}_0 = 0 \). Такой подход позволяет намного гибче задавать приращение статической нагрузки. Например, можно реализовать такое нагружение: \( {\bf p}(t) = p_0\cdot sin(πt/4) \), где \( t \in [0, 1] \). Тогда при равномерном шаге  \( \Delta t \) по квазивременному параметру \(t\) приращение нагрузки \( \Delta{\bf p} \) → 0 при \(t\) → 1, что может оказаться полезным при приближении к предельной точке кривой равновесных состояний.

Данная опция позволяет также моделировать циклическое нагружение.

Максимальное количество итераций метода Ньютона-Рафсона ограничивает количество итераций, выполняемых упомянутым методом, во избежание зацикливания на данном шаге по времени при отсутствии сходимости итерационного процесса. Если сходимость с заданной точностью, определяемой параметром Критерий сходимости итераций метода Ньютона-Рафсона, достигнута, то итерации прекращаются и происходит переход к следующему шагу интегрирования. Алгоритм решения нелинейной динамической задачи представлен в [1].

Кнопка Факторизовать матрицу жесткости на каждой итерации для задач упруго-пластического деформирования рекомендуется всегда включать, поскольку в противном случае мало шансов обеспечить сходимость метода Ньютона-Рафсона. Отключение этой кнопки приводит к тому, что матрица тангенциальной жесткости вычисляется только один раз в начале каждого шага по времени или по нагрузке.

 

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

 

Литература

[1] Фиалко С.Ю. Применение метода конечных элементов к анализу прочности и несущей способности тонкостенных железобетонных конструкций с учетом физической нелинейности, СКАД СОФТ, Издательский дом АСВ, Москва, 2018.

[2] S. Fialko, V. Karpilovskyi, Time history analysis formulation in SCAD FEA software, J. of Measurements in Engineering, 6 (4), 2018, 173 – 180.