CC BY
93
Серия «Математика»
2012. Т. 5, № 2. С. 69-80
Онлайн-доступ к журналу: http://isu.ru/izvestia
УДК 517.968.43, 532.135
Численное решение нелинейных интегральных уравнений Вольтерра с дробно-экспоненциальными ядрами реологических моделей вязкоупругой среды
А. Н. Тында
Пензенский государственный университет
А. Е. Романов
Самарский государственный университет
Аннотация. Статья посвящена численному исследованию реологических моделей в рамках теории нелинейной наследственной ползучести. Предложен численный метод решения нелинейных слабосингулярных интегральных уравнений Вольтерра с ядром Ржаницына, используемых в реологических моделях вязкоупругой среды. Сходимость метода проиллюстрирована на модельном уравнении, также приведены результаты моделирования ползучести на примере суглинка.
Ключевые слова: Нелинейные интегральные уравнения Вольтерра; ядро Ржаницына; вязкоупругость, метод Ньютона - Канторовича.
Известной проблемой теоретического исследования поведения системы контактирующих упругих и вязкоупругих тел является обоснование возможности получения адекватного решения в рамках математической модели контактного взаимодействия, построенной на основе интегрального соотношения между компонентами тензоров напряжения а и деформации е со сложными ядрами ползучести или релаксации. Реологические модели вязкоупругой среды на основе такого соотношения используются в постановке и решении задач контактной механики [1].
В теории наследственной ползучести известны различные виды ядер реологических моделей [7, 4], но далеко не для всех получены аналитические (численные, полуаналитические) решения интегральных уравнений связи напряжений, деформаций и времени Ь. Поэтому целью
1. Введение
статьи является построение метода решения интегральных уравнений с дробно-экспоненциальным ядром Ржаницына, используемым в ряде реологических моделей вязкоупругих сред [7, 9, 4].
2. Реологические уравнения в интегральной форме
Свойство вязкоупругости подразумевает, что, начиная с отличного от нуля момента времени, тело обладает ярко выраженными упругими свойствами на коротких промежутках времени и ярко выраженными вязкими свойствами на более длинных промежутках времени. На таких длинных промежутках времени Т вязкоупругие среды можно считать упругоползучими [4].
Основной закон нелинейного деформирования часто записывается либо в дифференциальной форме через реологические константы материала, либо в интегральной [7, 9]
Уравнение (2.1) является интегральным аналогом закона линейного деформирования, учитывающего наследственность деформаций тела. Здесь Е — мгновенный модуль упругости, первое слагаемое в левой части описывает изменение напряжения во времени за счет упругой деформации (мгновенно действующей по сравнению с длительностью Т), второе — за счет деформации ползучести, а правая часть характеризует нестационарный источник возникновения и развития напряжения в вязкоупругой среде. В (2.1) учитывается, что с течением времени свойства упругоползучих материалов способны меняться1, а также влияние наследственности, так как в начальный момент времени вязкоупругая среда может находиться в напряженно-деформированном состоянии (в результате осадки, уплотнения и т. п.). Поскольку наследственность проявляется в виде ранних изменений состояния тела, произошедших к моменту Ь, то в реологических моделях в результате декомпозиции интеграла в (2.1) на два, «наследуемое» состояние записывается в виде интеграла с пределами (—те, 0], а текущее — с пределами [0,Ь].
Подынтегральная функция f (а) в (2.1) обнаруживает свое влияние на формирование деформаций е возрастанием напряжений, но не ранее, чем с момента загрузки Ь = 0. Поэтому в нелинейной вязкоупругости согласно [9, 8, 6] применяются нелинейные зависимости. Например,
1 Материалы с такими свойствами относятся к категории «стареющих». Время играет роль возраста.
і
(2.1)
функция напряжения представима в степенном виде
/ (а) = Оа^, (2.2)
где О,у — эмпирические константы, /(0) = 0. При этом в (2.1) каждый из множителей произведения К/ может не иметь явного физического смысла, а являться чисто математической конструкцией. Например, в [8] при аппроксимации функции напряжения сдвига принималось 2 < ц < 5.
Простые наследственные функции некорректно описывают поведение тела в близкие к начальному моменты времени, если в эти моменты времени скорость деформации е! ^ те. С целью улучшения качества реологических моделей вязкоупругой среды используется слабосингулярное ядро Ржаницына [9, 4, 6]
К (Ь, т) = Ае—в(ь—т)(Ь — т)а—1. (2.3)
Оно описывает кривую скорости деформирования вязкоупругой среды при постоянной нагрузке единичной величины. Параметры ядра подбираются так, чтобы А, в > 0 и 0 < а < 1, и позволяют хорошо аппроксимировать множество полученных экспериментально кривых скорости деформирования [6]. Наличие в ядре К параметра в > 0 позволяет учитывать затухающую ползучесть, характерную, в частности, для водонасыщенных суглинков [4, 8].
Уравнение (2.1) согласно [9] может быть обобщено для «стареющих» вязкоупругих сред следующим образом:
т' N
а(Ь)+ ЕКп (Ь,т) /п [а(т)] йт = Е(Ь)е(Ь). (2.4)
^ п=1
—те
Поскольку при произвольных кривых деформирования решения
(2.1) и (2.4) не могут быть получены аналитически, то с целью снижения вычислительных трудностей используется ряд упрощающих допущений, индивидуальных для выбранной вязкоупругой среды. Например, для бетона при напряжениях, не превышающих половину предела прочности, предложено использовать в качестве (2.2) квадратный полином. В [8] зависимость обобщается до уровня /(а) = + фа^ при описании
реологии глинистых грунтов. В работе Чжао-Цзу-у согласно [9] упрощения вида (2.4) касаются не только функции напряжения, но и ядер: /1(а) = а, К2 (Ь,т) = 1, N = 2.
Резюмируя вышеизложенное, для уравнений (2.1) и (2.4) актуальной является разработка численных методов, позволяющих получить достаточно точные решения при произвольных кривых деформирования, модулях упругости, параметрах ядра К(Ь, т) и функции напряжения /(а).
Ниже предлагается эффективный численный метод решения уравнений вида (2.1) без учета унаследованного напряженно-деформированного состояния к моменту загрузки.
3. Исследование и решение основного уравнения
Рассмотрим уравнение
t
Ee(t) = a(t) + У K(t,r)f (a(r))dr, t e [0,T], (3.1)
0
где E — мгновенный модуль упругости, а ядро имеет следующий вид K(t, т) = Ae-e(t-T)(t — т)а-1, 0 < а < 1, A > 0, в > 0. (3.2)
3.1. ГЛАДКОСТЬ РЕШЕНИЙ
Уравнение (3.1) с ядром (3.2) представляют собой нелинейное слабосингулярное интегральное уравнение Вольтерра второго рода. Известно, что точное решение этого уравнения может иметь неограниченные производные в точке t = 0. Для того чтобы установить гладкостные свойства последнего дадим определение используемого класса функций [12]:
Определение 1. Функция g(t) принадлежит классу Cr,Y(t0,T] если имеет непрерывные производные до порядка r включительно при t e (t0,T], удовлетворяющие неравенствам
{1, при k < 1 — y,
1 + |log(t — to)\, при k = 1 — y, t e (to,T],
(t-to)\-l+Y , при k> 1 — ^(,
при k = 0,1, ...,r, y e (—<x>, 1).
Используя результаты работы [12], где рассмотрено слабосингулярное ядро общего вида K{Ь,т,а(т)), приходим к следующей теореме, устанавливающей существование, единственность и характер решения уравнения (3.1):
Теорема 1. Пусть функция f (а) имеет непрерывную производную при а e R+. Пусть также e(t) e Cr,l-a(0, T], 0 < а < 1. Тогда уравнение (3.1) с ядром (3.2) имеет единственное решение a(t) e Cr,l-a(0,T].
4. Численный метод
Существуют различные подходы к решению уравнений вида (3.1). В случае разностных ядер помимо коллокационных и итерационных методов применяются средства операционного исчисления. Однако численное вычисление с достаточной точностью обратного преобразования Лапласа во многих случаях является затруднительным. Для оптимальной по точности и сложности аппроксимации точного решения применяются различные неравномерные сетки. В работах H. Brunner (см., например, [11]) показано, что точность решения слабосингулярного уравнения методом коллокации на равномерных и квазиравномерных сетках не лучше O(N-а) для любого порядка гладкости r ^ 1 функций ядер и правой части. Поэтому для восстановления оптимального порядка
O(N-r) многие авторы используют существенно неравномерные сетки
r
(в лит. «graded mesh») вида tk = T (jj) a . Однако при этом метод колло-кации начинается на очень маленьком подинтервале, что может привести к значительным ошибкам округления и ограничивает применимость метода [14].
В работах [13, 10, 16] предложены некоторые альтернативные подходы к аппроксимации решений слабосингулярных уравнений. В [11, 12] для дискретизации слабосингулярных уравнений Вольтерра применяется метод сплайн-коллокации, приводящий к системам нелинейных алгебраических уравнений относительно неизвестных коэффициентов сплайна, а для решения таких систем используются итерационные алгоритмы.
В настоящей работе предлагается сначала линеаризировать интегральный оператор, а возникающие при этом линейные уравнения дискретизировать с помощью полиномиальной сплайн-аппроксимации. Подобный подход был применен в [15].
4.1. ЛИНЕАРИЗАЦИЯ
Воспользуемся модифицированным методом Ньютона - Канторовича. Предположим дополнительно, что функция f (а) имеет непрерывную вторую производную в R+ (данное условие не является ограничительным, так как в настоящей работе зависимость f от а имеет степенной характер). Определим далее нелинейный интегральный оператор
t
(Fa)(t) = a(t) — Ee(t) + A f в—тщ! т) 1-a . (4Л)
0
Модифицированный метод Ньютона - Канторовича, примененный к уравнению Fa = 0, примет вид
am+i = ат — [F'(ao)]-1(F(ат)), m = 0,1,...., (4.2)
где а0(Ь) — начальное приближение, а производная Г ,(ст0) оператора Г, вычисленная в точке сто, определяется следующим образом:
ґ = „ш Г(сто + Шст) - Гм =
ш—^0 Ш
і
Иш —
ш—о Ш
шмг) + а [ / Мг) + Ш*(т)\- {Ыт)|іт
7 Єв(і-т)(І - т)1-а
Применяя к подынтегральной функции формулу конечных приращений и учитывая, что f /(а0 + вша) — f /(а0) ^ 0 при ш ^ 0 равномерно относительно а, ||а|| = 1, 0 < в < 1, получаем
і
[Р(а) = ^ + А / в— Щ ^^-а йТ
0
Следовательно, для поправки Аа1 = а1(г) — а0(г) на первой итерации имеем уравнение
[Г/(ао)](Аа1) = —Г (ао). (4.3)
Продолжая этот процесс, получаем последовательность
приближенных решений, определяемых из уравнений
і
дт т , л / ГЫт)]Аат+1(т) Аат+1(і) + ^ ев(і-т)(і - т)і-а Лт =
о
і
_ (Ь) л [ f [ат (т)\(1т , Еє(ь)
-°т(І) - А] ев(і-т) и - т)1-а + Еє(І) •
(4.4)
где Аат+1(г) = ат+1(г) — ат(г), т = 0,1,....
Уравнение (4.4) линейное относительно функции Аат+1(г), определив которую, находим ат+1 (г).
Преобразовав (4.4), имеем окончательно
Ь
ат+1(г) + А! ^-^щ^+у-а ^ = фт(г), т = 0, 1, 2,..., (4.5)
0
где
Ь
f /[ао(т)]ат (т) — f [ат(т)]
Ф (Ь) Еє(Ь) , л І Ґітттт(т) - f [от
Фт(і) = Еє(і) + АУ ------е/3(і-т)(і __ т)1-а
ту,/ Чт.
4.2. Теоремы сходимости
На основании общих теорем о сходимости метода Ньютона - Канторовича [5] сформулируем следующую теорему сходимости итерационного процесса (4.5) в банаховом пространстве C[О, T] со стандартной
нормой ||ж(.)||с[0т] = max b(t)|.
’ te[0,т ]
Теорема 2. Пусть оператор F имеет непрерывную вторую производную в шаре Q0 (||a — aollcp,!1] ^ Р) и выполнены следующие условия:
1) уравнение (4-3) имеет единственное решение на [О, T], т. е. существует Г0 = [F/(a0)]-1;
2) ||Aa 1 |о[0,т] < П;
3) ||ГoF//(a) ||с[0,т] < L, a є Qo.
Если также
, r 1 1 — /1 — 2h 1 + y/1 — 2h
h = Ln К - и--------------n ^ p ^ -------------n,
2 h h
то уравнение (3.1) имеет единственное решение a* в Q0, процесс (4.5) сходится к a*, а скорость сходимости оценивается неравенством
11 am — ^Нср.т] ^ h (1 — Z1 — 2h)m+l, m = 0, 1,....
Существование единственного непрерывного решения каждого из уравнений (4.5), m = 0,1, 2,... (а уравнение (4.3) его частный случай при m = 0) является следствием теоремы 1, так как в ходе итераций мы не выходим за пределы класса Cr,1-a(0,T] С C[0, T]. Таким образом, условие 1 теоремы выполняется.
Нетрудно также проверить, что оператор F имеет в точке a0 є C[0, T] непрерывную вторую производную, которая является билинейным интегральным оператором следующего вида
t
\F//(a у|Г_ =л _ л ( f//[aoiт)]aiт)aiт) , lF (a0)\(a,a) = AJ ев(-т)(t — т) 1-а ГІт.
0
Для оценки скорости сходимости итерационного процесса (4.5) можно предложить альтернативную теорему, не требующую существования второй производной [F//(ao), что может быть полезно при использовании функции f (a) более общего вида, чем (2.2), рассматриваемого в данной задаче:
Теорема 3 ([2]). Пусть ||Aai||c[o,T] ^ Bq, F(хо)
C[0,T ]
= п0- Пусть
также в шаре Q jа : ||а — &0Цс[0,т] ^ (q < 1) выполняется сле
дующее условие
Fl(a(1)) — Fl(a(2)) ^ q/B0. Тогда уравнение (3.1)
с[о,т ]
имеет единственное решение а* в приближения (4-5) сходятся к а* и справедлива оценка погрешности
дтПоВо
1 — q '
4.3. Решение линейного уравнения
Применение описанного выше метода подразумевает необходимость решения на каждой итерации линейного слабосингулярного интегрального уравнения Вольтерра 2-го рода вида (4.5). Для получения приближенного решения будем использовать метод сплайн-коллокации с использованием существенно неравномерных сеток. При этом для минимизации ошибок округления на малых начальных интервалах применим подход, позволяющий уменьшить при практических вычислениях показатель неравномерности сетки посредством явного выделения особенности в аппроксимирующей функции. Подробно такой подход описан в работе [10]. В работах [3] и [16] доказана оптимальность по порядку по точности и сложности метода сплайн-коллокации для линейных слабосингулярных уравнений Вольтерра на различных классах функций.
5. Численные эксперименты
5.1. Сходимость на примере модельного уравнения
Проиллюстрируем сходимость предложенного алгоритма на следующем модельном уравнении
t
, . . f аа(т)+ba2(т) , , , 2. Г (a + 2bc)nt 4b з г п
a(t)+A —w dT = 2(ac + bc2)Vt—~-----------------—+—12, t e [0,1],
J V t — т 2 3
о
(5.1)
точным решением которого является функция a*(t) = c — л/t.
Ниже, в таблице 1, приведены результаты работы алгоритма. При этом используются следующие обозначения: N — число подобластей разбиения отрезка [0,Т], m — число итераций модифицированного метода Ньютона-Канторовича; r — порядок сплайна; £\ = max\am(t%k) —
k,i
a*(tk)\ — погрешность в узлах сетки; е2 = Цат(t) — a*(t)||C[0;T].
Таблица 1. Результаты решения уравнения (5.1).
т N г £1 £2
1 10 2 0.0078 0.0082
2 10 2 6.49 • 10~4 8.67 • 10~4
3 10 2 9.73 • 10~Б 2.84 • 10~4
5 10 2 4.33 • 10~б 7.61 • 10~б
2 20 2 6.29 • 10-Б 1.13 • 10~4
5 20 2 1.13 • 10~Б 6.32 • 10~б
2 20 3 8.26 • 10~6 2.92 • 10~Б
5 20 3 1.14 • 10~6 4.26 • 10~6
5 50 4 6.37 • 10-7 1.95 • 10~6
В расчетах параметры принимали следующие значения: А = 1, а = 0.015, Ь = 0.399, с = 1.
5.2. Решение уравнения (3.1) для водонасыщенного суглинка
Деформация одномерного уплотнения (компрессии) глинистых грунтов под действием постоянной нагрузки состоит из мгновенной компрессионной деформации ео и деформации компрессионной ползучести
е(г) = ео^) + еь(г), ег(0) = 0. (5.2)
Для примера из работы [8] возьмем кривые одноосного сжатия суглинка (с.185), полученные экспериментально в течение Т ~ 257 дней при различных постоянных нагрузках pj £ [0.025, 0.15] кг/мм2, ] = 1, 5. Начальная влажность равна 34.4% и находится за пределами зоны пластичности. Здесь модуль Е = 6.6 кг/мм2 является модулем упругости Юнга и определяет мгновенную компрессионную деформацию как
[ео Щ = Е. (5.3)
При полиномиальной аппроксимации функции /(а) с помощью метода наименьших квадратов получена следующая зависимость:
/(а) = 0.022а + 5.283а3 - 80.016а4'5 + 488.588а6. (5.4)
Максимальная погрешность такой аппроксимации достигается на начальном участке и составляет 1.5%, в среднем вдоль всей кривой — менее 0.01%.
Рис. 1. Динамика изменения напряжений суглинка при рх = 0.025 кг/мм2 0.03
0 Об
1-1
к
Н
I
I
И
0.04
0.02
\ 2
1—*
0.1 0.2 0.3 0.4 0.5 0.6
Время [б/р]
0.7
0.8
0.9
Рис. 2. Динамика изменения напряжений суглинка при р5 = 0.15 кг/мм2
На Рис.1,2 приведены результаты решения интегрального уравнения
(3.1) для суглинка при значениях р\ = 0.025 и ръ = 0.15 кг/мм2. При этом использовались следующие параметры ядра Ржаницына, соответствующие исходным кривым ползучести е^Ь) и еъ(Ь):
А\ = Аъ = 1, ^1 = 0.35, аъ = 0.04, в1 = въ = 2.
На указанных рисунках линия 1 описывает приближенное решение (Ь), а линия 2 — свободный член Ее^ (Ь) (] = 1 или ] = 5) уравнения (3.1). С целью уменьшения ошибок округления здесь используется нормированный временной интервал [0,1].
Следует отметить, что за время Т относительные изменения пористости, пропорциональные е(Ь) и равные 0.26% и 2.5% соответственно при р1 и ръ, свидетельствуют о том, что сжатие суглинка приводит к
незначительному изменению объема порового пространства. Из рис.1 можно сделать вывод, что при небольших нагрузках закон деформирования почти линеен (а(Ь) ~ Ее(Ь)), однако с ростом нагрузки (рис. 2) нелинейность процесса проявляется гораздо более существенно.
Среднее интегральное значение вклада ползучести на всем временном интервале при р = р1 возрастает медленнее, чем при р = р5, где ползучесть быстро возрастает на начальном этапе и составляет на всем интервале времени значение, близкое к максимальному. Таким образом, следует ожидать, что при нагрузках р > ръ, не приводящих за время Т к разрушению, вклад ползучести в деформацию суглинка будет еще более значительным и на длительных интервалах времени может превалировать над составляющей мгновенной деформации ео. Подобные сведения востребованы при скважинной нефтегазодобыче, когда продуктивный пласт контактирует с пластами, имеющими на интервалах времени эксплуатации скважины ярко выраженную ползучесть, усиливающуюся с глубиной.
Как видно из рисунков 1 и 2, найденное приближенное решение ат(Ь) уравнения (3.1) имеет особенности производных в точке Ь = 0, зависящие от параметра а ядра К и ведет себя как Ьа при Ь ^ 0, что подтверждает выводы теоремы 1.
Список литературы
1. Александров В. М. Аналитические методы в контактных задачах теории упругости / В. М. Александров, М. И. Чебаков. - М. : Физматлит, 2004. - 304 с.
2. Бойков И. В. Некоторые вопросы приближенного решения нелинейных операторных уравнений методом Ньютона - Канторовича / И. В. Бойков // Сб. аспирант. работ. Точные науки. - Казань : Изд-во Казан. гос. ун-та, 1970. - С. 82-94.
3. Бойков И. В. Оптимальные по точности приближенные методы решения интегральных уравнений Вольтерра / И. В. Бойков, А. Н. Тында // Дифференц. уравнения. - 2002. - Т. 38, N 9. - С. 1225-1232.
4. Вялов С. С. Реологические основы механики грунтов / С. С. Вялов. — М. : Высш. шк., 1978. - 447 с.
5. Канторович Л. В. Функциональный анализ / Л. В. Канторович, Г. П. Акилов.
- М. : Наука, 1979. - 744 с.
6. Колтунов М. А. Ползучесть и релаксация / М. А. Колтунов. - М. : Высш. шк., 1976. - 277 с.
7. Малкин А. Я. Реология: концепции, методы, приложения / А. Я. Малкин, А. И. Исаев. - СПб. : Профессия, 2007. - 560 с.
8. Месчян С. Р. Экспериментальная реология глинистых грунтов / С. Р. Месчян.
- М. : Недра, 1985. - 342 с.
9. Ржаницын А. Р. Теория ползучести / А. Р. Ржаницын. - М. : Стройиздат, 1968. - 418 с.
10. Тында А. Н. Смешанный сплайн-коллокационный метод решения слабосингулярных интегральных уравнений Вольтерра / А. Н. Тында // Тр. СВМО. -2005. - Т. 7, N 1. - С. 351-358.
11. Brunner H. Collocation methods for Volterra integral and related functional differential equations / H. Brunner. - Cambridge University Press, Cambridge, 2004.
12. Brunner H. The piecewise polynomial collocation method for nonlinear weakly singular Volterra equation / H. Brunner, A. Pedas, G. Vainikko // Math. Comp.
- 1999. - Vol. 68, N 227. — Р. 1079-1095.
13. Diogo T. Numerical solution of a nonlinear Abel type Volterra integral equation / T. Diogo, P. Lima, M. Rebelo // Commun. Pure Appl. Analysis. - 2006. - Vol. 5.
- P. 277-288.
14. Diogo T. Collocation methods for second-kind Volterra integral equations with weakly singular kernels / T. Diogo, S. McKee, T. Tang // Proc. Roy. Soc. Edin. -1994 - 124A. - P. 199-210.
15. Tynda A. N. On Numerical Solution of a Nonlinear Volterra Integral Equation of Fractional Order / A. N. Tynda // Тр. СВМО. - 2009. - Vol. 11, N 1. - P. 71-79.
16. Tynda A. N. Numerical algorithms of optimal complexity for weakly singular Volterra integral equations / A. N. Tynda // Computational Methods in Applied Mathematics. - 2006. - Vol. 6, N 4. - P. 436-442.
A. N. Tynda, A. E. Romanov
Numerical solution of nonlinear Volterra integral equations with fractionally-exponential kernels of rheological models of viscoelastic continuum
Abstract. This paper is devoted to numerical treatment of rheological models in the context of nonlinear heritable creep theory. An approximate method for nonlinear weakly singular Volterra integral equations with Rzhanitsyn’s kernel used in rheological models of viscoelastic continuum is suggested. In conclusion we adduce some numerical results demonstrating the convergence of this method and describing the deformation of loamy soil.
Keywords: Nonlinear Volterra integral equations, Rzhanitsyn’s kernel, viscoelasticity, Newton-Kantorovich method.
Тында Александр Николаевич, кандидат физико-математических наук, доцент, Пензенский государственный университет, 440026, Пенза, ул. Красная, 40, тел.: (8412)368277 (tynda@pnzgu.ru)
Романов Алексей Евгеньевич, кандидат физико-математических наук, Самарский государственный университет, 443011, Самара, ул. Академика Павлова, 1, тел.: (846)3345436 (rom-alex@mail.ru)
Tynda Aleksandr, Penza State University, 40, Krasnaya St., Penza, 440026, assistant professor, Phone: (8412)368277 (tynda@pnzgu.ru)
Romanov Alexey, Samara State University, 1, Academic Pavlov St., Samara, 443011, assistant professor, Phone: (846)3345436 (rom-alex@mail.ru)
