Метод граничных интегральных уравнений в моделировании нелинейного деформирования и разрушения трехмерных неоднородных сред Текст научной статьи по специальности «Физика»

Вестн. Сам. гос. техн. ун-та. Сер. Физ.-мат. науки. 2014. № 2(35). С.96—114

УДК 539.42

МЕТОД ГРАНИЧНЫХ ИНТЕГРАЛЬНЫХ УРАВНЕНИЙ В МОДЕЛИРОВАНИИ НЕЛИНЕЙНОГО ДЕФОРМИРОВАНИЯ И РАЗРУШЕНИЯ ТРЕХМЕРНЫХ НЕОДНОРОДНЫХ СРЕД

В. А. Петушков

Институт машиноведения им. А. А. Благонравова РАН,

Россия, 101990, Москва, М. Харитоньевский пер., 4.

Метод граничных интегральных уравнений применяется для решения нелинейных задач термо-упругопластического деформирования и разрушения неоднородных объёмных тел сложной формы. Решение строится на основе обобщённого тождества Сомильяны и метода последовательной линеаризации в форме начальных пластических деформаций. Приращения пластической деформации определяются на основе теории течения упрочняющихся упругопластических сред с использованием модифицированных соотношений Прандтля—Рейсса. Рассмотрены случаи сложного термосилового нагружения находящихся в контакте составных кусочно-однородных сред. Для описания процессов нелинейного деформирования и разрушения тел сложной геометрии с локальными особенностями используется разработанный ранее метод дискретных областей. Представлены решения трёхмерных нелинейных задач механики деформирования и разрушения, имеющих практическое значение.

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

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

Из-за малости размеров дефектов и других концентраторов по сравнению с размерами самих деформируемых тел последние приходится рассматривать в виде полубесконечных или бесконечных сред. Соответствующие краевые задачи переходят в разряд внешних, и требуются дополнительные подходы для учёта граничных условий «на бесконечности».

В методе граничных интегральных уравнений (МГИУ), получившего в по-

ISSN: 2310-7081 (online), 1991-8615 (print); doi: http://dx.doi.org/10.14498/vsgtu1292 © 2014 Самарский государственный технический университет.

Образец цитирования: В. А. Петушков, “Метод граничных интегральных уравнений в моделировании нелинейного деформирования и разрушения трехмерных неоднородных сред” // Вестн. Сам. гос. техн. ун-та. Сер. Физ.-мат. науки, 2014. № 2 (35). С. 96-114. doi: 10.14498/vsgtu1292.

Сведения об авторе: Владимир Алексеевич Петушков (д.ф.-м.н.), профессор, лаборатории математического моделирования.

E-mail address: pva_ras@inbox.com

96

Метод граничных интегральных уравнений в моделировании ...

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

Применение МГИУ основано на представлении решения в виде потенциалов простого и двойного слоя с последующим сведением краевой задачи к граничным интегральным уравнениям [1,2]. Основные затруднения связаны с необходимостью вычисления граничных почти сингулярных интегралов, которые возникают при вычислениях значений физического поля вблизи границы или изучении тонких пограничных структур, и решением несимметричных плотно упакованных систем разрешающих алгебраических уравнений большого порядка.

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

С вычислительной точки зрения все эти схемы могут быть более эффективными, если их использовать совместно с так называемыми быстрыми методами формирования и решения больших систем уравнений МГИУ [5]. Когда эти методы применимы, число операций O(N2), обычно необходимых на формирование, и O(N3) на решение систем уравнений, удаётся снизить до O(N), где N — размер линейной системы, и тем самым решать задачи практически любого размера.

Обстоятельный обзор современного состояния метода и направлений его развития дан в недавней работе [6]. В предлагаемой статье, основанной на полученных ранее результатах [4,7,8], применён так называемый прямой подход, когда в качестве потенциалов используются функции перемещений и усилий на границе тела. Интегральное уравнение строится на основе тождества Сомильяны с использованием фундаментальных решений Кельвина. Такой подход позволяет учитывать неоднородные материалы в теле, составные тела, рассматривать произвольные краевые условия.

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

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

97

В. А. Петушков

1. Постановка задачи и основные соотношения. Рассматривается нелинейная изотропная упругопластическая среда объёма V, занимающая в R3 произвольную, в общем случае, многосвязную область D, ограниченную поверхностью Липшица

n+1

dD = U Sr,

r

где Sn+1 —поверхность, внешняя по отношению ко всем остальным. Область D (тело, конструкция), рис. 1, может быть неоднородной, составленной из различных материалов, и включать в себя поверхности разрыва — трещины Qc(Xa), где ха, а = 1, 2 — локальная система координат, связанная с поверхностью разрыва. Здесь и далее используются обозначения и соглашения, принятые в тензорном исчислении.

Рис. 1. Схема к постановке задачи [Figure 1. Scheme to the problem formulation]

Пусть среда деформируется под действием массовых Xj(xl,t), хг £ D, и поверхностных pj(xi,t), хг £ dDa, сил, а также перемещений Uj(xi,t), заданных на части поверхности c хг £ dDu, и температурного поля T(хг,Ь), хг £ D. Следовательно, dD = dDu U 3Da, при этом dDu П 3Da = 0, а t £ Dt = = (0, т), где т — продолжительность истории нагружения. Состояние среды в момент времени t = 0 может быть как недеформированным, так и включать в себя начальные или остаточные напряжения и деформации, в том числе технологического происхождения.

