Boundary integral equation method in the modeling of nonlinear deformation and failure of the 3D inhomogeneous media

Full Text

Введение. Решение нелинейных задач механики деформирования и разрушения для трёхмерных тел (элементов конструкций, деталей машин) с реальными геометрией и размерами остаётся сложной проблемой даже при наличии современных мощных ЭВМ и таких универсальных методов, как метод конечных элементов (МКЭ). Это связано с высокими требованиями, предъявляемыми к объёму памяти и быстродействию ЭВМ, большой трудоёмкостью метода и особенно при описании локальных особенностей решения, связанных с условиями нагружения, наличием контактных разрывов, дефектов типа трещин и др. Из-за малости размеров дефектов и других концентраторов по сравнению с размерами самих деформируемых тел последние приходится рассматривать в виде полубесконечных или бесконечных сред. Соответствующие краевые задачи переходят в разряд внешних, и требуются дополнительные подходы для учёта граничных условий «на бесконечности». В методе граничных интегральных уравнений (МГИУ), получившего в поISSN: 2310-7081 (online), 1991-8615 (print); doi: http://dx.doi.org/10.14498/vsgtu1292 © 2014 Самарский государственный технический университет. Образец цитирования: В. А. П е т у ш к о в, “Метод граничных интегральных уравнений в моделировании нелинейного деформирования и разрушения трехмерных неоднородных сред” // Вестн. Сам. гос. техн. ун-та. Сер. Физ.-мат. науки, 2014. № 2 (35). С. 96-114. 96 Метод граничных интегральных уравнений в моделировании . . . следние годы большое развитие, размерность рассматриваемых задач понижается на единицу, поскольку аппроксимация решения соответствующих интегральных уравнений строится только на границе области. Благодаря этому уменьшается объём исходной информации и появляется возможность описания с высокой точностью локальных особенностей решения внутри области. При этом в отличие от МКЭ автоматически моделируется (описывается) поведение деформируемой среды на бесконечности. Применение МГИУ основано на представлении решения в виде потенциалов простого и двойного слоя с последующим сведением краевой задачи к граничным интегральным уравнениям [1, 2]. Основные затруднения связаны с необходимостью вычисления граничных почти сингулярных интегралов, которые возникают при вычислениях значений физического поля вблизи границы или изучении тонких пограничных структур, и решением несимметричных плотно упакованных систем разрешающих алгебраических уравнений большого порядка. Для вычисления граничных интегралов к настоящему времени разработан ряд подходов, основанных или на адаптивных схемах с заданной точностью и сложных квадратурах переменного порядка или методах преобразований и квадратуре Гаусса, переходе к полярным координатам с применением аналитического и численного интегрирования [3, 4] и другие. На гладких явно параметризованных границах обычно применяются численные схемы с квадратурами высокого порядка. С вычислительной точки зрения все эти схемы могут быть более эффективными, если их использовать совместно с так называемыми быстрыми методами формирования и решения больших систем уравнений МГИУ [5]. Когда эти методы применимы, число операций O(N 2 ), обычно необходимых на формирование, и O(N 3 ) на решение систем уравнений, удаётся снизить до O(N ), где N - размер линейной системы, и тем самым решать задачи практически любого размера. Обстоятельный обзор современного состояния метода и направлений его развития дан в недавней работе [6]. В предлагаемой статье, основанной на полученных ранее результатах [4, 7, 8], применён так называемый прямой подход, когда в качестве потенциалов используются функции перемещений и усилий на границе тела. Интегральное уравнение строится на основе тождества Сомильяны с использованием фундаментальных решений Кельвина. Такой подход позволяет учитывать неоднородные материалы в теле, составные тела, рассматривать произвольные краевые условия. Для расширения его возможностей и повышения эффективности применительно к задачам нелинейного деформирования и разрушения трёхмерных тел используется метод дискретных областей, позволяющий рассматривать тела сложной геометрии с контрастными по размерам элементами и с локальными особенностями в геометрии и режимах нагружения. Такой подход органически вписывается в современные вычислительные технологии для многопроцессорных систем [9]. Выполнены решения известных трёхмерных задач нелинейного деформирования и разрушения упругопластических сред, характеризующие возможности метода. Представлены результаты моделирования нелинейных полей деформаций и напряжений, возникающих в двухслойном неоднородном соединении с поверхностной трещиной. 97 В. А. П е т у ш к о в 1. Постановка задачи и основные соотношения. Рассматривается нелинейная изотропная упругопластическая среда объёма V , занимающая в R3 произвольную, в общем случае, многосвязную область D, ограниченную поверхностью Липшица n+1 Sr, ∂D = r S n+1 где - поверхность, внешняя по отношению ко всем остальным. Область D (тело, конструкция), рис. 1, может быть неоднородной, составленной из различных материалов, и включать в себя поверхности разрыва - трещины Ωc (ˆα ), где xα , α = 1, 2 - локальная система координат, связанная с поверхx ˆ ностью разрыва. Здесь и далее используются обозначения и соглашения, принятые в тензорном исчислении. Рис. 1. Схема к постановке задачи [Figure 1. Scheme to the problem formulation] Пусть среда деформируется под действием массовых Xj (xi , t), xi ∈ D, и поверхностных pj (xi , t), xi ∈ ∂Dσ , сил, а также перемещений uj (xi , t), заданных на части поверхности c xi ∈ ∂Du , и температурного поля T (xi , t), xi ∈ D. Следовательно, ∂D = ∂Du ∪ ∂Dσ , при этом ∂Du ∩ ∂Dσ = ∅, а t ∈ Dt = = (0, τ ), где τ - продолжительность истории нагружения. Состояние среды в момент времени t = 0 может быть как недеформированным, так и включать в себя начальные или остаточные напряжения и деформации, в том числе технологического происхождения. Поведение среды при указанных воздействиях определяется полями перемещений uj (xi , t), деформаций εjk (xi , t) и напряжений σjk (xi , t) или их скоростями, соответственно uj , εjk , σjk . Полагая деформации малыми, а процесс ˙ ˙ ˙ нагружения достаточно медленным (квазистатическим), можем для любой переменной физического поля g(xi , t) записать g(xi , t) = ˙ 98 dg ∂g = , ∂t dt ∆g = g∆t, ˙ Метод граничных интегральных уравнений в моделировании . . . и тогда t играет роль индекса, характеризующего изменения в истории нагружения. Тензор скорости деформаций εjk представим в виде суммы упругой εe и ˙ ˙jk 0 составляющих: начальной εjk ˙ εjk = ˙ 1 (uj ,k + uk,j ) = εe + ε0 , ˙ ˙ ˙jk ˙jk 2 (1) причём ε0 в свою очередь может включать в себя как пластические деформа˙jk p ции εjk , так и любые другие, например, технологические ε∗ ; отнесём к числу ˙ ˙jk θ = αδ θ, где α - коэффици˙ перечисленных и температурные деформации εjk ˙ jk ˙ ент температурного расширения, δij - символ Кронекера, θ - разность температур. Соответствующее поле скоростей напряжений определяется обобщённым законом Гука σjk = Λjklm εlm - ε0 , ˙ ˙ ˙lm (2) где тензор упругости для изотропной среды Λjklm = λδij δkl + µ(δik δjl + δil δjk ), λ = Eν/(1 + ν)(1 - 2ν) и µ = E/2(1 + ν) - зависящие от температуры параметры изотропной среды, E - модуль Юнга, ν - коэффициент Пуассона. Для вычисления скоростей пластических деформаций εp , входящих в (2), ˙jk удобно использовать модифицированные соотношения теории течения упругопластических сред с анизотропным упрочнением [10]. С учётом классических условий течения Мизеса, разгрузки и нагружения полная деформация εjk представляется в виде εjk = εe + εp + ∆εp + εθ , jk jk jk jk и приращения пластической деформации ∆εp определяются из соотношения jk ∆εp = jk ∆εp i Ejk , Ei (3) где 1 Ejk = ejk - µρjk = (Sjk - ρjk ) + ∆εp jk 2 - девиатор, а 1/2 2 Ejk Ejk 3 - интенсивность так называемой смещенной деформации, причём Ei = 1 ejk = εjk - εll δjk , 3 ρjk = H εp , ˙ ˙jk H= dσi dεp i 1 Sjk = σjk - σll δjk , 3 , ∗ εi =εi σi = 3 Sjk Sjk 2 1/2 , 99 В. А. П е т у ш к о в ∆εp - интенсивность приращения пластических деформаций определяется в i этом случае на каждом n-ном шаге нагружения деформируемой среды диаграммой деформирования σi ∼ εp , полученной для соответствующей темпеi ратуры, и выражением n-1 ∆εp = 3µEi - σi i 3µ + H n-1 -1 . Деформирование рассматриваемой изотропной среды, как известно, может быть описано обобщённым уравнением Ламе ˙ ˙kj µuj ,kk + (λ + µ)uk,kj - 2µε0 ,j + λε0 ,j + Xj = 0, ˙ ˙ ˙kk (xi , t) ∈ D × Dt (4) с граничными условиями Неймана и Дирихле соответственно в форме pj = σjk nk = pb ˙ ˙ ˙j uj = ub ˙ ˙j на ∂Dσ × Dt , на ∂Du × Dt . (5) Здесь nj - компонента вектора внешней нормали к поверхности ∂D. Запятая перед индексом обозначает производную по соответствующей координате. Существование и однозначность решения краевой задачи (4), (5) u ∈ C(D) ∩ C 2 (D) следует из классических результатов [11]. При наличии в области D внутренних границ Γ (рис. 1), обусловленных соединением разнородных материалов или составных тел, условия (2) дополняются следующим: u1 = u2 , ˙j ˙j p1 = p2 ˙j ˙j на Γ × Dt , (6) где цифрами 1, 2 обозначены тела (материалы), находящиеся по обе стороны от границы Γ, которая в данном случае не является обычной физической. Если соединение допускает относительное смещение тел, то внутренняя граница Γ между ними рассматривается в виде первоначально совпадающих поверхностей Γ1 , Γ2 тел в зоне контакта. Отрыву (расслоению) общих точек этих поверхностей соответствуют условия равенства в них нулю усилий, т. е. p1 = p2 = 0, а скольжению с трением - ˙j ˙j следующие условия: u1 = u2 , ˙n ˙n p1 + p2 = 0, ˙n ˙n p1 + p2 = 0, ˙τ ˙τ -κp1,2 ˙n p1 , τ p2 ˙τ κp1,2 ˙n на Γ1 , Γ2 , (7) где n и τ - нормальное и касательные направления к поверхностям Γ1 , Γ2 в точке контакта, κ - коэффициент трения. В этом случае задача (4), (5) оказывается нелинейной и для упругого поведения деформируемой среды. 2. Интегральные представления краевой задачи. Краевой задаче (4)-(7) могут быть поставлены в соответствие различные граничные интегральные уравнения, имеющие различные свойства погружения, и, как следствие, различные свойства используемых для их решения численных схем. Все они основаны на представлении решения uj (xi , t) краевой задачи в виде поверхностного потенциала с заданной плотностью и свойствах непрерывности подобных потенциалов [1, 2, 12]. 100 Метод граничных интегральных уравнений в моделировании . . . Одно из таких представлений, полученное на основе обобщённого тождества Сомильяны, ведёт к следующему интегральному сингулярному уравнению: Cij (x)uj (x) = - ˙ Uij (x, y)pj (y)dsy + ˙ Pij (x, y)uj (y)dsy + ˙ ∂D ∂D ˙ Uij (x, y)Xj (y)dv - + Σjki (x, y)ε0 dv, (8) ˙jk Vp V где y ∈ ∂D - точка поверхности, x ∈ R3 \ ∂D - точка (источник) внутри области D; V , Vp - соответственно области с отличными от нуля объёмными силами и пластического деформирования и/или ненулевой разности температур; Cij (x) = δij , если x ∈ D, и Cij (x) = 0.5δij , если x ∈ ∂D, при условии её гладкости; Uij , Pij , Σjki - тензоры, соответствующие фундаментальному решению Кельвина-Сомильяны, определяемые выражениями Uij = [16πµ(1 - ν)r]-1 [(3 - 4ν)δij + r,i r,j ] , -1 ∂r Pij = - 8π(1 - ν)r2 [(1 - 2ν)δij + 3r,i r,j ] - ∂n - (1 - 2ν)(r,i nj - r,j ni ) , Σjki = - 8π(1 - ν)r2 -1 (9) (1 - 2ν)(δij r,k + δik r,j - δjk r,i )+ + 3r,i r,j r,k , r = x-y и r,i = ∂r = r-1 (yi - xi ). ∂xi Существование и единственность решений уравнения (8) с условиями (5) для задач упругости установлено в классе непрерывных по Гельдеру функ1 1 ций C 1,α (D) ⊂ W2 (D), где W2 (D) - пространство Соболева с параметром 0 < α 1 [1, 2, 12]. Для нерегулярных областей, в том числе включающих поверхности разрыва (трещины), такие исследования выполнены, например в [13], на основе теории псевдодифференциальных операторов и метода Винера-Хопфа. Для рассматриваемых неоднородных (кусочно-однородных) сред уравнение (8) записывается отдельно для каждой из подобластей, входящих в область D и имеющих общие границы. В целом для области соответствующие уравнения объединяются с помощью условий (6), (7) на этих границах. Уравнение (8) применимо и для внешних краевых задач, когда область D дополняется до всего трехмерного пространства, если из его левой части вычесть интеграл с Pij по области дополнения. Выражение для напряжений σjk (xi , t) может быть получено дифферен˙ цированием (8) с учётом (1)-(3), и для точек внутри области D имеет вид uijk pk ds - ˙ σij = ˙ ∂D + Vp ˙ uijk Xk dV + pkij uk ds + ˙ ∂D V ˙ Σijkl εp + αδkl θ dV - ˙kl 7 - 5ν 4µ(1 + ν) ˙ 2µεp - ˙ij δij αθ, (10) 15(1 - ν) 3(1 - ν) 101 В. А. П е т у ш к о в где, принимая во внимание (9), ukij = -Σkij = λδij ukl,l , pkij = λδij pkl,l . Тензор Σijkl определяется соотношением Σijkl = 2µ 3(1 - 2ν)(δij r,k r,l + δkl r,i r,j )+ + 3ν(δli r,j r,k + δjk r,l r,i + δik r,l r,j + δjl r,i r,k ) - 15r,i r,j r,k r,l + + (1 - 2ν)(δik δlj + δjk δli ) + (1 - 4ν)δij δkl (8π(1 - ν)r3 )-1 . Для вычисления напряжений на границе области D выражение (10) не применимо из-за сингулярного характера входящих в него интегралов. Как будет показано, напряжения в точках xi ∈ ∂D могут быть получены через граничные усилия и производные от граничных перемещений по касательным направлениям. Для этого необходима параметризация геометрии границы ∂D и аппроксимации решения на её поверхности. Аппроксимация решения также необходима и внутри области D для определения последних двух интегралов по внутренним объемам V , Vp в выражениях (8), (10) и параметров механики разрушения в окрестности трещины Ωc (ˆα ). x Полагая рассматриваемые поверхности области D кусочно-гладкими, численное решение уравнения (8) выполним с использованием M граничных элементов ∆k так, чтобы M h ∂D ≈ ∂D = ∆k и ∆k ∩ ∆l = ∅, k = l. k=1 Параметрическое представление геометрии и искомых функций ui (xj , t) и ˙ pi (xj , t) на каждом из них выберем подобно МКЭ в следующем виде: ˙ xh := Nr (ζγ )xr , i i uh := Nr (ζγ )ur , ˙i ˙i ph := Nr (ζγ )pr , ˙i ˙i ζγ ∈ ∆k ⊂ R2 , (11) где Nr (ζγ ), r = 1, 2, . . . , 8 - квадратичная функция формы локальных координат ζγ ∈ (-1,1), принадлежащая к Лагранжеву семейству конечных элеk,m ментов Sh (∂Dh ) в смысле [14]. Параметрическое представление (11) используется для вычисления поверхностных интегралов, входящих в (8), (10). Вычисление объёмных интегралов для массовых сил и по области начальных деформаций ε0 осуществ˙jk ляется с использованием объёмных 8 или 20-узловых конечных элементов ∆k того же семейства с квадратичным изменением сил и деформаций, аналогично (11), но с функцией формы вида L k,m Nr (ζγ , η) ∈ Sh (D), D ≈ Dh = ∆k , k=1 η ∈ (-1, 1) - третья локальная координата, нормальная к двум первым ζγ , r = 1, 2, . . . , 8, . . . , 20. Отметим, что размер области с пластическими деформациями заранее не известен и устанавливается в процессе решения. 102 Метод граничных интегральных уравнений в моделировании . . . При наличии особенностей решения внутри и на границе области с известным характером сингулярности, её можно встроить в пространство пробных k,m функций Sh (∂Dh ) подобно МКЭ или за счёт подходящего построения неоднородной сетки ∂Dh или расширения самого пространства путём введения специальных сингулярных функций, т. е. использования сингулярных элементов ∆k [4, 8]. Заменяя далее интегралы по границе и объёму соответствующими суммами по граничным и объёмным элементам, вместо (8), (10) получаем их численные аналоги: ˙ ˙ [A]X = F + [D]{ε0 }, ˙ ˙ + F + [D ]{ε0 } + [G]ε0 . ˙ σ = -[A ]X ˙ ˙ ˙ (12) (13) Здесь [A], [A ] - матрицы коэффициентов при векторе граничных неизвест˙ ных; X - вектор, содержащий неизвестные усилия и перемещения на границе; ˙ ˙ векторы F , F содержат члены, определяемые отличными от нуля заданными значениями граничных перемещений и усилий, а также действием объёмных сил. Матрицы [D], [D ] включают в себя объёмные интегралы в (8), (10), квазидиагональная матрица [G] определяется наличием свободных от интеграла членов в (10). Для ускорения матрично-векторных перемножений могут быть использованы выше упомянутые быстрые методы [5]. Выражения для матриц и векторов в уравнения (12), (13) подробно представлены в [4], и здесь ввиду их громоздкости не приводятся. Вычисления входящих в них интегралов основано на численных квадратурах Гаусса. Объём и точность вычислений при этом зависят от порядка формулы и минимального расстояния между точкой наблюдения и узлом рассматриваемого элемента ∆k . Если точка наблюдения совпадает с одним из узлов граничного или объёмного элементов, т. е. r → 0, интегралы в (8), (10) становятся сингулярными. Для их вычисления используется переход соответственно к полярным и сферическим координатам и квадратуры Гаусса по радиальным и угловым координатам. Сильно сингулярные интегралы с особенностью O(r-3 ) рассматриваются в смысле главного значения Коши. Их вычисление по объёмным элементам строится путём перехода к сферическим координатам и выделением локальной зоны вблизи сингулярной точки, по которой проводится расчёт аналитически по радиусу и численно квадратурой Гаусса по угловым координатам. Тем самым получаем дополнительный вклад в диагональные блоки матрицы [D ] в выражении (13). Для кусочно-однородных (составных) сред уравнения (12), (13) записываются отдельно для каждой из сред, и на общей поверхности между ними Γ = Γ1 ∪ Γ2 вводится два слоя узлов и элементов. В случае идеального контакта на первом из них в качестве граничных условий задаются неизвестные перемещения и нулевые усилия, на втором - неизвестные усилия и нулевые перемещения. В случае проскальзывания с трением на первом слое задаются неизвестные компоненты перемещений - нормальные и касательные, а на втором - не известные нормальные усилия и касательные перемещения. При этом условия (6), (7) учитываются при формировании матриц [A] и [A ]. 103 В. А. П е т у ш к о в ˙ Исключив из (12), (13) вектор граничных неизвестных X, получаем окончательные выражения МГИУ для выбранного интегрального представления краевой задачи (4)-(7): ˙ X = [K]{ε0 } + {m}, ˙ ˙ (14) σ = [B]{ε0 } + {n} + [G]ε0 . ˙ ˙ ˙ ˙ (15) Здесь [K] = [A]-1 [D], {m} = [A]-1 F, ˙ [B] = [D ] - [A ][K], {n} = F - [A ]{m}. ˙ ˙ Отметим, что {m} и {n} представляют собой решение для упругой среды ˙ ˙ соответственно для вектора граничных неизвестных и вектора напряжений. Для определения обратной матрицы [A]-1 в её элиминативной форме и вектора {m} используется алгоритм гауссова типа подобно [15] с блочным ис˙ ключением и выбором главного элемента в блоке. При этом нулевые блоки матрицы, появляющиеся при анализе составных сред, не формируются, не хранятся и не участвуют в процессе исключения, что позволяет решать прикладные задачи большой размерности. Сложные истории нагружения деформируемой среды на временном слое Dt = (0, τ ) представим в виде совокупности простых нагрузок с помощью зависящих от времени параметров силового λp (t) и температурного λθ (t) нагружений. Разрешающие уравнения (14), (15) в этом случае приобретают вид k ˙ X tk N p = λj (tk ){mj }, p [K]{ε (tj )} + λθ (tk )[K]α{θ} + ˙ j=1 (16) j=1 k σ tk ([B] + [Gε ]){εp (tj )} + λθ (tj )([B] + [Gθ ])α{θ}+ ˙ = j=1 N λpj (tj ){nj }, + (17) j=1 где t ∈ Dt , {mj } и {nj } для любого шага j = 1, 2, . . . , N полного нагружения определяются, как и в (14), (15). Такой подход позволяет рассматривать как повторяющиеся статические, так и циклические условия нагружения нелинейно деформируемой среды с анизотропным упрочнением, используемым в (3). Моделирование нелинейного деформирования среды под действием j-того уровня нагружения осуществляется на основе классической схемы метода последовательной линеаризации в форме начальных деформаций [8, 10]. Приращения пластической деформации определяются на каждой итерации из соотношений (9) с использованием реальных диаграмм деформирования σ∼εp ∼θ, построенных для каждой из рассматриваемых сред. 3. Вычисление напряжений на границе тела. Интегральное выражение (10) не применимо для вычисления напряжений на границе нелинейно деформируемого тела из-за сингулярного характера входящих в него интегралов с особенностью O(r-2 , r-3 ). Поэтому воспользуемся вычисленными на ней усилиями и касательными производными от граничных перемещений. 104 Метод граничных интегральных уравнений в моделировании . . . Запишем выражения для напряжений и усилий через производные от перемещений, используя (1), (2) и первую формулу (5), σjk = λδjk ui,i + µ (uj ,k + uk,j ) , ˙ ˙ ˙ ˙ pj = σjk nk = λδjk nk ui,i + µ (uj ,k nk + uk,j nk ) . ˙ ˙ ˙ ˙ ˙ (18) Граничные усилия pj (xi , t), а также перемещения uj (xi , t) и их производные ˙ ˙ по локальным координатам ζγ с учётом аппроксимации (11) могут быть получены из решения системы уравнений (12). Затем через них с использованием (18) определяются напряжения на границе ∂D рассматриваемой области D. Предлагаемый подход оказывается весьма эффективным с вычислительной точки зрения, поскольку обладает высокой точностью и не требует дополнительного вычисления производных от интегралов по области пластического деформирования, нет необходимости и вычисления матрицы D в выражении (15) в процессе решения нелинейной задачи. Эти особенности подхода делают целесообразным его обобщение на нелинейно деформируемые среды, расчётные области которых состоят из элементов с сильно различающимися размерами или включают в себя произвольно ориентированные дефекты - трещины Ωc (ˆα ) с особым характером нагружеx ния и взаимодействия их берегов. В этом случае область, занимаемая деформируемой средой, рассматривается в виде дискретных подобластей, образованных соответствующими произвольно ориентированными фиктивными границами, например Γ2 на рис. 1, с заданными на них условиями типа (6) или (7), если необходимо учитывать взаимодействия берегов трещины. 4. Приложение к задачам механики разрушения. Процессы деформирования и разрушения рассматриваемых сред взаимосвязаны и протекают одновременно на различных уровнях их структуры. Предметом изучения классической линейной или нелинейной механики разрушения являются так называемые макротрещины с размерами, обнаруживаемыми существующими методами неразрушающего контроля. Выявляемые этими методами внутренние или поверхностные дефекты принимаются в виде эквивалентных по площади круглой или эллиптической формы трещин - поверхностей разрыва деформируемой среды. Наряду с изучением поведения таких трещин необходимо, очевидно, учитывать высокую вероятность наличия не обнаруживаемых трещин меньших размеров особенно в зонах высокой концентрации напряжений и деформаций. Поэтому разработка эффективных методов моделирования поведения реальных конструкций с подобными трещинами остаётся актуальной задачей. К преимуществам используемого в этом случае МГИУ следует отнести, прежде всего, возможность сколь угодно точного определения напряжённых и деформированных состояний в окрестности трещины, при наличии решения (14), (15) на поверхности рассматриваемого тела. Более того, анализ подрастания трещины в соответствии с историей нагружения также не вызывает затруднений. В практических приложениях для широкого круга задач механики разрушения пластические деформации обычно локализованы в окрестности трещин. Имеет место так называемое маломасштабное течение. Поэтому применимы подходы линейной механики разрушения, основанные на уравнениях 105 В. А. П е т у ш к о в теории упругости. Чаще всего используются коэффициенты интенсивности напряжений (КИН), которые являются мерой доминирующих напряжений в окрестности фронта трещины, определяющих её поведение. В случае объёмного напряженного состояния реализуется смешанный режим разрушения - сдвиговый и нормального отрыва, и необходимо вычислять все три коэффициента в зависимости от положения вдоль фронта трещины (рис. 2). Рис. 2. Схема разрушения вдоль фронта трещины при объёмном напряжённом состоянии [Figure 2. Scheme of destruction along the crack front under triaxial stress conditions] На границе трещины Ωc (ˆα ) искомые решения (14), (15) претерпевают x разрыв, а вблизи её фронта имеют особенность типа r1/2 и r-1/2 соответственно для перемещений и усилий. Для учёта этих особенностей используется то же представление (11), но со специальным набором базисных функций Nr (ζγ ), как, например в [7, 16]. КИН вычисляются непосредственно через перемещения ϕ 3ϕ - cos + 2 2 3ϕ 1 ϕ , + K2 (9 - 8ν) sin + sin 2 2 4E 2r 1/2 ϕ 3ϕ u2 = (1 + ν) K1 (7 - 8ν) sin - sin + π 2 2 ϕ 3ϕ 1 , + K2 (3 - 8ν) cos + cos 2 2 4E 2(1 + ν)(2r/π)1/2 ϕ u3 = K3 sin , E 2 которые получены вблизи фронта трещины с помощью МГИУ. Здесь индексы 1, 2, 3 соответствуют направлениям I, II, III на рис. 2. Полагая φ = π, как на рис. 2, получаем выражения для КИН, соответствующие трем типам упругого или квазиупругого (маломасштабное течение) u1 = 106 2r π 1/2 (1 + ν) K1 (5 - 8ν) cos Метод граничных интегральных уравнений в моделировании . . . разрушения: K1 = µu2 (π/2r)1/2 /(1 - ν), K2 = µu1 (π/2r)1/2 /(1 - ν), (19) K3 = µu3 (π/2r)1/2 . Наличие развитых зон пластического деформирования в окрестности трещин с размером, равным или большим характерного размера трещины, делает соотношения (19) не применимыми. В этом случае решение краевых задач нелинейной механики разрушения строится средствами МГИУ в соответствии с (8), (17). В качестве параметров разрушения, контролирующих поведение трещины, используются энергетические G- или J-интегралы и деформационные [16, 17], при этом наиболее подходящим, имеющим ясный физический смысл и доступным измерению является δ-раскрытие трещины. Подход, развиваемый здесь для описания составных кусочно-однородных ограниченных сред, может быть легко обобщён на решение трёхмерных задач механики разрушения с произвольно ориентированными в пространстве трещинами. Тело с трещиной рассматривается в этом случае в виде составного тела, внутренняя (фиктивная) граница которого включает в себя поверхность разрыва - трещину Ωc (ˆα ) с граничными условиями (6) и (10) с учётом x нагруженности берегов трещины и их взаимодействия, как, например, в когезионных моделях разрушения [18]. 5. Численные результаты. Важным вопросом при моделировании (и особенно для нелинейных процессов деформирования и разрушения сред) является математическое обоснование точности и вычислительной устойчивости получаемых результатов. Применительно к МГИУ такое обоснование может быть выполнено на основе интенсивно развиваемой теории псевдодифференциальных операторов [19]. Однако для трёхмерных задач и используемого коллокационного приближения до сих пор отсутствуют соответствующие асимптотические оценки. Поэтому анализ качества получаемых решений будем оценивать на примерах решения модельных задач с известными результатами, полученными другими методами. Вначале приведём решение известной задачи Киpша в объёмной постановке с диаметром отверстия, равным толщине пластины. Размеры и история нагружения показаны на рис. 3 a. В силу симметрии рассматривается 1/8 часть полосы. Материал полосы - аустенитная сталь (E = 2.0 · 105 МПа, ν = 0.3, σ0.2 = m = 250 МПа, параметр упрочнения m = 0.25; σi /σ0.2 = εp E/σ0.2 - диаi грамма деформирования). Однородные растягивающие напряжения σ0 приложены к двум противоположным граням пластины. Максимальная величина прикладываемой нагрузки σ0 = 0.6σ0.2 . История нагружения пластины на Dt = (0, τ ) также приведена на рис. 3 a. Результаты моделирования представлены на рис. 3 и 4 для двух моментов времени, соответствующих максимальной нагрузке (рис. 3 a, b) и полной разгрузке (рис. 4). В случае упругого поведения материала задача имеет аналитическое решение. Для нелинейного деформирования воспользуемся сравнительными результатами, полученными МКЭ [17]. На рис. 3 a результаты решения МГИУ (линии 1) показаны вместе с решением, полученным МКЭ (линии 2), и соот107 108 a деление деформаций b (BIEM), 2 - Finite Element Method (FEM)); b) the strain distribution] [Figure 3. 3D Kirsch’s problem: a) the geometry, the area of plastic deformation, and the stress distribution (1 - Boundary Integral Equation Method Рис. 3. Объёмная задача Кирша: a) геометрия, зона пластического деформирования и распределение напряжений (1 - МГИУ, 2 - МКЭ); b) распре- В. А. П е т у ш к о в Метод граничных интегральных уравнений в моделировании . . . a b Рис. 4. Распределение остаточных напряжений [Figure 4. Residual stresses distribution] ветствуют максимальной нагрузке в истории нагружения. Выполненное моделирование процесса деформирования пластины позволило выявить особенности формирования вблизи отверстия зон нелинейных деформаций (рис. 3, b) и остаточных напряжений (рис. 4), где приведены распределения остаточных напряжений в номинальном сечении на срединной плоскости (рис. 4, a) и вдоль образующей контура отверстия по толщине пластины (рис. 4, b). Аппроксимация поверхности и внутренней пластической области в пластине выполнена с использованием 20 поверхностных восьми узловых граничных элементов (62 узла) и 12 объёмных 20-узловых элементов. Для сравнения укажем, что в решении, полученном МКЭ, использовалось 26 объёмных 20-узловых элементов (900 уравнений) с тем, чтобы получить решение для коэффициента концентрации напряжений как на поверхности, так и на срединной плоскости вблизи отверстия пластины с погрешностью менее 3 % по сравнению с аналитическим. Погрешность используемого здесь МГИУ,в том числе вместе с методом дискретных областей, не превышает 1 %. Хорошее совпадение результатов вблизи отверстия получено и для задачи пластичности. В остальной области более точными следует считать, по-видимому, результаты, полученные МГИУ. Увеличение числа граничных элементов вдвое не оказало заметного влияния на точность полученных результатов при 16-точечной (по 4 в каждом направлении) квадратуре Гаусса вычисления интегралов. Отметим, что применение интегральных представлений решения в этом случае позволило с высокой точностью выявить роль объемности напряженных и деформированных состояний при оценке коэффициентов их концентрации вблизи отверстия. Максимальные значения напряжений и размеры зон пластических деформаций достигаются на срединной плоскости, откуда, как показывают эксперименты, и начинается разрушение полосы. В качестве другого модельного примера, имеющего прикладное значение и полученные ранее результаты, рассмотрим объёмную задачу механики разрушения для трещины нормального разрыва полуэллиптической формы на поверхности полосы. Геометрия полосы и трещины вместе с расчётной схемой приведены на рис. 5; здесь B/H = 2; H/c = 2.5. Максимальная нагрузка и механические свойства материала такие же, что и в предыдущей задаче. Рассматривались два случая полуэллиптической трещины с соотношени109 В. А. П е т у ш к о в ем полуосей a/c = 1/2 и более нагруженного с соотношением a/c = 2/3, полученные результаты для которых отмечены на рис. 5 цифрами 1 и 3 соответственно. Аппроксимация геометрии полосы с трещиной выполнена с использованием 24 восьми узловых граничных элементов (74 узла) и десяти объемных 20узловых элементов для последующего вычисления размеров предполагаемой зоны пластического деформирования в окрестности трещины. Для описания особенностей решения вблизи фронта трещины использовано 4 сингулярных граничных элемента с представлением [4]. Для трещины с размерами a/c = 0.5 и упругого поведения материала выполнен сравнительный анализ полученного распределения K1 вдоль фронта трещины (линия 1, рис. 5) с аналогичным распределением, заимствованным из [20] (линия 2, рис. 5). Столь же хорошее совпадение результатов с решением МКЭ [17] получено и для нелинейной задачи для случая a/c = 2/3. Максимальные значения величины КИН на графиках соответствуют условию плоской деформации, реализуемому в самой глубокой точке трещин. На поверхности пластины реализуется условие плоского напряженного состояния. Существенное различие результатов в обоих случаях наблюдается лишь в достаточно тонком слое в окрестности трещины вблизи поверхности полосы при φ = 90◦ . Выявленное здесь повышение КИН проявляется обычно при экспериментальном изучении роста трещины в условиях циклического упругопластического деформирования, обнаружить которое другими методами, например МКЭ, оказывается затруднительно. Зона пластических деформаций в окрестности полуэллиптической трещины a/c = 2/3 оказалась весьма локализованной. Пространственный характер её расположения выделен цветом (заливкой) на рис. 5. В качестве примера неоднородного (составного) тела рассмотрим задачу о растяжении двухслойной полосы с полуэллиптической трещиной при наличии идеального контакта между слоями, описываемого условиями (6). Полоса изготовлена из высоколегированной стали с антикоррозионной наплавкой из нержавеющей стали. Подобного рода плакированные соединения широко используются в химии и ядерной энергетике. Характеристики основного металла (материал 1 - реакторная сталь): E1 = 2.15 · 105 МПа, ν1 = Рис. 5. Полоса с поверхностной трещиной [Figure 5. The strip with a surface crack] 110 Метод граничных интегральных уравнений в моделировании . . . Рис. 6. Плакированное соединение с трещиной [Figure 6. Plating compound with crack] Рис. 7. Распределение напряжений и деформаций в окрестности трещины [Figure 7. Stress and strain distribution in the vicinity of the crack] 111 В. А. П е т у ш к о в 1 = 0.3, σy = 430 MПa; характеристики нержавеющей наплавки (материал 2): 2 E2 = 2.05 · 105 МПа, ν2 = 0.27, σy = 180 MПa. Размеры полосы с дефектом (трещиной) приведены на рис. 6. Здесь, как и в предыдущем примере, B/H = 2, H/c = 2.5. Толщина плакированного слоя h = 0.1H, размеры трещины a/c = 0.9, где a = 1.2h, приняты эквивалентными по площади дефекту, выявленному дефектоскопией в реальной конструкции. Зоны пластических деформаций, полученные в окрестности трещины для 1 1 двух случаев нагружения полосы σ0 = 0.3σy и 0.4σy , показаны на рис. 6 и отмечены цифрами 1 и 2 соответственно. Основной металл (материал 1) при этих нагрузках работает в упругой области. В качестве параметра разрушения выбрано раскрытие трещины δ, его значения, вычисленные для указанных случаев нагружения, также приведены на рис. 6. Трещина начинает распространяться при условии δ δc , где критическое значение δc определяется экспериментально. Представляющие практический интерес распределения напряжений и деформаций в окрестности трещины вдоль осей a и c, отражающие влияние неоднородности свойств плакированного соединения и пластического течения в наплавке, представлены на рис. 7. Этот пример кроме его практической значимости является хорошей демонстрацией возможностей МГИУ в описании тонких локальных особенностей процессов деформирования и разрушения неоднородных сред. Заключение. Разработан метод интегрального представления решения трёхмерных задач термо-упругопластического деформирования и разрушения неоднородных (составных) тел сложной формы с локальными особенностями типа трещин. Обобщением подхода на расчётные области с произвольно ориентированными дефектами, сильно различающимися размерами отдельных элементов и/или свойствами материалов, является предлагаемый метод дискретных областей. Учитываются сложные истории термо-силового нагружения составных тел и возможность их относительного смещения на границах в зоне контакта. Получены решения нелинейных задач деформирования и разрушения трёхмерных тел и дано их сравнение с известными численными или аналитическими решениями. Во всех случаях точность полученных результатов достаточно хорошая.

About the authors

Vladimir A Petushkov

A. A. Blagonravov Mechanical Engineering Institute RAS

Email: pva_ras@inbox.com
(Dr. Phys. & Math. Sci.), Professor, Lab. of Mathematics Modeling. 4, M. Khariton’evskii per., Moscow, 101990, Russian Federation

References

Supplementary files