Введение

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

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

В настоящее время известно множество критериев разрушения, однако лишь некоторые из них получили широкое практическое применение [1–6]. Среди наиболее признанных можно выделить следующие: коэффициент интенсивности напряжений (K), инвариантный энергетический интеграл (J-интеграл) и параметр раскрытия трещины в ее вершине (δ).

В данной работе проведен расчет L-образной панели,изготовленной из квазихрупкого материала и закрепленной в шарнирном четырехзвеннике, при действии вертикальной растягивающей нагрузки F. Использовали метод конечных элементов (КЭ) [7–13] при определении предельных напряженных состояний в дискретизированной модели конструкции.

Методика оценки трещиностойкости конструкции

Методика, описанная в работе [14], предназначена для оценки прочности и трещиностойкости конструкций на стадии проектирования, что позволяет значительно снизить затраты, связанные с проведением натурных испытаний [15–23].

На начальном этапе оценки прочности конструкции используется теория формоизменения (хорошо согласуется с квазихрупкими материалами [14]), которая позволяет определить, выдержит ли конструкция приложенные нагрузки. Критерий прочности записывается в следующем виде:

                               (1)

где σ₁, σ₂, σ₃ – главные напряжения, МПа.

Трещина зарождается, когда удовлетворяется следующее условие [14]:

σ > σ_в,                                                                   (2)

где σ_в – предел прочности при растяжении, МПа.

В области вершины трещины возникают концентрации напряжений, достигающие критических значений, значительно превышающих предел прочности материала. Для оценки поведения материала в таких условиях и прогнозирования роста трещины применяются методы линейной упругой механики разрушения. В работе [14] приведены следующие параметры процесса:

– критическое значение интенсивности выделения упругой энергии с учетом трех типов трещин (G_c);

– комбинированный параметр (критерий Гриффитса и энергия формоизменения) с учетом трех типов трещин (U_c).

Методы линейной упругой механики разрушения применимы при соблюдении следующих ограничений [14]:

– размер зоны пластичности в окрестности вершины трещины составляет ≤20 % от ее длины, что соответствует условию маломасштабной текучести;

– относительное удлинение при разрушении должно составлять ≤20 %;

– разрушающее напряжение в нетто-сечении образца должно быть ≤80 % от предела текучести материала.

Далее рассчитываются критические значения условных коэффициентов интенсивности напряжений для всех трех мод разрушения [14]:

где σ_t – предел текучести при растяжении, МПа; τ_t – предел текучести при сдвиге, МПа; A – площадь плоскости α, заданная единичными векторами  χ¯ и γ‾  (рис. 1), мм².

Рис. 1. Плоскость α в произвольной трещине: σ₀ – номинальное напряжение [14]

На следующем этапе рассчитывается критическое значение интенсивности выделения упругой энергии G_с, которое учитывает влияние всех трех типов трещин. Для его определения используется следующая формула [14]:

                                        (6)

где μ – коэффициент Пуассона; Е – модуль упругости, МПа.

Энергия деформации в вершине трещины определяется посредством комбинированного применения теорий Гриффитса и энергии формоизменения [14]:

                     (7)

где V_КЭ – объем КЭ, мм³; A_КЭ – площадь геометрической фигуры, полученной пересечением ребер КЭ с плоскостью α, мм².

Множитель π вводится для корректировки неоднородности механических свойств материала и неравномерности распределения напряжений в окрестности трещины [14].

На рис. 2 представлена геометрическая фигура, образованная пересечением плоскости α с ребрами КЭ (в качестве примера взят четырехузловой тетраэдрический элемент). Площадь этой фигуры определяется как площадь треугольника ABC и обозначается как A_КЭ = S_ABC [14].

Рис. 2. Геометрическая фигура в форме треугольника ABC, полученная пересечением плоскости α с ребрами конечного элемента:  –векторы главных напряжений [14]

В данном случае ориентация плоскости α определяется двумя векторами главных напряжений (σ¯₂ и σ¯₃ ), действующими в узловой точке вершины трещины. Сама трещина, проходящая через четырехузловой тетраэдрический элемент, целиком лежит в этой плоскости α. При выполнении условия

U_c/G_c > 1,                                                               (8)

КЭ ослабляется на величину 10^–4 [14]. Это свидетельствует о начале распространения трещины, после чего соответствующие КЭ исключаются из расчета за счет уменьшения их локальной матрицы жесткости в 10^–4 раза.

Процесс анализа трещиностойкости продолжается до тех пор, пока выполняются следующие условия [14]:

– суммарные перемещения рассматриваемой итерации не превышают десятикратные суммарные перемещения предыдущей итерации

U_i ≤ 10U_i–1;                                                              (9)

– зарождение или распространение трещины соответствует формулам (2) и (8).

Ослабленные КЭ конструкции соответствуют предполагаемому пути распространения трещины, определяемому по предельным напряженным состояниям этих элементов.

