Mathematical modeling of hereditary elastically deformable body on the basis of structural models and fractional integrodifferentiation Riemann-Liouville apparatus

Full Text

Введение. Основной задачей математического моделирования реологических сред является установление функциональной связи между напряжениями и деформациями, возникающими в изучаемом объекте. В одномерном случае эта связь может быть задана посредством некоторого непрерывного, но необязательно линейного, оператора A следующим образом: (1) σ(t) = Aε(t), где σ(t) и ε(t) - напряжение и деформация тела в момент времени t. При заданной зависимости деформации от времени равенство (1) должно позволять находить зависимость напряжения как функции времени, в частности при ε = ε0 = const находить решение задачи о релаксации напряжений в вязкоупругом теле; предполагая обратимость оператора A, определять зависимость деформации от времени при заданном законе изменения напряжения, в частности, при σ(t) = σ0 = const находить решение задачи ползучести [1-3]. В феноменологических моделях вид, структура и свойства оператора A принимаются за исходный, первичный элемент и в конечном итоге устанавливаются экспериментальным путём. Другой подход к моделированию заключается в замене реального деформируемого твердого тела некоторой механической структурной моделью, представляющей собой определённые комбинации из идеальных упругих и вязких элементов, подчиняющихся соответственно законам Гука и Ньютона для вязкой жидкости. Известно, что в рамках структурного моделирования стандартную одномерную обобщённую модель вязкоупругого тела можно записать в виде n m bk Dk σ = Eε + σ(t) + k=1 ak Dk ε, (2) k=1 где Dk = (d/dt)k ; ak , bk , E0 - постоянные величины, обычно заранее не известные; n = m или m = n + 1. Частными случаями модели (2) являются собственно законы Гука σ(t) = E0 ε(t), и Ньютона σ(t) = η ε(t), ˙ а также следующие модели: 168 Математическое моделирование наследственно упругого деформируемого тела . . . - модель Фойхта σ(t) = E0 ε(t) + η ε(t), ˙ (3) E0 σ(t) + η σ(t) = ηE0 ε(t), ˙ ˙ (4) (E1 + E2 ) σ(t) + η σ(t) = E1 E2 ε(t) + ηE2 ε(t), ˙ ˙ (5) - модель Максвелла - модель Кельвина - модель Зенера E1 σ(t) + η σ(t) = E1 E2 ε(t) + η(E1 + E2 )ε(t), ˙ ˙ (6) где E0 , E1 , E2 - модули соответствующих упругих элементов, η - коэффициент демпфирования. Заметим, что всё многообразие реологических моделей можно разбить на два типа: модели, описывающие явление мгновенной упругой деформации в момент приложения нагрузки (модели типа Максвелла, Кельвина, Зенера) и модели типа Фойхта, не наделённые (не описывающие) мгновенной упругой деформацией. Очевидно, что решения задачи о ползучести при постоянно действующей нагрузке σ(t) = σ0 = const, получающиеся из определяющего соотношения (2), будут носить в основном экспоненциальный характер. Однако известно, что многие среды (полимеры, резины и т.п.) в испытаниях на ползучесть демонстрируют степенную зависимость деформации от времени. Оказалось, что подобное поведение характерно не только для механических сред, ряда физических явлений и биологических процессов, но вообще для любых динамических систем с памятью. Это обстоятельство дало толчок применению понятий и методов дробного анализа в задачах моделирования вязкоупругого поведения механических сред и наследственно упругого деформируемого твердого тела [4-7]. В. Вольтерра предложил рассматривать законы наследственной упругости, связывающие между собой напряжение σ и деформацию ε, в форме равенства [8], [9, с. 223] ε= 1 σ+ E t K(t - τ )σ(τ )dτ . (7) -∞ В многочисленных исследованиях по наследственной механике, берущих начало с работ Х. Больцмана [10] и продолженных в работах Г. Дюффинга [11], А. Джемента [12], А. П. Бронского [13], Г. Л. Слонимского [14], А. Ю. Ишлинского [15], А. Р. Ржаницына [16], Ю. Н. Работнова [17] и других авторов (см. библиографический список в [18-21]), показано и обосновано, что для целого ряда физических сред ядра ползучести в интегральном операторе (7) являются ядрами абелевского типа, а основные соотношения между напряжениями и деформациями выражаются через дробные интегралы Римана-Лиувилля [22]. Производные по времени дробных интегралов будут являться дробными производными Римана-Лиувилля. На это обстоятельство мы ещё не раз обратим внимание. Отметим, что, по-видимому, А. Джемент [23] впервые использовал при моделировании механических свойств вязкоупругих материалов дифференцирование дробного порядка. 169 О г о р о д н и к о в Е. Н., Р а д ч е н к о В. П., У н г а р о в а Л. Г. Отталкиваясь от экспериментальных результатов Дюффинга [11], в работе [17] Ю. Н. Работнов предложил модель вязкоупругой среды на основе интегрального оператора с ядром Абеля в следующей форме: Iα = (t - τ )α , Γ(1 + α) (8) α > -1. Интегральные операторы с ядрами вида (8) ˜ Iα ϕ = t 1 Γ(1 + α) ϕ(τ )dτ , (t - τ )α 0 t > 0, обладают следующими композиционными свойствами [4, с. 28], [17]: ˜ ˜ ˜ Iα Iβ = Iα+β+1 , ˜n ˜ Iα = Inα+n-1 . Разрешая интегральное уравнение Вольтерры второго рода с операто˜ ˜ ром Iα путём обращения оператора I - λIα (I - тождественный оператор), ˜ Ю. Н. Работнов определил резольвентный оператор Эα (λ), ядро которого определяется рядом ∞ Эα (λ, z) = z α (λz 1+α )n , Γ (n + 1)(1 + α) n=0 α > -1, и получил определяющее соотношение t σ(t) = E ε(t) - λ Эα (-λ, t - τ )ε(τ )dτ , 0 позволяющее вычислять напряжение σ(t), если известен закон изменения деформации ε(t). Если вместо (8) записывать ядра интегральных операторов в виде K(t - τ ) = (t - τ )α-1 , Γ(α) α > 0, (9) то соответствующий интегральный оператор α I0t ϕ = 1 Γ(α) t 0 ϕ(τ )dτ , (t - τ )1-α t > 0, (10) представляет собой левосторонний интеграл Римана-Лиувилля порядка α, который обладает естественным полугрупповым свойством [20, с. 73], [22, с. 42] β α α+β α β I0t I0t ϕ = I0t I0t ϕ = I0t ϕ (11) для любой суммируемой функции ϕ ∈ L(0, T ), T > 0. Решение интегрального уравнения Вольтерры второго рода с оператором (10) α ϕ(t) - λI0t ϕ(t) = f (t) 170 (12) Математическое моделирование наследственно упругого деформируемого тела . . . записывается в терминах резольвентного оператора µ,α E0t;λ f = t Exp(α, µ; λ; t - τ )f (τ )dτ, (13) 0 где Exp(α, µ; λ; t) = tµ-1 Eα (λtα ; µ) - обобщённая (двупараметрическая) дробная экспоненциальная функция [24], Eα (z; µ) - функция типа Миттаг-Леффлера [25]. Свойства оператора (13) и функции Exp(α, µ; λ; z) рассматривались во многих публикациях авторов настоящей работы [см., например, 26- 28]. Решение интегрального уравнения (12) записывается в виде α ϕ(t) = I - λI0t -1 α,α f = I + λE0t;λ f. К настоящему моменту времени опубликовано значительное количество работ, посвящённых изучению дробных реологических моделей (см. обзор литературы в [21], а также [29-31]), причём в большинстве источников одномерная обобщённая модель вязкоупругого тела с памятью (наследственно упругого тела) записывается в виде [32, с. 202], [33, с. 742] n m bk Dβk σ = E0 ε (t) + σ (t) + k=1 ak Dαk ε, (14) k=1 причём, как отмечается в монографии [20, с. 443], многие авторы используют в равенстве (14) «подходящие дробные производные», лишь бы их преобразование Фурье удовлетворяло условию ˆ F Dα f (t) (s) = (is)α f (is), ˆ где s - частота, f (is) = F (f (t)) (s) - преобразование Фурье функции f . Однако этому условию, вообще говоря, могут удовлетворять дробные производные и Грюнвальда-Летникова, и Лиувилля, и Капуто, и Римана-Лиувилля на оси или полуоси, если f (t) ≡ 0 при t < 0. Для нахождения явных решений дифференциальных уравнений дробного порядка, возникающих из определяющего соотношения (14), используется в основном преобразование Лапласа, в то время как они могут быть найдены собственно методами дробного анализа, что не требует обязательного использования интегральных преобразований. В качестве «удобных для анализа и приложений» дробных производных некоторые авторы используют лиувиллевскую форму дробного интегро-дифференцирования благодаря замечательному свойству этих операторов сохранять экспоненциальную функцию [34, с. 2] (см., также [20, с. 81], [22, с. 87] и [19, с. 311]). Известно, что для любого α > 0 α I+ eατ (t) = α D+ eατ (t) = 1 Γ(α) t -∞ 1 d Γ(n - α) dt n eατ dτ = λ-α eλt , (t - τ )1-α t -∞ eατ dτ = λα eλt , (t - τ )α-n+1 171 О г о р о д н и к о в Е. Н., Р а д ч е н к о В. П., У н г а р о в а Л. Г. где n = [α] + 1, [α] - целая часть числа α ∈ R. Отметим, что использование лиувиллевской формы дробного интегро-дифференцирования в математических моделях динамических систем с памятью действительно оправдано, если за исходный постулат брать соотношение (7), в котором K(x - t) - ядро Абеля. На самом деле «история деформирования» тела начинается с некоторого конечного момента времени, например, времени первого нагружения, которое принимают за начало отсчёта времени t = 0. Таким образом, в основном соотношении (7) возникает интегральный оператор с началом в точке t = 0. В частности, если в (7) ядро интегрального оператора пропорционально абелевому ядру (9) с некоторым коэффициентом β > 0, то соотношение (7) будет представлять собой дробный аналог модельного определяющего соотношения Максвелла в интегральной форме α Eε(t) = σ(t) + βI0t σ(t). (15) Применим к левой и правой частям равенства (15) оператор дифференα цирования D0t и воспользуемся тождеством [20, с. 74], [22, с. 50] α α D0t I0t ϕ(t) = ϕ(t) (α > 0), справедливым для любой функции ϕ(t) ∈ L(0, T ) (T > 0). Получим α α ED0t ε(t) = D0t σ(t) + βσ(t). (16) Очевидно, что в классе абсолютно непрерывных функций AC[0, T ] при α ∈ (0, 1) равенства (15) и (16) эквивалентны. Действительно, если α ∈ (0, 1), а функция ϕ(t) ∈ AC[0, T ], то α α I0t D0t ε(t) = ϕ(t). (17) Тогда, применяя к левой и правой частям равенства (16) оператор интегриα рования I0t , в силу тождества (17) получаем (15). Отметим здесь, что тождество (17) при α 1 несправедливо даже в классе AC n [0, T ]. Известно, что для функций f (t), принадлежащих классу n-α Lα (0, T ) = {f (t) ∈ L(0, T ) : I0t f (t) ∈ AC n [0, T ]}, α который обеспечивает существование суммируемой дробной производной D0t f , справедлив аналог формулы Ньютона-Лейбница [20, с. 75], [22, с. 50]: n α α I0t D0t f = f (t) - k=1 где (n-k) (n-k) fn-α (0+) α-k t , Γ(α - k + 1) α-k fn-α (0+) = lim D0t f, t→0+ k = 1, 2, . . . , n. Упомянутые выше классы абсолютно непрерывных функций, без которых невозможно построение полноценной теории дробного интегро-дифференцирования, основанного на лебеговской конструкции интеграла, можно заменить несколько более узкими классами функций, разрывных в конечном 172 Математическое моделирование наследственно упругого деформируемого тела . . . множестве точек, достаточное число раз дифференцируемых в интервалах своей непрерывности с суммируемой старшей производной. Действительно, в задачах одноосной ползучести особый интерес представляют функции напряжения, имеющие разрыв 1-го рода в конечном множестве точек, которые мы будем считать непрерывными слева в точках разрыва. Например, функция σ(t) = σ0 (t) H(t) - H(t - t1 ) (18) описывает нагружение для t : 0 < t t1 и разгрузку для t : t1 < t T деформируемого тела. Здесь σ0 (t) - некоторая неотрицательная непрерывно дифференцируемая на [0, T ] функция, σ0 (0) = σ0 > 0, H(t) = 0, t 0, 1, t > 0 - функция Хевисайда. Из (18) ясно, что σ(0) = 0, в то время как σ(0+) = lim σ(t) = σ0 (0) = σ0 . t→0+ Отметим, что ε(t) ∈ C(0, t1 ] ∪ C(t1 , T ], так как σ(t) ∈ C(0, t1 ] ∪ C(t1 , T ], а σ(t) ∈ C(0, t1 ) ∪ C(t1 , T ) ∩ L(0, T ). Кроме этого, из уравнения (15) следует ˙ равенство Eε(0+) = σ(0+) = σ0 , (19) которое описывает мгновенную упругую деформацию в момент приложения нагрузки. Вернёмся ещё раз к уравнению (16). Относительно функции u(t) = Eε(t) - σ(t) это уравнение является простейшим дифференциальным уравнением с дробной производной Римана-Лиувилля порядка α ∈ (0, 1): α D0t u = βσ(t). (20) Известно, что естественным начальным условием для такого дифференциального уравнения будет начальное условие типа Коши: lim I 1-α u t→0+ 0t = u0 . (21) В классе кусочно-непрерывных ограниченных функций начальное условие (21) будет однородным, что обеспечивает возможность применения тождества (17) для решения дифференциального уравнения (16). Возникает вопрос: каким может быть начальное условие для деформации ε(0), если (16) рассматривать как дифференциальное уравнение относительно деформации? Этот вопрос равносилен следующему: существует ли кусочно-дифференцируемое решение u(t) дифференциального уравнения (20) с начальным условием u(0) = u0 , если σ(t) - кусочно-непрерывная ограниченная функция? 173 О г о р о д н и к о в Е. Н., Р а д ч е н к о В. П., У н г а р о в а Л. Г. Хорошо известно [20, с. 73], [22, с. 46], что для функции u(t) ∈ AC[0, T ] справедлив достаточный признак существования дробной производной Римана-Лиувилля порядка α ∈ (0, 1): α D0t u(t) = u(0) -α 1-α t + I0t u(t). ˙ Γ(1 - α) (22) Это тождество тем более справедливо в классе кусочно-дифференцируα емых ограниченных функций. Из (22) видно, что производная D0t u имеет α интегрируемую особенность в нуле. Подставляя выражение производной D0t из (22) в левую часть уравнения (20), заключаем, что равенство (20) возможно, если u(0) = Eε(0) - σ(0) = 0. Отсюда следует, что начальное условие ε(0) для искомой функции деформации ε(t) не может быть произвольным, а зависит от значения нагрузки σ(0) в начальный момент времени: (23) ε(0) = σ(0)/E, что физически проявляется в форме скачка мгновенной упругой деформации. Формальная редукция дифференциального уравнения (16) или (20) к интегральному уравнению (15) очевидна. Действительно, запишем (20) с учетом (22) и (23). Применим к левой и α правой частям полученного равенства оператор интегрирования I0t : α α 1-α ˙ I0t I0t u(t) = βI0t σ(t). Используя полугрупповое свойство дробных интегралов (11) и равенство (23), находим 1 α 1-α ˙ ˙ I0t I0t u(t) = I0t u(t) = u(t) - u(0) = u(t) и окончательно получим равенство (15). 1-α Заметим, что возникающий в равенстве (22) оператор I0t u является ˙ α дробной производной Капуто (см., например, [5, 20, 21]) и обозначается C D0t u. В подавляющем большинстве публикаций по дробным реологическим моделям авторы используют в основном определяющем соотношении (2) именно дробную производную Капуто (см., например, [35, 36]), которую определяют как «нормализованную производную Римана- Лиувилля» равенством n-1 C α D0t f = α D0t f - k=0 f (k) (0) α tk-α = D0t f (t) - Γ(k - α + 1) n-1 k=0 f (k) (0) k t , (24) k! где n = [α]+1. Затем путём использования определения дробной производной Римана-Лиувилля из (24) выводится достаточный признак её существования в классе функций AC n [0, T ] и обнаруживается равенство C 174 n-α α D0t f = I0t f (n) . Математическое моделирование наследственно упругого деформируемого тела . . . В действительности же известно, что для функций f (t) ∈ AC n [0, T ] дробная производная Римана-Лиувилля порядка α > 0, n - 1 α < n, представима в виде [22, с. 46] n-1 α D0t f = k=0 откуда следует, что f k (0) n-α tk-α + I0t f (n) , Γ(k - α + 1) n-α α I0t f (n) = C D0t f, причем для функции f (t) при f (k) (0) = 0 (k = 0, 1, . . . , n - 1) дробные производные Римана-Лиувилля и Капуто совпадают. 1. От структурных механических моделей - к математическим моделям наследственно упругого тела. В качестве исходного постулата в настоящей работе взята феноменологическая модель памяти В. Вольтерры. Применительно к задачам о вязкоупругом поведении деформируемого твёрдого тела эта математическая модель сводится к равенству (15), позволяющему сразу найти закон изменения деформации ε(t), если зависимость напряжения σ(t) от времени задана. Во введении показано, что определяющее соотношение (15) эквивалентно соотношению (16) в дифференциальной форме для достаточно гладких функций с условием (19) или (23). Будем считать, что параметр β в уравнении (15) при α = 1 совпадает с величиной E/η, где E - модуль упругости, η - коэффициент демпфирования в классической модели Максвелла. Тогда видно, что если равенство (16) записывать в виде α α EηD0t ε(t) = ηD0t σ(t) + Eσ(t), (25) то при α → 1 - 0 оно переходит в классическое определяющее соотношение (4) и поэтому названо во введении дробным аналогом модели Максвелла. Запишем равенство (25) в виде η α α (26) ηD0t ε(t) = D0t σ(t) + σ(t) E и устремим E → +∞ в (26). Это равносильно замене упругого элемента в механической модели абсолютно твёрдым стержнем. В итоге приходим к определяющему соотношению α σ(t) = ηD0t ε(t), 0 < α < 1, (27) предложенному в 1947 году Скотт-Блэром [37, 38]. А. Н. Герасимов в 1948 году [39] предложил аналогичное (27) соотношение с использованием производной α α D+ (t > 0) на всей числовой оси [20, 22] вместо производной D0t на отрезке. Г. Л. Слонимский [14] рассматривал закон, устанавливающий связь между напряжением и деформацией в интегральной форме ε(t) = 1 -α 1 α D σ(t) = I0t σ(t). η 0t η (28) Ясно, что при условии ε(0) = 0 законы Скотт-Блэра и Слонимского эквивалентны. Они оба будут эквивалентны соотношению Герасимова, если проα изводная D+ рассматривается на функциях, равных нулю тождественно при t < 0. 175 О г о р о д н и к о в Е. Н., Р а д ч е н к о В. П., У н г а р о в а Л. Г. Определяющие соотношения (27) или (28) могут быть взяты в качестве простейших математических моделей вязкоупругого демпфирования (наследственно-упругого тела). Более сложные, зависящие от трёх и более параметров модельные уравнения могут получаться путём структурного моделирования, рассматривая различные комбинации и соединения идеально упругих элементов с элементом Скотт-Блэра. Действительно, параллельное соединение упругого элемента с модулем E0 и элемента Скотт-Блэра (27) приводит к дробному аналогу модели Фойхта: α σ(t) = E0 ε(t) + ηD0t ε(t). (29) Соединяя последовательно дробный аналог элемента Фойхта с упругим модулем E1 с идеальным упругим элементом, модуль которого равен E2 , приходим к дробному аналогу определяющего соотношения Кельвина: α α (E1 + E2 )σ(t) + ηD0t σ(t) = E1 E2 ε(t) + ηE2 D0t ε(t). (30) Если взять соотношение (25) с E = E1 , описывающее дробный аналог тела Максвелла, и присоединить параллельно ещё один упругий элемент с модулем E2 , придём к определяющему соотношению, которое ассоциируется с дробным аналогом модели (тела) Зенера: α α E1 σ(t) + ηD0t σ(t) = E1 E2 ε(t) + η(E1 + E2 )D0t ε(t). (31) Мы сохраняем в равенствах (27)-(31) обозначения постоянных коэффициентов в том же виде, что и в классических определяющих соотношениях (3)-(6), полагая, что при α = 1 они должны совпадать. 2. Аналитические решения задачи о ползучести для некоторых дробных аналогов классических реологических моделей. Пусть α ∈ (0, 1), σ = σ(t) - заданная суммируемая функция. 2.1. Модель Скотт-Блэра-Герасимова. Пусть η - постоянная величина. Определяющее соотношение в данной модели в соответствии с (27) записывается в виде 1 α D0t ε = σ(t). (32) η Ясно, что при α → 0 соотношение (32) переходит в закон Гука с заменой константы η на E0 , а при α → 1 - в закон Ньютона. Найдем решение этого простейшего дифференциального уравнения с начальным условием lim ε(t) = ε(0+) = 0 t→0+ в классе функций ε(t), таких, что ε(t) ∈ C(0, T ) ∩ L(0, T ). Учитывая тожде˙ α ство (17), применим к левой и правой частям уравнения (32) оператор I0t : ε(t) = 1 α I σ(t). η 0t Если σ(t) изменяется по закону (18), в котором σ0 (t) = σ0 > 0: σ(t) = σ0 H(t) - H(t - t1 ) , 176 (33) Математическое моделирование наследственно упругого деформируемого тела . . . то σ0 α I H(t) - H(t - t1 ) = η 0t t t t1 σ0 (t - τ )α-1 dτ = + H(t) (t - τ )α-1 dτ - H(t - t1 ) = ηΓ(α) t1 0 0 σ0 = tα H(t) - (t - t1 )α H(t - t1 ) . (34) ηΓ(α + 1) ε(t) = Отметим, что при α → 0 решение (34) сводится к закону Гука, если положить η = E0 , а при α → 1 определяет линейный закон деформирования. 2.2. Дробный аналог модели Фойхта. Пусть E0 и η - постоянные величины. Определяющее соотношение для модели Фойхта (29) записывается в виде α ηD0t ε + E0 ε = σ(t). (35) Видно, что при α → 1 соотношение (35) переходит в классическое реологическое соотношение модели Фойхта (3). Найдем решение дифференциального уравнения (35) с начальным условием ε(0+) = 0 в классе функций ε(t), таких что ε(t) ∈ C(0, T ) ∩ L(0, T ), T > 0. Учитывая ˙ α тождество (17), применим интегральный оператор I0t к (35): α ε(t) - λI0t ε = 1 α I σ, η 0t (36) где λ = -E0 /η. Уравнение (36) - интегральное уравнение Вольтерры второго рода, решение которого имеет вид ε(t) = 1 α I - λI0t η -1 α I0t σ = 1 1 α α,α α α,α α I + λE0t;λ I0t σ = I0t σ + λE0t;λ I0t σ = η η 1 α 1 α,α 2α,α = I0t σ + λE0t;λ σ = E0t;λ σ. (37) η η µ,α Здесь использованы некоторые свойства оператора E0x;λ и тождество [26] 2α,α α,α α λE0t;λ f = E0t;λ f - I0t f. Для σ(t), определённой в (33), получим равенство ε(t) = σ0 α,α σ0 α,α α,α α,α E0t;λ H(t) - H(t - t1 ) = E0t;λ H(t) - (E0t1 ;λ + Et1 t;λ )H(t - t1 ) . η η Используя тождества [27] µ;σ Eat;λ (t - a)α-1 = Γ(α)(t - a)µ+α-1 Eσ (λ(t - a)σ ; µ + α) = = Γ(α) Exp(σ, µ + α; λ; t - a), 177 О г о р о д н и к о в Е. Н., Р а д ч е н к о В. П., У н г а р о в а Л. Г. зависимости деформации от времени придадим вид ε(t) = σ0 H(t) Exp(α, α + 1; λ; t) - H(t - t1 ) Exp(α, α + 1; λ; t - t1 ) . η Преобразуем полученное выражение. Используя некоторые свойства функции Миттаг-Леффлера [25], найдём 1 Exp(α, α + 1; λ; t) = tα Eα (λtα ; α + 1) = (λtα )Eα (λtα ; α + 1) = λ η η α Eα (λt ; 1) - 1 = 1 - Exp(α,1; λ; t) . =- E0 E0 Тогда ε(t) = σ0 H(t) 1 - Exp(α,1; λ; t) - H(t - t1 ) 1 - Exp(α,1; λ; t - t1 ) E0 . (38) Нетрудно видеть, что lim ε(t) = α→1- = σ0 H(t) 1 - Exp(1, 1; λ; t) - H(t - t1 ) 1 - Exp(1, 1; λ; t - t1 ) = E0 σ0 = H(t)(1 - e-βt ) - H(t - t1 )(1 - e-β(t-t1 ) ) , E0 где β = -λ. Это следует из того, что lim Exp(α,1; λ; t) = E(1, 1; λ; t) = eλt . α→1- Таким образом, при α → 1 решение (38) принимает вид классического решения для модели Фойхта (3) . В работе [26] рассмотрены некоторые иные аналоги дробной модели Фойхта. Например, при α ∈ (0, 1/2) было рассмотрено определяющее соотношение следующего вида: 1-α α E0 D0t ε + ηD0t ε = σ(t). (39) По отношению к искомой функции ε(t) уравнение (39) является обыкновенным дифференциальным уравнением дробного порядка (1 - α), причем при α → 0+ уравнение (39) принимает вид классического реологического соотношения для модели Фойхта (3). Его решение для σ(t) = σ0 = const найдено в [26] в виде ε(t) = σ0 tα - Exp 1 - 2α, 1 + α; -E0 /η; t E0 Γ(1 + α) . В работе [26] также рассмотрен аналог дробной реологической модели Фойхта 1+α α E0 D0t ε + ηD0t ε = σ(t), 178 (40) Математическое моделирование наследственно упругого деформируемого тела . . . который обладает тем же свойством, что и определяющее соотношение (39), а именно при α → 0 уравнение (40) принимает вид классического соотношения (3). Решение (40) для σ = σ0 = const найдено в виде ε(t) = tα σ0 - Exp 1, 1 + α; -E0 /η; t E0 Γ(1 + α) . Приведенные дробные аналоги модели Фойхта, явные решения соответствующих дифференциальных уравнений, рассматривались в работе [26] в качестве примеров в стадии нагружения, а в работах [41-43] был проведен детальный анализ решений задачи ползучести при различных зависимостях напряжения от времени, доказана непрерывная зависимость найденных решений от параметра α при его устремлении к единице или к нулю, изучена асимптотика этих решений. 2.3. Дробный аналог модели Максвелла. Пусть E0 , η - постоянные величины. Рассмотрим определяющее соотношение (25) в виде, аналогичном (20): σ 1 α D0t ε - = σ(t). (41) E0 η Найдем решение этого дифференциального уравнения с начальным условием ε(0+) = σ(0+)/E0 . α Применяя к левой и правой частям уравнения (41) оператор I0t , получим равенство σ(t) 1 α ε(t) = + I0t σ(t). (42) E0 η Если σ(t) изменяется по закону (33), то решение (42) принимает вид ε(t) = σ0 α σ0 α α H(t) - H(t - t1 ) + I H(t) - (I0t1 + It1 t )H(t - t1 ) = E0 η 0t βtα β(t - t1 )α σ0 1+ H(t) - 1 + H(t - t1 ) , (43) = E0 Γ(α + 1) Γ(α + 1) где β = E0 /η. Очевидно, что существует предел σ0 lim ε(t) = (1 + βt)H(t) - 1 + β(t - t1 ) H(t - t1 ) , α→1-0 E0 который является классическим решением задачи о ползучести для модели Максвелла (4). Хорошо видно, что, как и в классическом случае, существует предел lim ε(t) = σ0 /E0 , t→0+ характеризующий величину мгновенной упругой деформации, совпадающую с мгновенной упругой разгрузкой. Действительно, значение деформации при t = t1 : σ0 βtα 1 ε(t1 ) = 1+ . E0 Γ(α + 1) 179 О г о р о д н и к о в Е. Н., Р а д ч е н к о В. П., У н г а р о в а Л. Г. С другой стороны, существует предел lim ε(t) = t→t1 +0 σ0 βtα 1 . E0 Γ(α + 1) Тогда ε(t1 ) - lim ε(t) = t→t1 +0 σ0 . E0 2.4. Дробные аналоги моделей Кельвина и Зенера. Рассмотрим определяющие соотношения (30) и (31) одновременно, записывая их следующим образом: α α σ(t) + bD0t σ = Eε(t) + aD0t ε. (44) Соотношению (44) можно придать такой вид: b E b E 1 b α D0t ε - σ + ε- σ = - σ(t). a a a a E a Относительно функции b u(t) = ε(t) - σ(t) a получим дифференциальное уравнение Барретта [20, 40]: α D0t u + βu = 1 (1 - bβ)σ(t), a (45) где β = E/a. Решение дифференциального уравнения (45) будем искать в классе функций u(t), таких, что u(t) ∈ C(0, T ) ∩ L(0, T ), T > 0, с начальным условием ˙ lim u(t) = u(0+) = 0. t→0+ (46) Для таких функций u(t) справедливо тождество (17): α α I0t D0t u = u(t), если α ∈ (0, 1). Тогда дифференциальное уравнение (45) с условием (46) эквивалентно редуцируется к интегральному уравнению типа Вольтерры второго рода 1 α α u(t) - λI0t u = (1 - bβ)I0t σ, (47) a где λ = -β = -E/a. Интегральное уравнение (47) имеет ту же структуру, что и уравнение (36), и, следовательно, его решение записывается в виде, аналогичном формуле (37): 1 α,α u(t) = (1 - bβ)E0t;λ σ. (48) a 180 Математическое моделирование наследственно упругого деформируемого тела . . . Для закона изменения напряжения (33) решение (48) записывается в виде σ0 (1 - bβ) H(t) Exp(α, α + 1; λ; t) - H(t - t1 ) Exp(α, α + 1; λ; t - t1 ) . a u(t) = Отметим, что u(t) = σ0 (1 - bβ) H(t) 1 - Exp(α, 1; -β; t) - E - H(t - t1 ) 1 - Exp(α, 1; -β; t - t1 ) , (49) так как Exp(α, α + 1; λ; t) = Exp α, α + 1; -E/a; t = a 1 - Exp(α,1; λ; t) . E Найдём законы деформирования для моделей Кельвина и Зенера в терминах параметров этих моделей: - модель Кельвина: a= ηE2 , E1 + E2 b= η , E1 + E2 β= E1 E = , a η b 1 = , a E2 E= E1 E2 , E1 + E2 1 - βb = E2 ; E1 + E2 из равенства (49) имеем ε(t) = σ(t) 1 α,α + (1 - bβ)E0t;λ σ = E2 a σ0 E2 = H(t) 1 + 1 - Exp α, 1; -E1 /η; t - E2 E1 E2 - H(t - t1 ) 1 + 1 - Exp α, 1; -E1 /η; t - t1 E1 ; (50) - модель Зенера: a= β= ε(t) = η(E1 + E2 ) , E1 E E0 = , a η E0 = b= η , E1 E1 E2 , E1 + E2 E = E2 , 1 - βb = E2 ; E1 + E2 σ0 E1 H(t) 1 + 1 - Exp α, 1; -E0 /η; t - E1 + E2 E2 E1 - H(t - t1 ) 1 + 1 - Exp α, 1; -E0 /η; t - t1 E2 . (51) 181 О г о р о д н и к о в Е. Н., Р а д ч е н к о В. П., У н г а р о в а Л. Г. Полученные в формулах (50) и (51) законы деформации в стадии нагружения и разгрузки совпадают с приведенными в работах [41-43]. А в настоящей работе изучена асимптотика найденных решений и показано равенство мгновенных упругих деформаций в стадии нагружения и разгрузки и их совпадение с классическими значениями. Прежде всего отметим непрерывную зависимость решений (50) и (51) от параметра α в том смысле, что при α → 1-0 эти решения переходят в хорошо известные решения для классических моделей Кельвина (5) и Зенера (6). Для дробного аналога модели Кельвина существует предел lim ε(t) = σ0 /E2 , α→0+ а в момент разгрузки значение деформации примет вид ε(t1 ) = σ0 E2 1+ 1 - Exp α, 1; -E1 /η; t1 E2 E1 . С другой стороны, lim ε(t) = t→t1 +0 σ0 1 - Exp α, 1; -E1 /η; t1 E1 . Тогда разность ε(t1 ) - lim ε(t) = σ0 /E2 . t→t1 +0 Тем самым показано, что величины мгновенной упругой деформации для дробного аналога модели Кельвина в стадии нагружения и разгрузки совпадают с аналогичным значением в классическом случае. Для дробного аналога модели Зенера имеет место аналогичный результат: lim ε(t) = ε(t1 ) - lim ε(t) = t→0+ t→t1 +0 σ0 . E1 + E2 3. Асимптотика найденных решений. Основные результаты этого раздела сформулированы в нижеследующих теоремах о явных решениях задачи ползучести при законе нагружения, определённом формулой (33), и их асимптотических свойствах. Используется следующий факт, касающийся асимптотического поведения дробной (двупараметрической) экспоненциальной функции [24]: Exp(α, µ; -β; t) = tµ-1 Eα (-βtα ; µ) = tµ-1-α + O(tµ-1-2α ), βΓ(µ - α) который следует из хорошо известной асимптотической формулы [25]: p Eα (z; µ) = - k=1 если α0 | arg z| z -k + O |z|-1-p , Γ(µ - kα) π, где πα < α0 < min {π, απ} . 2 182 Математическое моделирование наследственно упругого деформируемого тела . . . В нашем случае arg(-βtα ) = π. Во всех нижеследующих теоремах идёт речь о функции деформации ε(t), такой, что ε(t) ∈ C(0, t1 ] ∪ C(t1 , T ) ∩ L(0, T ). ˙ Теорема 1. Решение задачи ползучести для дробной математической модели Скотт-Блэра-Герасимова (32) с начальным условием lim ε(t) = 0 t→0+ существует в виде (34), единственно и обладает следующими свойствами: lim ε(t) = +∞, t→+∞ если σ(t) = σ0 для любого t ∈ (0, +∞), и lim ε(t) = 0, t→+∞ если σ(t) определяется формулой (33). Теорема 2. Решение задачи ползучести для дробного аналога модели Фойхта (35) с начальным условием lim ε(t) = 0 t→0+ существует в виде (38), единственно и обладает следующими свойствами: lim ε(t) = σ0 /E0 , t→+∞ если σ(t) = σ0 для любого t ∈ (0, +∞), и lim ε(t) = 0, t→+∞ если σ(t) определяется формулой (33). Замечание. Сравнительный анализ показывает, что в качественном отношении дробные модели Скотт-Блэра-Герасимова и Фойхта ведут себя подобным образом, обладая пределом lim ε(t) = 0 t→+∞ в стадии разгрузки. Однако модель Скотт-Блэра-Герасимова при постоянно действующей нагрузке демонстрирует неограниченный рост деформации, в то время как дробный аналог модели Фойхта обладает асимптотой при t → +∞, как и в классическом случае. Следующие три теоремы касаются решений задачи ползучести в математических моделях второго типа. Теорема 3. Решение задачи ползучести для дробного аналога модели Максвелла (41) с начальным условием lim ε(t) = (1/E0 ) lim σ(t) = σ0 /E0 t→0+ t→0+ 183 О г о р о д н и к о в Е. Н., Р а д ч е н к о В. П., У н г а р о в а Л. Г. существует в виде (43), единственно и обладает следующими свойствами: lim ε(t) = +∞, t→+∞ если σ(t) = σ0 для любого t ∈ (0, +∞), и lim ε(t) = 0, t→+∞ если σ(t) определяется формулой (33). Теорема 4. Решение задачи ползучести для дробного аналога модели Кельвина с начальным условием lim ε(t) = (1/E2 ) lim σ(t) = σ0 /E2 t→0+ t→0+ существует в виде (50), единственно и обладает следующими свойствами: lim ε(t) = σ0 /E, t→+∞ где E = E1 E2 /(E1 + E2 ), если σ(t) = σ0 для любого t ∈ (0, +∞), и lim ε(t) = 0, t→+∞ если σ(t) определяется формулой (33). Теорема 5. Решение задачи ползучести для дробного аналога модели Зенера с начальным условием lim ε(t) = 1/(E1 + E2 ) lim σ(t) = σ0 /(E1 + E2 ) t→0+ t→0+ существует в виде (51), единственно и обладает следующими свойствами: lim ε(t) = σ0 /E2 , t→+∞ если σ(t) = σ0 для любого t ∈ (0, +∞), и lim ε(t) = 0, t→+∞ если σ(t) определяется формулой (33). Доказательство существования решений во всех приведенных выше теоремах осуществляется непосредственной подстановкой решений в определяющие соотношения; единственность является следствием единственности решений интегральных уравнений вольтерровского типа. 4. Пример использования дробной модели Фойхта для описания экспериментальных данных. Одной из задач является решение проблемы применимости дробных моделей вязкоупругого поведения материалов для построения соответствующих феноменологических теорий на основе известных экспериментальных данных. В качестве примера рассмотрим дробный аналог модели Фойхта, задаваемой дифференциальным уравнением (35). Это типичная 184 Математическое моделирование наследственно упругого деформируемого тела . . . трёхпараметрическая модель относительно величин α, E0 и η. При кусочнопостоянном напряжении σ = σ0 = const (см. формулу (33)) решение (35) имеет вид (38), где λ = -E0 /η. Таким образом, для идентификации параметров α, E0 и η достаточно иметь кривые прямой и обратной деформации вязкоупругости при напряжении, изменяющемся по закону (33). Для построения соответствующей дробной феноменологической модели (35) воспользуемся экспериментальными данными для вязкоупругой деформации поливинилхлоридного пластиката, которые приведены в работе [44] и были получены при растяжении поливинилхлоридных трубок площадью поперечного сечения 1.2 мм2 при температуре 20 ℃ и постоянных напряжениях σ0 = 4.655; 6.288; 8.738; 10.372; 12.005 МПа, действующих на образцы в течение 8 часов, после чего производилась полная разгрузка образцов (σ0 = 0 при t > 8 часов). При каждом напряжении испытывалось от 3 до 5 образцов, а затем результаты испытаний осреднялись. Экспериментальные кривые вязкоупругой деформации поливинилхлоридного пластиката имеют мгновенно упругую компоненту и собственно вязкоупругую составляющую. Поэтому мгновенно упругая деформация в момент приложения нагрузки t = 0 + 0 при σ0 = const и в момент разгрузки t = 8 + 0 часов не учитывалась, т. е. из общей деформации вычиталась чисто упругая деформация, и в дальнейшем учитывалась только наследственная компонента деформации. На рис. 1 и 2 сплошными линиями приведены соответствующие осреднённые экспериментальные данные наследственной (вязкоупругой) деформации. Рассмотрим методику идентификации параметров дробной модели Фойхта (35). С точки зрения регрессионного анализа решение (38) представляет нелинейную зависимость относительно параметров α, E0 и λ, поэтому применение классического метода наименьших квадратов даёт нелинейную систему алгебраических уравнений, решение которой (с учётом того, что функция Exp(α, 1; λ, t) задаётся рядом) связано с непреодолимыми математическими трудностями. В связи с этим разработан итерационный численный метод нахождения этих параметров. Для получения нулевого приближения функция Exp(α, 1; λ, t) раскладывается в ряд, в котором используются лишь первые два члена и применяется метод наименьших квадратов к полученной линейной регрессионной модели. Дальнейшее уточнение параметров α, E0 и λ (а значит, и η) осуществляется методом координатного спуска. В результате первичной обработки каждой кривой вязкоупругого деформирования при всех пяти напряжениях σ0 и и времени от 0 до 12 часов найдеРис. 1. Экспериментальная (сплошная линия) и расчетная по модели (38) (штриховая линия) кривые вязкоупругого деформирования поливинлхлоридного пластиката при напряжении σ0 = 6.288 МПа с последующей разгрузкой [Figure 1. Experimental (solid line) and calculated by the model (38) (dashed line) viscoelastic creep curves of the flexible PVC at the stress σ0 = 6.288 MPa with subsequent unloading] 185 О г о р о д н и к о в Е. Н., Р а д ч е н к о В. П., У н г а р о в а Л. Г. a b c d Рис. 2. Экспериментальная (сплошная линия) и расчетная по модели (38) с осреднёнными значениями параметров (штриховая линия) кривые вязкоупругого деформирования поливинлхлоридного пластиката при напряжениях σ0 = 4.655 МПа (a), σ0 = 6.288 МПа (b), σ0 = 8.738 МПа (c), σ0 = 10.372 МПа (d) и σ0 = 12.005 МПа (e) с последующей разгрузкой e [Figure 2. Experimental (solid line) and calculated by the model (38) with averaged values of parameters (dashed line) viscoelastic creep curves of the flexible PVC at the stresses σ0 = 4.655 MPa (a), σ0 = 6.288 MPa (b), σ0 = 8.738 MPa (c), σ0 = 10.372 MPa (d), and σ0 = 12.005 MPa (e) with subsequent unloading] Значения параметров аппроксимации (38) и погрешности [The values of the approximation parameters for Eq. (38) and measure of inaccuracy] ¯ σ0 , MPa α E0 , MPa η ∆1 , % α ¯ E0 , MPa η ¯ 4.655 6.288 8.738 10.372 12.005 186 0.355 0.326 0.394 0.318 0.400 145.525 133.601 140.194 97.571 127.962 190.067 158.793 128.956 158.840 113.856 3.102 5.078 3.768 4.163 3.700 0.3846 137.155 140.191 ∆2 , % 18.580 8.892 4.841 4.841 12.273 Математическое моделирование наследственно упругого деформируемого тела . . . ны значения параметров модели (35), которые приведены в таблице. В колонке ∆1 приведены отклонения теоретической зависимости ε(tk ), вычисленной по формуле (38) при каждом значении σ0 с найденными параметрами α, E0 и η, от экспериментальных данных ε(tk ): ˜ n n 2 ε(tk ) - ε(tk ) / ˜ ∆= k=1 ε(tk ) ˜ 2 · 100 %. (52) k=1 Здесь tk - точки дискретизации времени, n - множество точек дискретизации. В качестве примера на рис. 1 штриховой линией приведена расчётная зависимость для вязкоупругой деформации ε(t), вычисленная по формуле (38) при σ0 = 6.288 МПа с параметрами из второй строки таблицы. Вычисленное по формуле (52) отклонение расчётных данных вязкоупругой деформации от экспериментальных значений в этом случае имеет величину ∆1 = 5.078 %. Аналогичный расчёт выполнен и для других уровней напряжения, полученные значения параметров α, E0 , η и величины ∆1 приведены в таблице. Из результатов, представленных в таблице, следует, что имеется существенный разброс относительно значений параметров α, E0 и η. Поэтому для построения модели (35) выполнено осреднение величин α, E0 и η по пяти реализациям (для всех пяти уровней напряжения σ0 ) и получены значения усреднённых параметров α, E0 , η , приведённые в таблице. ¯ ¯ ¯ На рис. 2 штриховыми линиями показаны расчётные значения вязкоупругой деформации для всех пяти уровней напряжения, построенных при осреднённых значениях параметров. В этом случае погрешность отклонения расчётных значений от экспериментальных для каждого уровня напряжений, вычисленная по формуле (52), приведена в таблице в колонке ∆2 . Средняя погрешность ∆2 по всем пяти реализациям составляет величину 11 %. Отметим, что для вычислений дробной экспоненциальной функции Exp(α, µ; λ; t) использовался «Автоматизированный исследовательский комплекс «MitLef».1 Анализ данных расчёта для дробной модели (35) и экспериментальных данных свидетельствует, что, по всей видимости, кривые стационарной вязкоупругости поливинилхлоридного пластиката обладают слабой нелинейностью связи «напряжение - деформация», а это не заложено в модель (35). Поэтому дальнейшее уточнение дробных аналогов модели Фойхта применительно к моделированию реологического деформирования поливинилхлоридного пластиката связано с введением нелинейности, но это не являлось целью настоящей работы. Резюме. В работе рассмотрены математические модели наследственно упругого тела: модель Скотт-Блэра-Герасимова, дробные аналоги моделей Фойхта, Максвелла, Кельвина и Зенера. Показан переход от структурных реологических моделей к математическим моделям наследственно упругого тела на основе гипотезы В. Вольтерры (7) с использованием в определяющем соотношении ядра Абелевского типа и применения оператора дробного дифференцирования Римана-Лиувилля с началом в нуле. Найдены явные решения задачи ползучести для перечисленных дробных дифференциальных уравнений и обоснована корректность постановки задачи Коши при использовании дробной производной Римана-Лиувилля. Отмечено, что многие ав187 О г о р о д н и к о в Е. Н., Р а д ч е н к о В. П., У н г а р о в а Л. Г. торы применяют иные виды дробных производных (Грюнвальда-Летникова, Лиувилля, Капуто), упрощающих постановку и решение начальных задач соответствующих интегро-дифференциальных уравнений, но искажающих суть моделируемых реальных физических процессов. Показано, что в предельном переходе при устремлении параметра дробности α к единице решения задач ползучести для дробных аналогов реологических моделей переходят в хорошо известные решения для классических реологических моделей. Сформулированы и обоснованы теоремы о существовании и асимптотических свойствах найденных решений задач ползучести на основе дробных моделей. В качестве иллюстрации приведён пример построения феноменологической модели поливинилхлоридного пластиката на основе дробной модели Фойхта (40). Разработана и реализована методика идентификации параметров этой модели. Приведены расчётные данные для деформации ползучести в режиме «нагрузка - разгрузка». Наблюдается удовлетворительное соответствие расчетных и экспериментальных данных. ORCIDs

About the authors

Eugeniy N Ogorodnikov

Samara State Technical University

Email: eugen.ogo@gmail.com
(Cand. Phys. & Math. Sci.; eugen.ogo@gmail.com; Corresponding Author), Dept. of Applied Mathematics & Computer Science 244, Molodogvardeyskaya st., Samara, 443100, Russian Federation

Vladimir P Radchenko

Samara State Technical University

Email: radch@samgtu.ru
(Dr. Phys. & Math. Sci.; radch@samgtu.ru), Head of Dept., Dept. of Applied Mathematics & Computer Science 244, Molodogvardeyskaya st., Samara, 443100, Russian Federation

Luiza G Ungarova

Samara State Technical University

Email: algluiza@gmail.com
Postgraduate Student, Dept. of Applied Mathematics & Computer Science 244, Molodogvardeyskaya st., Samara, 443100, Russian Federation

References

Supplementary files