МОДЕЛИРОВАНИЕ ВЯЗКОУПРУГИХ ХАРАКТЕРИСТИК МАТЕРИАЛОВ НА ОСНОВЕ ЧИСЛЕННОГО ОБРАЩЕНИЯ ПРЕОБРАЗОВАНИЯ ЛАПЛАСА Текст научной статьи по специальности «Математика»

УДК 539.3

DOI: 10.18698/2309-3684-2020-3-321

Моделирование вязкоупругих характеристик материалов на основе численного обращения преобразования Лапласа

© А.А. Валишин, М.А. Тиняев МГТУ им. Н.Э. Баумана, Москва, 105005, Россия

При проектировании изделий из композиционных материалов, предназначенных для эксплуатации в сложных условиях неоднородных деформаций и температур, важно учитывать вязкоупругие свойства связующего и наполнителей. В статье анализируется взаимосвязь характеристик релаксации и ползучести. Рассмотрены все известные в литературе функции ползучести и релаксации. Подробно обсуждается проблема преобразования характеристик ползучести в характеристики релаксации и наоборот. Между функциями ползучести и релаксации существует простая взаимосвязь в пространстве изображений по Лапласу. Однако при возвращении в пространство оригиналов во многих случаях возникают большие трудности при обращении преобразования Лапласа. Рассмотрены два численных метода обращения преобразования Лапласа: использование ряда Фурье по синусам и метод квадратурных формул. Составлены алгоритмы и компьютерные программы для реализации этих методов. Показано, что время работы компьютерной программы, реализующей метод Фурье по синусам, почти в 2 раза меньше времени работы компьютерной программы, реализующей метод квадратурных формул. Однако первый метод уступает второму методу в погрешности вычислений: функции релаксации и скорости релаксации целесообразно находить первым способом, поскольку погрешность вычислений почти неразличима, а функции ползучести и скорости ползучести — вторым способом, т.к. для большинства таких функций результат, полученный вторым методом, значительно точнее результата, полученного первым методом. Функцию ползучести и функцию скорости ползучести Гаврильяка-Негами не удалось построить в связи со сложной рекурсивной формулой для коэффициентов ряда, однако, используя оба метода, эти функции всё же можно получить и сравнить друг с другом.

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

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

Ржаницына, функция Больцмана-Слонимского, функция Гаврильяка-Негами и др. [7-9]. Для инженерных расчетов важно уметь преобразовывать друг в друга ядра ползучести и релаксации. Между функциями ползучести и релаксации существует простая взаимосвязь в пространстве изображений по Лапласу. Однако при возвращении в пространство оригиналов во многих случаях возникают большие трудности при обращении преобразования Лапласа. В этих случаях приходится обращаться к численным методам обращения. Статья посвящена анализу, адаптации и применению различных численных методов обращения преобразования Лапласа и их программной реализации применительно к проблеме расчета вязкоупругих функций.

Математическая постановка задачи, принятые допущения. Для изотропных вязкоупругих материалов определяющие соотношения между напряжением и деформацией при одноосном напряженно-деформированном состоянии записывается так [1-11]:

ст(0 = 0(гМО) +10 (г - г) дГ Жт,

дт

(1)

е(г) = 3 №(0) + Г 3 (г - т) д^(Т) dт.

дт

где в (г) и 3 (г) — функции релаксации и ползучести, соответственно.

Функция в (г) описывает изменение напряжения со временем при постоянной деформации. Такой процесс называется релаксацией напряжения [12], в этом процессе напряжение убывает со временем, т.е. функция в (г) — убывающая и, следовательно,

ЖМ < 0. (2)

сИ

Функция 3 (г) описывает изменение деформации со временем при постоянном напряжении. Такой процесс называется ползучестью [12], в процессе ползучести при постоянном напряжении деформация растёт, т.е. функция 3 (г) — возрастающая и, следовательно,

> о. (3)

Л

В теории упругости между упругим модулем 0 и упругой податливостью 3 есть простая связь [13-17]

03 = 1. (4)

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

0

0

Запишем формулы (1) в пространстве изображений по Лапласу, используя теорему свёртки

a(p) = pG( p)s(p\ s( p) = pJ (p)a(p).

где s и а — тензоры деформаций и напряжений соответственно. Из формул (5) следует

p2 G (p) J (p) = 1. (6)

Это выражение определяет зависимость между изображениями функций релаксации и ползучести.

В теории преобразования Лапласа имеют место следующие предельные соотношения [2-3]

lim G(t) J (t) = 1,

t i0+ ir7\

lim G(t) J (t) = 1. ()

Это означает, что соотношения типа (4) в теории вязкоупругости имеют место только в двух предельных случаях: при t i 0 + и при t i ю.

Удобно функцию релаксации G(t) и функцию ползучести J (t) представить в безразмерном виде. Для этого обозначим

G(t) = GoW(t) [G0] = HA,

/2м (8)

J (t)=J0 g (t), [ J0]=м/Н.

Функции y/(t) и g(t) безразмерны, функцию y/(t) называется безразмерной функцией релаксации, а функция g (t) — безразмерной функцией ползучести.

В силу (4) и (7) имеем

G0J0 = 1

lim0+ W(t) g (t) = 1, (9)

t i0+

lim Kt) g (t) = 1.

В пространстве изображений по Лапласу между изображениями безразмерных функций y/(t) и g (t) имеет место соотношение типа (6), т.е.

pMp)g(p) = 1. (10)

Известные типы функций релаксации и ползучести. Приведем наиболее известные функции релаксации и ползучести [1-12].

1) Функции Максвелла. Функции Максвелла являются самыми простыми из всех известных. В таблице 1 приведены функции Максвелла и их изображения.

Таблица 1

Функции Максвелла и их изображения

№ Оригинал Изображение

1 ехр т(1 + рт)-1, тReр > -1

2 1 + */ /т (1 + рт) р-/