В представленной работе для анализа зарождения и роста трещины используется комплексный подход, сочетающий силовой (формула (2)) и энергетический (формула (8)) критерии. Силовой и энергетический критерии дополняют друг друга, обеспечивая двухуровневую проверку на разрушение.

Проведение расчета трещиностойкости L-образной панели

Исходные данные L-образной панели

Для моделирования роста трещины выбрана L-образная панель толщиной 1 мм из квазихрупкого материала, которую подвергали воздействию вертикального усилия F (рис. 3). Параметры материала L-образной панели:

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

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

Рис. 3. L-образная панель

Разбиение на конечные элементы

Необходимо разбить L-образную панель на КЭ, т. е. дискретизировать рассматриваемую область Ω на некоторое количество подобластей Ω_е. В дальнейшем индекс e будет указывать на принадлежность той или иной величины конкретному КЭ. При этом КЭ не должны пересекаться между собой, т. е. не допускается их взаимное наложение, а также они должны полностью (без зазоров и разрывов) заполнять исходную область Ω.

Панель разбили на четырехузловые плоские изопараметрические КЭ размером 10×10 мм в локальной системе координат с осями ξ и η (рис. 4). Оси ξ и η локальной системы координат совпадали с осями x и y глобальной системы координат. Значения координат КЭ на осях ξ и η находились в диапазоне [–1; +1] (в безразмерных величинах). Это означает, что каждая координата изменялась в пределах от –1 на одной стороне элемента до +1 на противоположной стороне, принимая нулевое значение на медианах [10].

Выбор данного типа КЭ (четырехузловой изопараметрический элемент с билинейными функциями формы) обусловлен его эффективностью и устойчивостью при моделировании сингулярностей полей напряжений в окрестности вершины трещины. Линейные функции формы обеспечивают сходимость решения и адекватную аппроксимацию больших градиентов напряжений при умеренных вычислительных затратах. В качестве альтернативы можно рассмотреть восьмиузловые квадратичные элементы. Однако восьмиузловые элементы могут приводить к избыточной сложности сетки в зонах с высокими градиентами и повышению вычислительных затрат.

Рис. 4. Четырехузловой плоскодеформируемый элемент

В узлах КЭ были связаны между собой и взаимодействовали друг с другом как единое целое. Таким образом, КЭ взаимодействовали через «связи», установленные исключительно в узлах принятой конечно-элементной модели. В этих узловых точках прикладывали узловые силы и возникали узловые перемещения. Узловые перемещения представляют собой степени свободы узла – параметры, определяющие его положение в пространстве, а значит, положение всей L-образной панели в целом. Узлы располагались в вершинах КЭ (на рис. 4 обозначены черными точками). Вектор узловых перемещений элемента имеет вид

                                                 (10)

После разбиения сетка состояла из 1976 узлов и 1875 КЭ. Каждый элемент состоял из четырех узлов со своей глобальной нумерацией (рис. 5).

Рис. 5. Глобальная нумерация узлов

При построении конечно-элементной модели важным этапом является выбор базисных кусочно-определенных функций N₁–N₄ для подобластей Ω_е. Эти функции, известные как функции формы, предназначены для однозначного определения искомой функции в любой точке внутри элемента на основе дискретных значений этой функции, заданных в узловых точках. С их помощью описываются перемещения, деформации и напряжения в узлах КЭ через соответствующие дискретные параметры. Функции формы выражаются следующими соотношениями, записанными через локальные координаты [10]:

                               (11)

Метод построения локальной и глобальной матриц жесткости

Суть изопараметрической техники построения КЭ заключается в том, что для описания геометрических размеров КЭ используются те же функции формы, что и для описания перемещений узлов КЭ. Интерполирующие функции координат узлов через функции формы имеют вид:

                                                              (12)

где N_i(ξ, η) – функция формы i-го узла элемента; x_i, y_i – координаты i-го узла элемента.

Интерполирующие функции перемещений через функции формы, с помощью которых описывается потенциальное распределение перемещений по области КЭ вдоль осей координат:

                                                             (13)

где u_i, v_i – перемещение i-го узла элемента.

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

,                                           (14)

где  – матрица Якоби;  – деформации в глобальной системе координат;  – деформации в локальной системе координат.

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

Использованный изопараметрический четырехузловой КЭ должен быть применим для решения задач плоского напряженного состояния. В этом случае можно воспользоваться соотношениями между напряжениями и деформациями по закону Гука:

,                             (15)

где   – модуль сдвига.

Применено соотношение Коши для деформаций:

.                                                                (16)

Согласно закону Гука (15) выражена матрица упругих постоянных:

.                                                                   (17)

Из соотношений Коши (16) получена матрица дифференциальных операторов:

                                                                         (18)

Таким образом, можно записать следующее векторное уравнение:

{L_e}= [N]{u_e},                                                                          (19)