Поведение среды при указанных воздействиях определяется полями перемещений Uj (хг, t), деформаций ejk(хг, t) и напряжений ajk(хг, t) или их скоростями, соответственно Uj, j, &jk. Полагая деформации малыми, а процесс нагружения достаточно медленным (квазистатическим), можем для любой переменной физического поля д(хг,t) записать

д(х^)=I=!• Дд=дАг-

98

Метод граничных интегральных уравнений в моделировании ...

и тогда t играет роль индекса, характеризующего изменения в истории нагружения.

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

£ jk = ^ (uj,k + uk,j) = j + £°jk, (!)

причём £0k в свою очередь может включать в себя как пластические деформации £pk, так и любые другие, например, технологические £*k отнесём к числу перечисленных и температурные деформации j = abjk0, где а — коэффициент температурного расширения, bij — символ Кронекера, 0 — разность температур.

Соответствующее поле скоростей напряжений определяется обобщённым законом Гука

&jk — A-jklm i&lm £lm) , (2)

где тензор упругости для изотропной среды

Лjklm = Лbij bkl + Ц-(bik bjl + bilbjk),

Л = Ev/( 1 + v)(1 — 2v) и ^ = E/2(1 + v) — зависящие от температуры параметры изотропной среды, E — модуль Юнга, v — коэффициент Пуассона.

Для вычисления скоростей пластических деформаций j,, входящих в (2), удобно использовать модифицированные соотношения теории течения упругопластических сред с анизотропным упрочнением [10]. С учётом классических условий течения Мизеса, разгрузки и нагружения полная деформация £jk представляется в виде

£jk = £jjk + £jk + A£jk + £%,

и приращения пластической деформации Aj, определяются из соотношения

(3)

A£p

AeP. = A£i Ejk,

jk

Ei

где

Ejk ejk 2 №pjk (Sjk pjk) + A£jk

девиатор, а

Ei — ^3 Ejk EjkJ

\ 1/2

— интенсивность так называемой смещенной деформации, причём

ejk — £jk з £llbjk, Sjk — &jk з &llbjk,

Pjk = H£

p, jk,

H

daj

d£pp

&i — I - SjkSjk J

\ 1/2

£i=£

99

В. А. Петушков

Дер — интенсивность приращения пластических деформаций определяется в этом случае на каждом n-ном шаге нагружения деформируемой среды диаграммой деформирования а г ~ еР, полученной для соответствующей температуры, и выражением

дер = (3рЕ% - а™-1) (3^ + Hга-1)-1.

Деформирование рассматриваемой изотропной среды, как известно, может быть описано обобщённым уравнением Ламе

№uj,kk + (А + У)йk,kj — {2^e°kk,j + Aekj,j) + Xj = 0, (хг, t) G D x Dt (4)

с граничными условиями Неймана и Дирихле соответственно в форме

pj = ajknk = pj на 0Da x Dt,

, (5)

Uj = uj на dDu x Dt.

Здесь nj — компонента вектора внешней нормали к поверхности dD. Запятая перед индексом обозначает производную по соответствующей координате. Существование и однозначность решения краевой задачи (4), (5) и G C(D) П C2(D) следует из классических результатов [11].

При наличии в области D внутренних границ Г (рис. 1), обусловленных соединением разнородных материалов или составных тел, условия (2) дополняются следующим:

и] = и2, р] = р2 на Г x Dt,

(6)

где цифрами 1, 2 обозначены тела (материалы), находящиеся по обе стороны от границы Г, которая в данном случае не является обычной физической. Если соединение допускает относительное смещение тел, то внутренняя граница Г между ними рассматривается в виде первоначально совпадающих поверхностей Г], Г2 тел в зоне контакта.

Отрыву (расслоению) общих точек этих поверхностей соответствуют условия равенства в них нулю усилий, т. е. р] = р2 = 0, а скольжению с трением — следующие условия:

ип ип, рп + рп 0,

р] + Р2 = 0, -крП2 < р],

Pi < Кр>П’2

на Г1, Г2,

(7)

где n и т — нормальное и касательные направления к поверхностям Г], Г2 в точке контакта, к — коэффициент трения. В этом случае задача (4), (5) оказывается нелинейной и для упругого поведения деформируемой среды.

2. Интегральные представления краевой задачи. Краевой задаче (4)—(7) могут быть поставлены в соответствие различные граничные интегральные уравнения, имеющие различные свойства погружения, и, как следствие, различные свойства используемых для их решения численных схем. Все они основаны на представлении решения v,j(xi,t) краевой задачи в виде поверхностного потенциала с заданной плотностью и свойствах непрерывности подобных потенциалов [1,2,12].

100

Метод граничных интегральных уравнений в моделировании ...

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

Cij (x)i j (x) = - Pij (x,y)i j (y)dsy + / Uij (x,y)pj (y)dsy+

JdD JdD

+ Uij(x,y)Xj(y)dv - ^jki(x,y)e°jkdv, (8)

JV JVp