3 т 1 ехр [-(/ (1 + рт)"1

4 /-1 (т

2) Функции Абеля. В таблице 2 приведены функции Абеля и их изображения. Таблица 2 Функции Абеля и их изображения

№ Оригинал Изображение

1 1+х (-1)» [Г(3)]И (/Г р-1 [1+г(/)( рт)-3]-1

2 1 + 3- (/ р-1 [1 + Г(3)( рт)/

3 ,1X (-1)-+' [Г(3)Г [Г(»3)]-' (</)' —=1 {1+(рт)3[г(3)]-1)-1

4 Г(3)( рт)-3

В таблице 2 /3 принимает следующие значения 0 < /3 < 1. Заметим, что при / = 1 функции Абеля преобразуются в функции Максвелла.

3) Функции Работнова. В таблице 3 приведены функции Работнова и их изображения.

Таблица 3

Функции Работнова и их изображения

№ Оригинал Изображение

1 ! X (-1)» ( г/)7(—+" -=0 у(-+1)Г[г(-+ р-1 [1+(ю-7]-1

2 1+[Г(1+ 7)]-1 (ут)г р1 + (рт)7

3 г1х (-1)- Г[Г(-+1)]-1 (/т)г(п+г> —=0 1 + (рт)-7

4 г[Г(1+7)]-1т77-1 (рт)'7

В таблице 3 у принимает следующие значения 0 < у < 1. Заметим, что при у = 1 функции Работнова преобразуются в функции Максвелла.

4) Функции Ржаницына. В таблице 4 приведены функции Ржаницына и их изображения.

Таблица 4

Функции Ржаницына и их изображения

№ Оригинал Изображение

1 1-у/, Х)[Г(3)1 ^ р-1 [1 - (1 + />т)-3]

2 1 + £у( п/, /^[ПпЗ)]-1 п=1 р-1 [1 - (1 + рт)-3]-1

3 [Г(3)]-1т-3/3-1ехр [-(ут)] (1 + рту3

4 ' 1 ехр [-(>/)]]С[Г(п3)]-1 (%)" п=1 [(1 + рт)33-1]-1

В таблице 4 /3 принимает значения 0 < /3 < 1, а

У

у(/, Т) = . Заметим, что при ¡3 = 1 функции Ржаницына

0

преобразуются в функции Максвелла.

5) Функции Кольрауша. В таблице 5 приведены функции Кольрауша и их изображения.

Таблица 5

Функции Кольрауша и их изображения

№ Оригинал Изображение

1 ехр[-(Ут) \ р-1 {1+(1 (-1)п Г[(13П 1 п=1(п -1)! (рт)(1-3)и }

2 1+(1 -3)-1 X [ к ](/Г)п У И) пГ[(1 -Р)п]\/тТ да р 1 +т£ Ъп (рт)<3-1)"-1 п=1

3 (1 -3)т^3 ехр [-( /У] (1 3)Х (-!)" Г[(!-3)п] 1 р;х(п -1)! (рт)(1-3)и

4 т'-1'-3+т-1 X [ К ](УГ)п-1 Г[(1 -Р)п]\/т> да X к (^т)(/-1)и п=1

В таблице 5 / принимает значения 0 </< 1, Ь0 = 1,

^ к1 Г[(1 -3)к ] Ъп = У (-1)к-1^——1 Ъп к , п = 1,2,...,да . Заметим, что при / = 0 к=1 (к -1)!

функции Кольрауша преобразуются в функции Максвелла.

6) Функции Гаврильяка-Негами. В таблице 6 приведены функции Гаврильяка-Негами и их изображения.

Таблица 6

Функции Гаврильяка-Негами и их изображения

№ Оригинал Изображение

1 1 Е (-1)п г(п+3) /,/Ха(п+3) Е а(п + Р)п! Г[а(п + /)] Г//г> р-1 {1 -[1 + (ргТ]Р]

2 Б 1 + /\а(пР+т) 1+ЕЕ(-1)т Г[ ( / ) 1](г) П=1 т=0 Г[а(п/ + т) + 1\^Т> ,-1 (1 + [(,Г)а3-1]1 £(р^ )

3 Г1 ( ;/)а3Е (-1)П Г(п + 3) ( ,/)аП Чг> П=0 П! Г[а(п +3)]Г(/Л/г' [1 + (рг)аУ

4 Б 1 + /\а(пР+т) ЕЕ(-1)т Г[ ( 3+ )] (Г) П=1 т=0 Г[а(п3+ т)]\/Т> 1 Б [(1+рг)аа-1] Е<-1)т (Г-

В таблице 6а, /3 принимают значения 0 < а < 1,0 < /3 < 1,

т к1 к„-2 п-1 1 Г(п + 3)

Бт = ЕЕ......ЕП " k]+l), а(п) = ^^^П/З1, К = т , К = 0.

к1 =0 к2 =0 кп-1 ]=о п! 1 (Р)

Заметим, что при а = 1 функции Гаврильяка-Негами преобразуются в функции Ржаницына, а при а = /3 = 1 они преобразуются в функции Максвелла.

Во всех приведённых выше таблицах под номером 1 обозначена функция релаксации, под номером 2 — функция ползучести, под номером 3 — функция скорости релаксации и под номером 4 — функция скорости ползучести.

Отсутствие размерности у всех приведённых выше функций обусловлено параметром г, т.е. временем релаксации. Временем релаксации называется период времени, за который амплитудное значение возмущения в выведенной из равновесия физической системе уменьшается в е раз (е — основание натурального логарифма) [11-18].

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

методы обращения преобразования Лапласа. Далее рассмотрены два таких метода. [19]

Использование ряда Фурье по синусам. Обращение преобразования Лапласа с помощью ряда Фурье по синусам заключается в нахождении оригинала f (t) по значениям изображения F (p) в равноотстоящих точках на действительной оси. Этот способ основан на двух допущениях:

1) изображение F (p) существует и Rep > 0 (к этому всегда можно прийти, если рассматривать вместо изображения F (p) изображение F (p + а) при достаточно больших а е! );

2) f (0) = 0 (этого можно достичь, положив f1(t) = f (t)-f (0) exp(-t), что равносильно замене

F(p) ^ F(p - f <0)

(p +1)

Интеграл Лапласа

си

F (p) = J f (t)e-ptdt (11)

преобразуем, используя подстановку

е-а = cos 0, (12)

где а > 0 — произвольное действительное число. Из формулы (12) выразим г через 0 :

г = -—lncos0. (13)

а

Тогда

f(t) = f (-—IncosQ) = p(Q), (14)

о

где 0 <Q< л/2.

Интеграл (11) примет следующий вид:

л

1 } p-1

F(p) = —J (cosQ)— sin Qcp(Q)dQ. (15)

0 0

Функцию (р{в) разложим в ряд Фурье по синусам нечётных кратных дуг

да

ср(в) = Х Ck sin [(2k +1)0], (16)

k=0

0

где коэффициенты ck определяются обычным способом:

ck =-

4_

ж

sin [(2k + 1)в]в. (17)

Поскольку функция <р(в) неизвестна, а именно её необходимо найти, то поступим следующим образом: выразим значения коэффициентов ск через значения изображения F (р). В выражении (15) положим р = (2п + (п = 0,1,...), тогда получим

oF [(2n + 1)o] = J cos2n 6 sin 6cp(6)d6. (18)

0

Функция этого интеграла может быть представлена в виде линейной комбинации функций sin[(2k + 1)6]

1 n

cos2n в sin в = с\п - Ck;1) sin [[ 2(n - k) +1] в\, (19)

2 k=0

причём С2П = 0.

Выражения (16) и (19) подставим в формулу (18). В силу ортогональности тригонометрической системы функций,

ж

2

J sin [( 2/ + 1)в] sin [( 2v + 1)в\ de-

0, u^v,

ж (20)

-, u = v,

где Z.

Тогда при фиксированном n останутся только члены, для которых v = n - к ( к = 0,1,..., n ), т.е.

1 _ n

iск ск-Г

22n 4

oF [(2n + 1)o] = 4 £( C2n - С-1 )cn-k, (21)

Распишем ряд в (21) и, приведя подобные слагаемые, получим новое выражение

п 2к +1 4п+1

Е ^ С^с, = [ (2п + 1)а]. (22)

к=0 2п +1 -

Подставив в это равенство последовательно п = 0,1,..., получим неоднородную линейную относительно ск систему уравнений с треугольной матрицей для определения коэффициентов ск :

ж

0

0

4

Со = —oF (о-); к

42

Со + с =~°F (3о); (23)

43

2с0 + 3с1 + с2 = —oF (5о); к

Найденные ck подставляем в ряд Фурье (16) и из формулы (13) выражаем в через t :

в = arccos(eo). (24)

Таким образом, подставляя формулу (24) в выражение (16), получаем

ад

f (t) = 2 ck sin [(2k +1) arccos(e-ot)]. (25)

k=0

В зависимости от величины промежутка, для которого необходимо вычислить значение оригинала f (t) , о следует выбирать малым для

больших t и, наоборот, большим для малых t.

Метод квадратурных формул с равными коэффициентами.

Если задано изображение F (p), то соответствующий ему оригинал можно записать в виде интеграла Меллина:

1 с+гад

f (t) = — f eptF{p)dp. (26)

2m J

с-гад

По теореме о конечном значении

lim F(р) = 0, (27)

Re р^ад

поэтому можно предположить, что функция F (р) представима в виде

F (р) = А < р), (28)

Р

где s > 0, а функция <р{'р) регулярна в полуплоскости Re р > 0 и ограничена в полуплоскости Re р > с (с > 0). Тогда интеграл (25) запишем как

1 с+гад

f (t) = — Г е"р <Шр. (29)

2кг J

с-гад

В формуле (29) сделаем замену переменной р = р' /1:

/С) = — ^1 | ер'(р') ^(р'^р\ (30)

2жг

где 8 = е1 > 0 . Для удобства заменим р' на р и получим

л ь+гы

/С) = — ^ | е'р~Х»ф. (31)

-Г 2жг

Таким образом, задача обращения преобразования Лапласа сведена к вычислению интеграла

8+г'вд

-1 б^/1^

I = — | ерр-Х»ф. (32)

2жг J

Для его вычисления построим квадратурную формулу с равными коэффициентами

I - С ± ср(рк). (33)

к=1

Неизвестные величины Сп и рк (к = 0,1,...,п) в формуле (33) можно выбрать так, чтобы она была точной для любого многочлена степени п от переменной^р , т.е. (р(р) = 1/рк (к = 1,0,...,п). Множитель Сп определим так, чтобы выражение (33) было точным для функции

<Р( р 1:

I = — | е'р-'ф = пСп, (34)

тогда

11 Т 1

Сп = - — | ^-ф = — , (35)

П 8-г'» пГ(5)

где Г( ^ ) — гамма-функция.

С учётом (35), перепишем (34)

1 п

1 - "Т7Т ). (36)

пГ(5) к=1

Если перейти от переменной р к переменной х = 1/р и обозначить хк = 1/ рк то для неизвестных хк получится следующая система уравнений

х + х2 +... + х = пГ(5) Г ерр "р 1ёр = пГ(л) ; 1 2 " 2т Р Р Р Г(" +1)'

х2 + х22 +... + х2 = ^М Т еРр"р2 dp 1 2 " 2т/ * Р Р Р Г(" + 2У (37)

< + х2" +... + хпп = = - ^

2т Л Г(" + П)

Очевидно, что эта система нелинейна, поэтому её решение может вызвать вычислительные затруднения. Целесообразно найти другой способ нахождения хк. Введём полином (п (х) степени п, корнями которого будут числа хк:

(п (х) = (х - х1)(х - х2) . (х - хп ). (38)

Разложим (3 8) по степеням х :

( (х) = х" + А хп-1 + А2 хп-2 +... + Лп -1 х + Лп. (3 9)

Известно, что коэффициенты Лк будут элементарными симметрическими функциями корней (функцию /(х1,х2,...,хп) от п переменных называют симметрической, если она инвариантна относительно любой перестановки входящих в неё переменных х). Например, для

п = ^^ЛС^ х2) = -( х1 + х2 )= Л1 (х2 , X1), Л2 (х1, х2) = х1 х2 = Л2 (х2 , х1 ) .

Равенства (37) есть сумма степеней корней хк:

$ (/ =1п). (40)

к=1 Г(/ + ")

Из теории многочленов известны соотношения между элементарными симметрическими функциями корней Лк и функциями $. Эти соотношения называют соотношениями Ньютона, которые имеют следующий вид:

$1 + Л = 0;

$2 + ЛД + 2 Л2 = 0;

$ + Л1$2 + Л2 $ + 3 Л3 = 0; (41)

$п + ЛА-1 + Л2 $п-2 +... + пЛп = 0.

Подставив в систему (41) значения из (40), получим следующую систему уравнений:

1 4= о.

г(л +1) пг(^)

-+—1— а1 +—2— а = о,

Г(5 + 2) Г(5 +1) 1 пГ(5)

11 13

- +---А +---А А3 = 0; (42)

Г(5 + 3) Г(^ + 2) 1 Г(5 +1) 2 пГ(у)

1 1 1 п „ п

-- А1 + -^ АП-1 Ап = 0

Г(5 + п) Г(5 + п - 1) Г(5 + 1) пГ(,5)

Она линейна относительно Ак , причём её матрица является треугольной, определитель которой равен к !/[пк[Г(5)]к]. Найдя все Ак, построим полином а>п (х) и найдём его корни хк, а значит, определим все узлы хк = 1/рк квадратурной формулы (33). Однако для нахождения оригинала /(/) в формуле (36) необходимо заменить рк на рк^.

Наконец, формулы (31), с учётом вышеизложенного, запишем в окончательном виде:

/(/) V (р{^1 (43)

пГ(5)£Г I / ) К ]

Примеры численного решения задачи. Были написаны компьютерные программы, реализующие описанные выше методы. Приведем несколько расчетов, полученных с помощью данных программ. На рис. 1, 3, 5, 7 показаны графики, полученные для 1-го метода, на рис. 2, 4, 6, 8 — для 2-го метода. Во всех приведённых случаях параметры т = 10, а = 0,01, 3 = 0,5, к = 2,0, 5 = 1,0.

На рис. 1 и 2 приведены результаты работы программ, определяющих функцию ползучести Абеля, которую сравнивают с её аналитическим выражением g(/) = 1 + ¡т~3 .

На рис. 3 и 4 приведены результаты работы программ, определяющих скорость ползучести Работнова, которую сравнивают с её аналитическим выражением g1(/) = ¡[Г(1 + ¡)] 1 .

На рис. 5 и 6 приведены результаты работы программ, определяющих скорость релаксации Ржаницына, которую сравнивают с её аналитическим выражением у/1(/) = [Г(3)] 1 3-1 ехр [-(//т)] .

1

Моделирование вязкоупругих характеристик материалов на основе...

ЛО

7 6 5 4 3

2

0

20

40

60

100

Рис. 1. Функция ползучести Абеля при п = 20, Ag и 0,041

— численное решение, — — аналитическое решение

f С ) 6

5

4

3 2

0

Рис. 2. Функция ползучести Абеля при п = 100, Agcp и 0,0062 — — численное решение, — — аналитическое решение

200

Рис. 3. Функция скорости ползучести Работнова при п = 20, Ag

— — численное решение, — — аналитическое решение

1ср

: 0,0102

t

t

t

Рис. 4. Функция скорости ползучести Работнова при п = 100, Ag1 и 0,0034

— — численное решение, — аналитическое решение

Рис. 5. Функция скорости релаксации Ржаницына при п = 20, А^. и 0,0041

— численное решение, — — аналитическое решение

Рис. 6. Функция скорости релаксации Ржаницына при п = 100, А у/1ср и 0,0047

— — численное решение, — — аналитическое решение

На рис. 7 и 8 приведены результаты работы программ, определяющих функцию релаксации Кольрауша, которую сравнивают с её

аналитическим выражением) = ехр -(г/г)1 -

- г 100

Рис. 7. Функция релаксации Кольрауша при п = 20, Д ц/ср и 0,0001

— — численное решение, — — аналитическое решение

0,8 0,6 0,4 0,2

100

Рис. 8. Функция релаксации Кольрауша при п = 100, Д ц/ср и 0,0027 — — численное решение, — — аналитическое решение

Выводы. Проанализировав графики, делаем следующие выводы:

1) время работы компьютерной программы, реализующей метод Фурье по синусам (I метод), почти в 2 раза меньше времени работы компьютерной программы, реализующей метод квадратурных формул с равными коэффициентами (II метод);

2) однако I метод уступает II методу в погрешности вычислений: функции релаксации и скорости релаксации целесообразно находить I способом, поскольку погрешность вычислений почти неразличима, а функции ползучести и скорости ползучести — II способом, т.к. для большинства таких функций результат, полученный II методом, значительно точнее результата, полученного I методом;

г

3) функцию ползучести и функцию скорости ползучести Гаврильяка-Негами, определённые в таблице 6, не удалось построить в связи со сложной рекурсивной формулой для коэффициентов Dm ,

однако, используя оба метода, эти функции всё же можно получить и сравнить друг с другом.

Для всех случаев число слагаемых ряда Фурье n = 20 . Этот выбор обусловлен тем, что при росте n в таких случаях (кроме ядер Максвелла) значения коэффициентов ck быстро растут с ростом k, из-за чего ряд Фурье начнёт очень сильно осциллировать. Поэтому I метод следует применять для более простых случаев.

ЛИТЕРАТУРА

[1] Димитриенко Ю.И. Нелинейная механика сплошной среды. Москва, Физма-тлит, 2009, 624 с.

[2] Кристенсен Р. Введение в теорию вязкоупругости. Москва, Мир, 1974, 339 с.

[3] Бленд Д. Теория линейной вязкоупругости. Москва, Мир, 1965, 199 с.

[4] Биргер И.А., Пановко Я.Г. Прочность, устойчивость, колебания. Справочник в трёх томах, том 1. Москва, Машиностроение, 1968, 831 с.

[5] Димитриенко Ю.И. Механика сплошной среды. Т. 4. Основы механики твердых сред. Москва, Изд-во МГТУ им. Н.Э. Баумана, 2013, 624 с.

[6] Ильюшин А.А., Победря Б.Е. Основы математической теории термовязкоупругости. Москва, Наука, 1970, 356 с.

[7] Работнов Ю.Н. Ползучесть элементов конструкций. Москва, Наука, 1966, 752 с.

[8] Колтунов М.А. Ползучесть и релаксация. Москва, Высшая школа, 1976, 276 с.

[9] Москвитин Б.В. Сопротивление вязкоупругих материалов. Москва, Наука, 1972, 328 с.

[10] Ржаницын А.Р. Теория ползучести. Москва, Стройиздат, 1968, 416 с.

[11] Бартенев Г.М., Френкель С.Я. Физика полимеров. Ленинград, Химия, 1990, 432 с.

[12] Прохоров А.М. Физическая энциклопедия. Т.1. Аароново — Длинные. Москва, Советская энциклопедия, 1988, 703 с.

[13] Димитриенко Ю.И., Губарева Е.А., Сборщиков С.В. Конечно-элементное моделирование эффективных вязкоупругих свойств однонаправленных композиционных материалов. Математическое моделирование и численные методы, 2014, № 2, с. 28-48.

[14] Dimitrienko Yu.I., Minin V.V., Syzdykov E.K. Modeling of the thermomechani-cal processes in composite shells in local radiation heating. Composites: Mechanics, Computations, Applications, 2011, vol. 2, iss. 2, pp. 147-169.

[15] Dimitrienko Yu.I., Gubareva E.A., Pichugina A.E. Theory of composite cylindrical shells under quasistatic vibrations, based on an asymptotic analysis of the general viscoelasticity theory equations. IOP Conference Series: Material Science and Engineering, 2019, vol. 683, no. 012013. DOI: 10.1088/1757-899X/683/1/012013

[16] Димитриенко Ю. И., Коряков М. Н., Захаров А. А., Строганов А. С. Численное моделирование сопряженных аэрогазодинамических и термомеханических процессов в композитных конструкциях высокоскоростных летательных аппаратов. Математическое моделирование и численные методы, 2014, № 3, c. 3-24.

[17] Янке Е., Эмде Ф., Леш Ф. Специальные функции. Москва, Физматгиз, 1968, 344 с.

[18] Пагурова В.И. Таблицы неполной гамма-функции. Москва, Физматгиз, 1963, 235 с.

[19] Крылов В.И., Скобля Н.С. Методы приближенного преобразования Фурье и обращения преобразования Лапласа. Москва, Наука, 1974, 223 с.

Статья поступила в редакцию 14.10.2020

Ссылку на эту статью просим оформлять следующим образом: Валишин А.А., Тиняев М.А. Моделирование вязкоупругих характеристик материалов на основе численного обращения преобразования Лапласа. Математическое моделирование и численные методы, 2020, № 3, с. 3-21.

Валишин Анатолий Анатольевич — д-р. физ.-мат. наук, профессор кафедры «Вычислительная математика и математическая физика» МГТУ им. Н.Э. Баумана. e-mail: enf@mail.ru

Тиняев Михаил Александрович — магистр кафедры «Вычислительная математика и математическая физика» МГТУ им. Н.Э. Баумана. e-mail: infomaniak97@yandex.ru

Modeling the viscoelastic characteristics of materials based on the numerical inversion of the Laplace transform

© A.A. Valishin, M.A. Tinyaev Bauman Moscow State Technical University, Moscow, 105005, Russia

When designing products made of composite materials intended for operation in difficult conditions of inhomogeneous deformations and temperatures, it is important to take into account the viscoelastic properties of the binder and fillers. The article analyzes the relationship between relaxation and creep characteristics. All known creep and re-laxation kernels in the literature are considered. The problem of transformation of creep characteristics into relaxation characteristics and vice versa is discussed in de-tail. There is a simple relationship between the creep and relaxation functions in the Laplace image space. However, when returning to the space of the originals, in many cases there are great difficulties in reversing the Laplace transform. Two numerical methods for inverting the Laplace transform are considered: the use of the Fourier series in sine and the method of quadrature formulas. Algorithms and computer programs for realization of these methods are made. It is shown that the operating time of a computer program implementing the Fourier method by sine is almost 2 times less than the operating time of a computer program implementing the quadrature formula method. However, the first method is inferior to the latter method in accuracy of calculations: the relaxation functions and relaxation rates, it is advisable to find the first method, since the computational error is almost indiscernible, and the functions of creep and creep speed, the second way, because for most functions, the result obtained by the second method is much more accurate than the result obtained by the first method. The Gavrillac-Negami creep function and the Gavrillac-Negami creep rate function could not be constructed due to a complex recursive formula for the series coefficients, but using both methods, these functions can still be obtained and compared with each other.

Keywords: viscoelasticity, relaxation, creep, relaxation and creep nuclei, Laplace transform, Fourier method, quadrature formula method.

REFERENCES

[1] Dimitrienko Yu.I. Nelineynaya mekhanika sploshnoy sredy [Nonlinear continuum mechanics]. Moscow, Fizmatlit Publ., 2009, 624 p.

[2] Kristensen R. Vvedenie v teoriyu vyazkouprugosti sredy [sredy [Nonlinear continuum mechanics]. Moscow, Mir Publ., 1974, 339 p.

[3] Blend D. Teoriya linejnoj vyazkouprugosti [Theory of linear viscoelasticity]. Moscow, Mir Publ., 1965, 199 p.

[4] Birger I.A., Panovko Ya.G. Prochnost', ustojchivost', kolebaniya. Spravochnik v tryoh tomah, tom 1 [Strength, stability, oscillations. Handbook in three volumes, volume 1]. Moscow, Mashinostroenie Publ., 1968, 831 p.

[5] Dimitrienko Yu. I. Mekhanika sploshnoy sredy. Tom 4. Osnovy mekhaniki tver-dogo tela [Continuum Mechanics. Vol. 4. Fundamentals of solid mechanics]. Moscow, BMSTU Publ., 2013, 624 p.

[6] Il'yushin A.A., Pobedrya B.E. Osnovy matematicheskoj teorii termovyazkoupru-gosti [Fundamentals of the mathematical theory of thermovyazcoelasticity]. Moscow, Nauka Publ., 1970, 356 p.

[7] Rabotnov Yu.N. Polzuchest' elementov konstrukcij [Creep of structural elements]. Moscow, Nauka Publ., 1966, 752 p.

[8] Koltunov M.A. Polzuchest' i relaksaciya [Creep and relaxation]. Moscow, Vysshaya shkola Publ., 1976, 276 p.

[9] Moskvitin B.V. Soprotivlenie vyazkouprugih materialov [Resistance of viscoelas-tic materials]. Moscow, Nauka Publ., 1972, 328 p.

[10] Rzhanicyn A.R. Teoriyapolzuchesti [Theory of creep]. Moscow, Stroyizdat Publ., 1968, 416 p.

[11] Bartenev G.M., Frenkel S.Ya. Fizikapolimerov [Physics of polymers]. Leningrad, Chemistry Publ., 1990, 432 p.

[12] Prohorov A.M. Fizicheskaya enciklopediya. T.1. Aaronovo — Dlinnye [Physical Encyclopedia. Vol. 1. Aaronovo — Dlinnye]. Moscow, Soviet Encyclopedia Publ., 1988, 703 p.

[13] Dimitrienko Y.I., Gubareva E.A., Sborschikov S.V. Finite element modulation of effective viscoelastic properties of unilateral composite materials. Mathematical Modeling and Computational Methods, 2014, no. 2, pp. 28-48.

[14] Dimitrienko Yu.I., Minin V.V., Syzdykov E.K. Modeling of the thermomechani-cal processes in composite shells in local radiation heating. Composites: Mechanics, Computations, Applications, 2011, vol. 2, iss. 2, pp. 147-169.

[15] Dimitrienko Yu.I., Gubareva E.A., Pichugina A.E. Theory of composite cylindrical shells under quasistatic vibrations, based on an asymptotic analysis of the general viscoelasticity theory equations. IOP Conference Series: Material Science and Engineering, 2019, vol. 683, no. 012013. DOI: 10.1088/1757-899X/683/1/012013

[16] Dimitrienko Y.I., Koryakov M.N., Zakharov A.A., Stroganov A.S. Computational modeling of conjugated gasdynamic and thermomechanical processes in composite structures of high speed aircraft. Mathematical Modeling and Computational Methods, 2014, no. 3, pp. 3-24.

[17] Janke E., Emde F., Lesh F. Special'nye funkcii [Special functions]. Moscow, Fiz-matgiz Publ., 1968, 344 p.

[18] Pagurova V.I. Tablicy nepolnoj gamma-funkcii [Tables of incomplete gamma function]. Moscow, Fizmatgiz Publ., 1963, 235 p.

[19] Krylov V.I., Skoblya N.S. Metody priblizhennogo preobrazovaniya Fur'e i obrashcheniya preobrazovaniya Laplasa [Methods of the approximate Fourier transform and inversion of the Laplace transform]. Moscow, Nauka Publ., 1974, 223 p.

Valishin A.A., Dr. Sc. (Phys.-Math.), Professor, Department of Computational Mathematics and Mathematical Physics, Bauman Moscow State Technical University. e-mail: enf@mail.ru

Tinyaev M.A., Student of Department of Computational Mathematics and Mathematical Physics, Bauman Moscow State Technical University. e-mail: infomaniak97@yandex.ru