где {L_e} – вектор перемещений в текущей точке элемента;  – матрица базисных функций формы.

Деформация может быть выражена через узловые перемещения элемента:

{ε_e}= [C]{L_e}= [C][N]{u_e} = [B]{u_e},                                                    (20)

где [B] = [C][N] – матрица связи узловых деформаций и перемещений.

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

{σ_e}= [D]{ε_e}= [D][B]{u_e}.                                                           (21)

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

Рассмотрен вариационный подход получения матрицы жесткости КЭ, основанного на минимизации полной потенциальной энергии системы П_е [10, 14], при условии единичной толщины КЭ:

П_е = U_e + V_e,                                                                              (22)

где U_e – потенциальная энергия упругой деформации КЭ; V_e– работа внешних сил КЭ.

Потенциальная энергия упругой деформации КЭ определяется по формуле (интегрирование проводится по площади КЭ A)

              (23)

Выражение для работы внешних сил КЭ будет выглядеть следующим образом (интегрирование проводится по контуру поверхности Г в случае поверхностных сил и по площади КЭ А в случае объемных сил):

                   (24)

где {р} – вектор поверхностных сил; {f_e} – вектор объемных сил.

Следует отметить, что в выражениях (23) и (24) варьируемыми параметрами являются узловые перемещения {u_e}.

Вариация полной потенциальной энергии системы в рамках вариационного принципа по {u_e} с последующим приравниванием полученного выражения к нулю    позволяет определить вариацию потенциала Лагранжа:

           (25)

Следовательно, уравнение равновесия КЭ примет вид

[K_e]{u_e} = {F_e},                                                                                                           (26)

где  – локальная матрица жесткости КЭ;                                                                 (27)

 – вектор внешних сил, действующих на КЭ.                                                      (28)

Формулы (27) и (28) приведены к нормированным координатам (x, y) → (ξ, η):

Полученные матрицы жесткости для элементов и векторы ассамблируются в глобальную матрицу жесткости [K_глоб] и вектор внешних нагрузок {F}. Таким образом, можно записать основное уравнение метода КЭ, решение которого позволит определить вектор перемещений узлов {U}:

[K_глоб]{U} = {F}.                                                                        (31)

Преобразование узловых значений в элементарные

После нахождения узловых значений перемещений, деформаций и напряжений по описанному алгоритму следует перейти к определению напряжений каждого КЭ.

Процесс преобразования прост. Полученные значения для всех узлов, на которые ссылается элемент, усредняются по формуле (32) для вычисления итогового значения элемента σ^j_x :

                                                    (32)

где j– номер элемента;i– номер узла в соответствии с глобальной нумерацией узлов на рис. 5.

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

Рис. 6. Преобразование узловых напряжений в элементарные для первого элемента

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

                             (33)

Для проведения анализа на трещиностойкость определяли значения K_Ic, K_IIc, G_c, U_cв соответствии с выражениями (3), (4), (6) и (7).

Результаты и обсуждение

Для зарождения или распространения трещины необходимо выполнение условия U_c/G_c > 1 или    > σ_в, где  – значения главных напряжений элемента, вычисленные по формуле (33).

При выполнении одного из условий делается вывод, что в элементе образуется трещина. Это состояние КЭ считается предельным, локальная матрица жесткости элемента ослабляется на величину, равную 10^–4. Далее проводится повторный расчет с новой глобальной матрицей жесткости.

Предельные состояния КЭ образуются как от растягивающих напряжений (рис. 7, с 1-й по 13-ю итерации), так и от сжимающих напряжений (рис. 7, с 10-й по 13-ю итерации). В случае сжимающих напряжений локальные матрицы жесткости не ослабляли на величину 10^–4, так как трещина распространяется только от растягивающих напряжений в правой области. Предельные состояния от сжимающих напряжений показаны для демонстрации. Полное разрушение L-образной панели, соответствующее разделению конструкции на две части, зафиксировано на 13-й итерации расчета при соприкосновении левой и правой областей предельных состояний. На рис. 7 показано распределение предельных состояний на каждой итерации расчета.

Реализация описанного алгоритма, включая процедуру сборки глобальной матрицы жесткости, решение уравнения (31), расчет напряжений и критериев разрушения, выполнена на языке Python с использованием основных библиотек: sympy, scipy, math и numpy.

Распространение предельных состояний в описанном алгоритме хорошо согласуется с экспериментальными данными бетонной панели при визуальном наложении трещины на предельные состояния КЭ (рис. 8) [24].

Рис. 8. Экспериментальный фронт трещины бетонной панели [24]

Заключения

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

Работа выполнена в рамках реализации комплексного научного направления 2. «Фундаментально-ориентированные исследования, квалификация материалов, неразрушающий контроль» («Стратегические направления развития материалов и технологий их переработки на период до 2030 года») при поддержке ЦКП «Климатические испытания» НИЦ «Курчатовский институт» – ВИАМ.