где y £ dD — точка поверхности, x £ R3 \ dD — точка (источник) внутри области D; V, Vp — соответственно области с отличными от нуля объёмными силами и пластического деформирования и/или ненулевой разности температур; Cj(x) = bij, если x £ D, и Cij(x) = 0.5bj, если x £ dD, при условии её гладкости; Uij, Pij, 'Pjti — тензоры, соответствующие фундаментальному решению Кельвина—Сомильяны, определяемые выражениями

Uij = [16^(1 - v)r] 1 [(3 - Av)5ij + r,ir,j]

ij i j

i г dr

pj = - [8n(1 - v)т2]- ] — [(1 - 2v)5ij + 3r,ir,j]

3r,ir,j I -

j — j \gn^ ' j ' ’i

- (1 - )(r,inj - r,jm')},

Zjki = - [8n(1 - v)r2]-1 [(1 - 2v)(5ijr,k + bikr,j - bjkr,i) +

+ 3r,ir,j r,k],

dr

r = ||x - y\\ и r,i = — = r-1(yi - xi).

(9)

Существование и единственность решений уравнения (8) с условиями (5) для задач упругости установлено в классе непрерывных по Гельдеру функций C 1,a(D) С W2l(D), где W1(D) — пространство Соболева с параметром 0 < а ^ 1 [1,2,12]. Для нерегулярных областей, в том числе включающих поверхности разрыва (трещины), такие исследования выполнены, например в [13], на основе теории псевдодифференциальных операторов и метода Ви-нера—Хопфа.

Для рассматриваемых неоднородных (кусочно-однородных) сред уравнение (8) записывается отдельно для каждой из подобластей, входящих в область D и имеющих общие границы. В целом для области соответствующие уравнения объединяются с помощью условий (6), (7) на этих границах. Уравнение (8) применимо и для внешних краевых задач, когда область D дополняется до всего трехмерного пространства, если из его левой части вычесть интеграл с Pij по области дополнения.

Выражение для напряжений Ujk(xi,t) может быть получено дифференцированием (8) с учётом (1)—(3), и для точек внутри области D имеет вид

Cij — j iijk Pkds

Pkij 'Ck ds + I iijk Xk dV+ dD JdD JV

+ / j Upkl + abj) dV - 1 7 - 5VA 2^j - У + ^ bij ad, (10)

>Vv

15(1 - v)'

3(1 - v)

101

В. А. Петушков

где, принимая во внимание (9),

ukij £kij ^Sij uki,i, Pkij ^SijPki ,i.

Тензор £ijki определяется соотношением

£ijki = [3(1 - 2v)(5ijr,kr,i + 5kir,iT,j) +

+ 3v(5ur,jr,k + bjkr,ir,i + 6ikr,ir,j + Sjir,ir,k) - 15r,ir,jr,kr,i+

+ (1 - 2v)(5ik5ij + bjkбы) + (1 - 4v)6ijSki] (8n(1 - v)r3)-1.

Для вычисления напряжений на границе области D выражение (10) не применимо из-за сингулярного характера входящих в него интегралов. Как будет показано, напряжения в точках xi £ dD могут быть получены через граничные усилия и производные от граничных перемещений по касательным направлениям. Для этого необходима параметризация геометрии границы dD и аппроксимации решения на её поверхности. Аппроксимация решения также необходима и внутри области D для определения последних двух интегралов по внутренним объемам V, Vp в выражениях (8), (10) и параметров механики разрушения в окрестности трещины Qc(Xa).

Полагая рассматриваемые поверхности области D кусочно-гладкими, численное решение уравнения (8) выполним с использованием M граничных элементов Ak так, чтобы

м

dD « dDh = ^ Ak и Ak П Ai = 0, k = l.

k=1

Параметрическое представление геометрии и искомых функций Ui(xj,t) и Pi(xj,t) на каждом из них выберем подобно МКЭ в следующем виде:

Xi := Nr((j)x.

и „■

:= Nr (СУ )и

r,

Ph := Nr(ZyЖ, Zy £ Ak C R2

(11)

где Nr(Zy), r = 1, 2,..., 8 — квадратичная функция формы локальных координат Zy £ (-1,1), принадлежащая к Лагранжеву семейству конечных элементов Sh’m(dDh) в смысле [14].

Параметрическое представление (11) используется для вычисления поверхностных интегралов, входящих в (8), (10). Вычисление объёмных интегралов для массовых сил и по области начальных деформаций j, осуществляется с использованием объёмных 8 или 20-узловых конечных элементов Ak того же семейства с квадратичным изменением сил и деформаций, аналогично (11), но с функцией формы вида

Nr (Zy ,п) £ Skh,m(D), D « Dh

L

A

k=1

k,

П £ (-1,1) — третья локальная координата, нормальная к двум первым Zy, r = 1, 2,..., 8,..., 20. Отметим, что размер области с пластическими деформациями заранее не известен и устанавливается в процессе решения.

102

Метод граничных интегральных уравнений в моделировании ...

При наличии особенностей решения внутри и на границе области с известным характером сингулярности, её можно встроить в пространство пробных функций Sh’m(dDh) подобно МКЭ или за счёт подходящего построения неоднородной сетки dDh или расширения самого пространства путём введения специальных сингулярных функций, т. е. использования сингулярных элементов Afc [4,8].

Заменяя далее интегралы по границе и объёму соответствующими суммами по граничным и объёмным элементам, вместо (8), (10) получаем их численные аналоги:

[A]* = F+[D]{e0}, (12)

<т = -[A']X + F' + [D'}{s0} + [G]e°. (13)

Здесь [A], [A'] — матрицы коэффициентов при векторе граничных неизвестных; X — вектор, содержащий неизвестные усилия и перемещения на границе; векторы F, F содержат члены, определяемые отличными от нуля заданными значениями граничных перемещений и усилий, а также действием объёмных сил. Матрицы [D], [D'] включают в себя объёмные интегралы в (8), (10), квазидиагональная матрица [G] определяется наличием свободных от интеграла членов в (10). Для ускорения матрично-векторных перемножений могут быть использованы выше упомянутые быстрые методы [5].

Выражения для матриц и векторов в уравнения (12), (13) подробно представлены в [4], и здесь ввиду их громоздкости не приводятся. Вычисления входящих в них интегралов основано на численных квадратурах Гаусса. Объём и точность вычислений при этом зависят от порядка формулы и минимального расстояния между точкой наблюдения и узлом рассматриваемого элемента A к.

Если точка наблюдения совпадает с одним из узлов граничного или объёмного элементов, т. е. r ^ 0, интегралы в (8), (10) становятся сингулярными. Для их вычисления используется переход соответственно к полярным и сферическим координатам и квадратуры Гаусса по радиальным и угловым координатам.

Сильно сингулярные интегралы с особенностью O(r-3) рассматриваются в смысле главного значения Коши. Их вычисление по объёмным элементам строится путём перехода к сферическим координатам и выделением локальной зоны вблизи сингулярной точки, по которой проводится расчёт аналитически по радиусу и численно квадратурой Гаусса по угловым координатам. Тем самым получаем дополнительный вклад в диагональные блоки матрицы [D'] в выражении (13).

Для кусочно-однородных (составных) сред уравнения (12), (13) записываются отдельно для каждой из сред, и на общей поверхности между ними Г = Г1 U Г2 вводится два слоя узлов и элементов. В случае идеального контакта на первом из них в качестве граничных условий задаются неизвестные перемещения и нулевые усилия, на втором — неизвестные усилия и нулевые перемещения. В случае проскальзывания с трением на первом слое задаются неизвестные компоненты перемещений — нормальные и касательные, а на втором — не известные нормальные усилия и касательные перемещения. При этом условия (6), (7) учитываются при формировании матриц [A] и [A'].

103

В. А. Петушков

Исключив из (12), (13) вектор граничных неизвестных X, получаем окончательные выражения МГИУ для выбранного интегрального представления краевой задачи (4)-(7):

Х = [К]{ё° } + {m}, (14)

а = [В]{ё°} + {n} + [G]i°. (15)

Здесь

[K] = [A]-1[D], {m} = [A]-1F, [B] = [D'\ - [A'][K], {n} = F' - [A']{m}.

Отметим, что {m} и {n} представляют собой решение для упругой среды соответственно для вектора граничных неизвестных и вектора напряжений. Для определения обратной матрицы [A]-1 в её элиминативной форме и вектора {m} используется алгоритм гауссова типа подобно [15] с блочным исключением и выбором главного элемента в блоке. При этом нулевые блоки матрицы, появляющиеся при анализе составных сред, не формируются, не хранятся и не участвуют в процессе исключения, что позволяет решать прикладные задачи большой размерности.

Сложные истории нагружения деформируемой среды на временном слое Dt = (0, т) представим в виде совокупности простых нагрузок с помощью зависящих от времени параметров силового Лp(t) и температурного Л,(t) нагружений. Разрешающие уравнения (14), (15) в этом случае приобретают вид

k N

X\tk = SK]{^(tj)} + Л(>(tk)[K]a{0} + ^ XP(tk){mj}, (16)

j=1 j=1

\k

< = E([B ] + [G,]){ir(ti)} + Л, (tj )([B] + [a, ])a{9}+

j= 1

N

+ Fpj(tj){nj}, (17)

j=1

где t € Dt, {mj} и {nj} для любого шага j = 1,2 ,...,N полного нагружения определяются, как и в (14), (15). Такой подход позволяет рассматривать как повторяющиеся статические, так и циклические условия нагружения нелинейно деформируемой среды с анизотропным упрочнением, используемым в (3).

Моделирование нелинейного деформирования среды под действием j-того уровня нагружения осуществляется на основе классической схемы метода последовательной линеаризации в форме начальных деформаций [8,10]. Приращения пластической деформации определяются на каждой итерации из соотношений (9) с использованием реальных диаграмм деформирования a^ep^0, построенных для каждой из рассматриваемых сред.

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

104

Метод граничных интегральных уравнений в моделировании ...

Запишем выражения для напряжений и усилий через производные от перемещений, используя (1), (2) и первую формулу (5),

&jk — ui,i + p (Uj,k + Uk,j) , (18)

Pj — &jknk — ^&jknkUi,i + p (Uj,knk + Uk,jnk) .

Граничные усилия pj(xi,t), а также перемещения Uj(xi,t) и их производные по локальным координатам Д с учётом аппроксимации (11) могут быть получены из решения системы уравнений (12). Затем через них с использованием (18) определяются напряжения на границе dD рассматриваемой области D.

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

Эти особенности подхода делают целесообразным его обобщение на нелинейно деформируемые среды, расчётные области которых состоят из элементов с сильно различающимися размерами или включают в себя произвольно ориентированные дефекты — трещины Qc(Xa) с особым характером нагружения и взаимодействия их берегов.

В этом случае область, занимаемая деформируемой средой, рассматривается в виде дискретных подобластей, образованных соответствующими произвольно ориентированными фиктивными границами, например Г2 на рис. 1, с заданными на них условиями типа (6) или (7), если необходимо учитывать взаимодействия берегов трещины.

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

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

К преимуществам используемого в этом случае МГИУ следует отнести, прежде всего, возможность сколь угодно точного определения напряжённых и деформированных состояний в окрестности трещины, при наличии решения (14), (15) на поверхности рассматриваемого тела. Более того, анализ подрастания трещины в соответствии с историей нагружения также не вызывает затруднений.

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

105

В. А. Петушков

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

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

Фронт трещины

Рис. 2. Схема разрушения вдоль фронта трещины при объёмном напряжённом состоянии [Figure 2. Scheme of destruction along the crack front under triaxial stress conditions]

На границе трещины Qc(xa) искомые решения (14), (15) претерпевают разрыв, а вблизи её фронта имеют особенность типа r1/2 и r-1/2 соответственно для перемещений и усилий. Для учёта этих особенностей используется то же представление (11), но со специальным набором базисных функций Nr (Z7), как, например в [7,16].

КИН вычисляются непосредственно через перемещения

u1

2 r \1/2, .г [. . ю 3ю'

— ) (1 + v и K1 (5 - 8v) cos 77 - cos — +

п ' v 2 2 -

+ K2

. . ю 3ю 1

(9 — 8v) sin ^ + sin —

1

4Ё

U2

^\I/2(1 + v){K+7 — 8v)sin ю — sin + 1 +

п / 1 L 2 2 .

+ K2

ю 3ю1

(3 — 8v) cos ^ + cos —

1

4Ё

2(1 + v)(2r/n)1/2 K . ю u3 =---------e---------K3sin 2 ’

которые получены вблизи фронта трещины с помощью МГИУ. Здесь индексы 1, 2, 3 соответствуют направлениям I, II, III на рис. 2.

Полагая ф = п, как на рис. 2, получаем выражения для КИН, соответствующие трем типам упругого или квазиупругого (маломасштабное течение)

106

Метод граничных интегральных уравнений в моделировании ...

разрушения:

Ki = VU2(n/2r)1/2/(1 - V),

K2 = vui(n/2r)1/2/(1 - v), (19)

K3 = ^из(п/2г)1/2.

Наличие развитых зон пластического деформирования в окрестности трещин с размером, равным или большим характерного размера трещины, делает соотношения (19) не применимыми. В этом случае решение краевых задач нелинейной механики разрушения строится средствами МГИУ в соответствии с (8), (17). В качестве параметров разрушения, контролирующих поведение трещины, используются энергетические G- или J-интегралы и деформационные [16,17], при этом наиболее подходящим, имеющим ясный физический смысл и доступным измерению является ^-раскрытие трещины.

Подход, развиваемый здесь для описания составных кусочно-однородных ограниченных сред, может быть легко обобщён на решение трёхмерных задач механики разрушения с произвольно ориентированными в пространстве трещинами. Тело с трещиной рассматривается в этом случае в виде составного тела, внутренняя (фиктивная) граница которого включает в себя поверхность разрыва — трещину Qc(xa) с граничными условиями (6) и (10) с учётом нагруженности берегов трещины и их взаимодействия, как, например, в когезионных моделях разрушения [18].

5. Численные результаты. Важным вопросом при моделировании (и особенно для нелинейных процессов деформирования и разрушения сред) является математическое обоснование точности и вычислительной устойчивости получаемых результатов. Применительно к МГИУ такое обоснование может быть выполнено на основе интенсивно развиваемой теории псевдодифференциальных операторов [19].

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

Вначале приведём решение известной задачи Кирша в объёмной постановке с диаметром отверстия, равным толщине пластины. Размеры и история нагружения показаны на рис. 3 а. В силу симметрии рассматривается 1/8 часть полосы.

Материал полосы — аустенитная сталь (E = 2.0 ■ 105 МПа, V = 0.3, a0.2 = = 250 МПа, параметр упрочнения m = 0.25; a^a0.2 = (e?E/a0.2)m — диаграмма деформирования). Однородные растягивающие напряжения ao приложены к двум противоположным граням пластины. Максимальная величина прикладываемой нагрузки ao = 0.6ao.2. История нагружения пластины на Dt = (0, т) также приведена на рис. 3 а. Результаты моделирования представлены на рис. 3 и 4 для двух моментов времени, соответствующих максимальной нагрузке (рис. 3 а, b) и полной разгрузке (рис. 4).

В случае упругого поведения материала задача имеет аналитическое решение. Для нелинейного деформирования воспользуемся сравнительными результатами, полученными МКЭ [17]. На рис. 3 а результаты решения МГИУ (линии 1) показаны вместе с решением, полученным МКЭ (линии 2), и соот-

107

о

GO

Рис. 3. Объёмная задача Кирша: а) геометрия, зона пластического деформирования и распределение напряжений (1 — МГИУ, 2 — МКЭ); Ь) распределение деформаций

[Figure 3. 3D Kirsch’s problem: a) the geometry, the area of plastic deformation, and the stress distribution (1 — Boundary Integral Equation Method

(BIEM), 2 — Finite Element Method (FEM)); b) the strain distribution]

В. А. Петушков

Метод граничных интегральных уравнений в моделировании ...

Рис. 4. Распределение остаточных напряжений [Figure 4. Residual stresses distribution]

ветствуют максимальной нагрузке в истории нагружения. Выполненное моделирование процесса деформирования пластины позволило выявить особенности формирования вблизи отверстия зон нелинейных деформаций (рис. 3, b) и остаточных напряжений (рис. 4), где приведены распределения остаточных напряжений в номинальном сечении на срединной плоскости (рис. 4, а) и вдоль образующей контура отверстия по толщине пластины (рис. 4, b).

Аппроксимация поверхности и внутренней пластической области в пластине выполнена с использованием 20 поверхностных восьми узловых граничных элементов (62 узла) и 12 объёмных 20-узловых элементов. Для сравнения укажем, что в решении, полученном МКЭ, использовалось 26 объёмных 20-узловых элементов (900 уравнений) с тем, чтобы получить решение для коэффициента концентрации напряжений как на поверхности, так и на срединной плоскости вблизи отверстия пластины с погрешностью менее 3 % по сравнению с аналитическим.

Погрешность используемого здесь МГИУ,в том числе вместе с методом дискретных областей, не превышает 1 %. Хорошее совпадение результатов вблизи отверстия получено и для задачи пластичности. В остальной области более точными следует считать, по-видимому, результаты, полученные МГИУ. Увеличение числа граничных элементов вдвое не оказало заметного влияния на точность полученных результатов при 16-точечной (по 4 в каждом направлении) квадратуре Гаусса вычисления интегралов.

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

В качестве другого модельного примера, имеющего прикладное значение и полученные ранее результаты, рассмотрим объёмную задачу механики разрушения для трещины нормального разрыва полуэллиптической формы на поверхности полосы. Геометрия полосы и трещины вместе с расчётной схемой приведены на рис. 5; здесь B/H = 2; H/c = 2.5. Максимальная нагрузка и механические свойства материала такие же, что и в предыдущей задаче.

Рассматривались два случая полуэллиптической трещины с соотношени-

109

В. А. Петушков

ем полуосей а/c = 1/2 и более нагруженного с соотношением а/c = 2/3, полученные результаты для которых отмечены на рис. 5 цифрами 1 и 3 соответственно.

Аппроксимация геометрии полосы с трещиной выполнена с использованием 24 восьми узловых граничных элементов (74 узла) и десяти объемных 20-узловых элементов для последующего вычисления размеров предполагаемой зоны пластического деформирования в окрестности трещины. Для описания особенностей решения вблизи фронта трещины использовано 4 сингулярных граничных элемента с представлением [4].

Для трещины с размерами а/c = 0.5 и упругого поведения материала выполнен сравнительный анализ полученного распределения Ki вдоль фронта трещины (линия 1, рис. 5) с аналогичным распределением, заимствованным из [20] (линия 2, рис. 5). Столь же хорошее совпадение результатов с решением МКЭ [17] получено и для нелинейной задачи для случая а/c = 2/3. Максимальные значения величины КИН на графиках соответствуют условию плоской деформации, реализуемому в самой глубокой точке трещин. На поверхности пластины реализуется условие плоского напряженного состояния.

Существенное различие результатов в обоих случаях наблюдается лишь в достаточно тонком слое в окрестности трещины вблизи поверхности полосы при ф = 90°. Выявленное здесь повышение КИН проявляется обычно при экспериментальном изучении роста трещины в условиях циклического упругопластического деформирования, обнаружить которое другими методами, например МКЭ, оказывается затруднительно. Зона пластических деформаций в окрестности полуэллиптической трещины а/c = 2/3 оказалась весьма локализованной. Пространственный характер её расположения выделен цветом (заливкой) на рис. 5.

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

Полоса изготовлена из высоколегированной стали с антикоррозионной наплавкой из нержавеющей стали. Подобного рода плакированные соединения широко используются в химии и ядерной энергетике. Характеристики основного металла (материал 1 — реакторная сталь): Ei = 2.15 ■ 105 МПа, Vi =

Рис. 5. Полоса с поверхностной трещиной [Figure 5. The strip with a surface crack]

10° 30° 50° 70° tp

110

пиши

Метод граничных интегральных уравнений в моделировании ...

-1.0

-0.5

0.5 r/a

сЗ

Рч

I

о

b

Рис. 7. Распределение напряжений и деформаций в окрестности трещины [Figure 7. Stress and strain distribution in the vicinity of the crack]

111

В. А. Петушков

= 0.3, а^ = 430 МПа; характеристики нержавеющей наплавки (материал 2): Е2 = 2.05 ■ 105 МПа, v2 = 0.27, а2 = 180 MПa. Размеры полосы с дефектом (трещиной) приведены на рис. 6. Здесь, как и в предыдущем примере, B/H = 2, H/c = 2.5. Толщина плакированного слоя h = 0.1 H, размеры трещины а/c = 0.9, где а = 1.2h, приняты эквивалентными по площади дефекту, выявленному дефектоскопией в реальной конструкции.

Зоны пластических деформаций, полученные в окрестности трещины для двух случаев нагружения полосы ао = 0.3аУ и 0.4а1, показаны на рис. 6 и отмечены цифрами 1 и 2 соответственно. Основной металл (материал 1) при этих нагрузках работает в упругой области. В качестве параметра разрушения выбрано раскрытие трещины 5, его значения, вычисленные для указанных случаев нагружения, также приведены на рис. 6. Трещина начинает распространяться при условии 5 ^ 5c, где критическое значение 5С определяется экспериментально.

Представляющие практический интерес распределения напряжений и деформаций в окрестности трещины вдоль осей а и c, отражающие влияние неоднородности свойств плакированного соединения и пластического течения в наплавке, представлены на рис. 7.

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

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

СПИСОК ЛИТЕРАТУРЫ/ REFERENCES

1. С. Г. Михлин, Многомерные сингулярные интегралы и интегральные уравнения. М.: Физматгиз, 1962. 254 с.; S. G. Mikhlin, Multidimensional Singular Integrals and Integral Equations, International Series of Monographs in Pure and Applied Mathematics, vol. 83, New York, Pergamon Press, 1965, xii+259 pp.

2. В. Д. Купрадзе, Т. Г. Гегелия, М. О. Башелейшвили, Т. В. Бурчуладзе, Трехмерные задачи математической теории упругости и термоупругости. М.: Наука, 1986. 664 с.; V. D. Kupradze, T. G. Gegelia, M. O. Basheleishvili, T. V. Burchuladze, Three-Dimensional Problems of the Mathematical Theory of Elasticity and Thermoelasticity, New York, North-Holland, 1979, xix+929 pp.

3. K. Hayami, “Variable Transformations for Nearly Singular Integrals in the Boundary Element Method”, Publ. Res. Inst. Math. Sci., 2005, vol. 41, pp. 821-842 doi: 10.2977/ prims/1145474596.

4. В. А. Петушков, “Численная реализация метода граничных интегральных уравнений применительно к нелинейным задачам механики деформирования и разрушения объемных тел” // Моделирование в механике: Сб. научн. тр. ИТПМ СО АН СССР, Ново-

112

Метод граничных интегральных уравнений в моделировании ...

сибирск, 1989. Т. 3(30), №1. С. 133-156. [V. A. Petushkov, “Numerical realization of the method of boundary integral equations as applied to nonlinear problems of mechanics of deformation and fracture of volume bodies”, Model. Mekh., 1989, vol. 3(30), no. 1, pp. 133156 (In Russian)].

5. S. Rjasanow, O. Steinbach, The fast solution of boundary integral equations, Mathematical and Analytical Techniques with Applications to Engineering, Heidelberg, Springer, 2007, xi+279 pp. doi: 10.1007/0-387-34042-4.

6. Y. J. Liu, S. Mukherjee, N. Nishimura, M. Schanz, W. Ye, A. Sutradhar, E. Pan, N. A. Dumont, A. Frangi, A. Saez, “Recent Advances and Emerging Applications of the Boundary Element Method”, Appl. Mech. Rev., 2012, vol. 64, no. 3, 030802, 38 pp. doi: 10.1115/1.4005491.

7. N. A. Makhutov, V. A. Petushkov, V. I. Zysin, “Using the method of boundary integral equations for the numerical solution of volume problems of nonlinear fracture mechanics”, Strength of Materials, 1988, vol. 20, no. 7, pp. 833-841 doi: 10.1007/BF01528693.

8. В. А. Петушков, В. И. Зысин, “Пакет прикладных программ МЕГРЭ-3Д для численного моделирования нелинейных процессов деформирования и разрушения объемных тел. Алгоритмы и реализация в ОС ЕС” / Программное обеспечение математического моделирования, Сб. Пакеты прикладных программ. М.: Наука, 1992. С. 111-126. [V. A. Petushkov, V. I. Zysin, ““MEGRE-3D” software package for numerical simulation of nonlinear processes of deformation and fracture 3D bodies. Algorithms and implementation in OS ES”, Programmnoye obespecheniye matematicheskogo modelirovaniya [Software for Mathematical Simulation], Moscow, Nauka, 1992, pp. 111-126 (In Russian)].

9. G. C. Hsiao, O. Steinbach, W. L. Wendland, “Domain decomposition methods via boundary integral equations”, J. Comp. Appl. Math., 2000, vol. 125, no. 1-2, pp. 521-537 doi: 10. 1016/S0377-0427(00)00488-X.

10. В. А. Петушков, Р. М. Шнейдерович, “О термоупругопластическом деформировании гофрированных оболочек вращения при конечных смещениях” // Проблемы прочности, 1979. №6. С. 21-27; V. A. Petushkov, R. M. Shneiderovich, “Thermoelastic plastic deformation of corrugated shells of revolution at finite displacements”, Strength of Materials, 1979, vol. 11, no. 6, pp. 578-585 doi: 10.1007/BF00770100.

11. J. Necas, I. Hlavacek, Mathematical Theory of Elastic and Elasto-Plastic Bodies — An Introduction, Studies in Applied Mechanics, vol. 3, Amsterdam, New York, Elsevier Sci. Publ. Co., 1980, 342 pp. doi: 10.1016/B978-0-444-99754-8.50001-1.

12. W. L. Wendland, “On Some Mathematical Aspects of Boundary Element Methods for Elliptic Problems”, The mathematics of finite elements and applications. V. V, ed. J. R. Whiteman, New York, Acad. Press Inc., 1985, pp. 193-227 doi: 10.1016/ B978-0-12-747255-3.50019-8.

13. M. Costabel, “Boundary Integral Operators on Lipschitz Domains: Elementary Results”, SIAM J. Math. Anal., 1988, vol. 19, no. 3, pp. 613-626 doi: 10.1137/0519043.

14. G. Strang, G. Fix, An Analysis of the Finite Element Method, 2nd edition, Wellesley, MA, Wellesley-Cambridge Press, 2008, 400 pp.; Г. Стренг, Д. Фикс, Теория метода конечных элементов. М.: Мир, 1977. 294 с.

15. J. M. Crotty, “A block equation solver for large unsymmetric matrices arising in the boundary integral equation method”, Int. J. Numer. Meth. Engng., 1982, vol. 18, no. 7, pp. 997-1017 doi: 10.1002/nme.1620180705.

16. В. А. Петушков, А. Ф. Аникин, “Исследование разрушения трехмерных упругих тел с трещинами”// Прикладная механика, 1986. Т. 22, №9. С. 15-23; V. A. Petushkov, A. F. Anikin, “Investigation of the fracture of three-dimensional elastic bodies with cracks”, Soviet Applied Mechanics, 1986, vol. 22, no. 9, pp. 815-822 doi: 10.1007/BF00888886.

17. А. Ф. Аникин, В. А. Петушков, “О комплексе программ «САПР-82» и вычислительных аспектах моделирования на ЭВМ пространственных процессов деформирования и разрушения конструкций при повышенных температурах” // Проблемы прочности, 1987. №7. С. 62-67; A. F. Anikin, V. A. Petushkov, ““SAPR-82” software package and

113

В. А. Петушков

the computational aspects of the computer modeling of three-dimensional deformation and failure processes of structures at elevated temperatures”, Strength of Materials, 1987, vol. 19, no. 7, pp. 944-951 doi: 10.1007/BF01523534.

18. K. Park, G. H. Paulino, “Cohesive Zone Models: A Critical Review of Traction-Separation Relationships Across Fracture Surfaces”, Appl. Mech. Rev., 2013, vol. 64, no. 6, 060802, 20 pp. doi:10.1115/1.4023110.

19. G. C. Hsiao, W. L. Wendland, Boundary Integral Equations, Applied Mathematical Sciences, vol. 164, Berlin, Springer, 2008, xix+618 pp. doi: 10.1007/978-3-540-68545-6.

20. I. S. Raju, J. C. Newman Jr., “Stress-intensity factors for a wide range of semi-elliptical surface cracks in finite-thickness plates”, Engng. Fracture Mech., 1979, vol. 11, no. 4, pp. 817829 doi:10.1016/0013-7944(79)90139-5.

Поступила в редакцию 25/1/2014; в окончательном варианте — 17/IV/2014; принята в печать — 23/V/2014.

MSC: 74S30, 74R20

BOUNDARY INTEGRAL EQUATION METHOD IN THE MODELING OF NONLINEAR DEFORMATION AND FAILURE OF THE 3D INHOMOGENEOUS MEDIA

V. A. Petushkov

A. A. Blagonravov Mechanical Engineering Institute RAS,

4, M. Khariton’evskii per., Moscow, 101990, Russian Federation.

The method of boundary integral equations is applied, for solving the nonlinear problems of thermo-elastic-plastic deformation and fracture of inhomogeneous 3D bodies of the complex form. Solution is constructed on the basis of the generalized identity of Somigliana involving method of sequential linearization in the form of initial plastic deformations. The increments of plastic deformation are determined on the basis of the flow theory of hardening elastoplastic media with the use of modifed Prandtl-Reus’s relations. The cases of complex thermo mechanical loading of compound piecewise homogeneous media in contact are considered. For describing the processes of nonlinear deformation and fracture of the bodies with a complex geometry and local peculiarities a method of discrete domains (DDBIEM) is developed. The solutions of some practical significant 3D non-linear problems of the mechanics of deformation and fracture are presented.

Keywords: inhomogeneous 3D media, nonlinear deformation and fracture, boundary integral equation method, collocation approach, subdomains method.

Received 25/I/2014;

received in revised form 17/IV/2014;

accepted 23/V/2014.

ISSN: 2310-7081 (online), 1991-8615 (print); doi: http://dx.doi.org/10.14498/vsgtu1292 © 2014 Samara State Technical University.

Citation: V. A. Petushkov, “Boundary Integral Equation Method in the Modeling of Nonlinear Deformation and Failure of the 3D Inhomogeneous Media”, Vestn. Samar. Gos. Tekhn. Univ., Ser. Fiz.-Mat. Nauki [J. Samara State Tech. Univ., Ser. Phys. & Math. Sci.], 2014, no. 2(35), pp. 96-114. doi: 10.14498/vsgtu1292. (In Russian)

Author Details: Vladimir A. Petushkov (Dr. Phys. & Math. Sci.), Professor, Lab. of Mathematics Modeling.

E-mail address: pva_ras@inbox.com

114