The use of linear fractional analogues rheological models in the problem of approximating the experimental data on the stretch polyvinylchloride elastron

Full Text

Введение. В работе [1] на основе гипотезы В. Вольтерры о наследственно упругом деформируемом твердом теле и метода структурного моделирования рассмотрены дробные аналоги некоторых классических реологических моделей: модель Скотт Блэра как аналог и обобщение модели Ньютона для вязкой жидкости, а также дробные аналоги моделей Фойхта, Максвелла, Кельвина и Зенера. В указанной работе было отмечено, что если в исходном определяющем соотношении В. Вольтерры использовано ядро Абеля, то возникающие в определяющих соотношениях дробные производные будут являться производными Римана-Лиувилля на отрезке временной оси. Показана корректность классической задачи Коши для возникающих из определяющих соотношений дробных дифференциальных уравнений относительно некоторых линейных комбинаций функций напряжений и деформаций. Найдены явные решения задачи о ползучести при постоянном напряжении на стадиях нагружения и разгрузки и отмечено, что они могут быть использованы для аппроксимации экспериментальных зависимостей деформации от времени, полученных в опытах по одноосному растяжению деформируемых тел (физических сред) под действием постоянных нагрузок. Успешное решение задачи аппроксимации позволяет идентифицировать параметры реологической модели и тем самым определить вид и структуру математической модели вязкоупругого поведения изучаемой среды в форме определяющего соотношения между функциями напряжения и деформации. В работе [1] авторы ограничились примером использования дробного аналога модели Фойхта для описания данных эксперимента по нагружению образцов из поливинилхлоридного пластиката, которые приведены в работе [2]. В настоящей работе решается задача идентификации параметров математической модели вязкоупругого деформирования указанных выше математических объектов и выполняется оценка погрешности построенных моделей. Прежде чем перейти к описанию алгоритма построения математической модели одноосного вязкоупругого деформирования, отметим, что идея замены классических вязкоупругих моделей деформируемого твердого тела наследственно упругими моделями имеет давнюю историю, восходящую к работам В. Вольтерры [3, 4] и Х. Больцмана [5]. Она получила мощное развитие в работе Ю. Н. Работнова [6], которая, по существу, заложила начало школы по механике наследственно упругих твердых тел. Отметим здесь также несколько предшествующих работ [7-12]. К настоящему моменту опубликовано огромное количество работ, посвященных дробным реологическим моделям (см., например, библиографические списки в работах [13-17]). Как отмечено в монографии [15, с. 443], большинство авторов используют в определяющих соотношениях дробные производные Капуто, производные Грюнвальда-Летникова [18] или лиувиллевскую форму [19] дробного интегро-дифференцирования на всей числовой оси, в то время как дробные реологические модели выражаются через дробные производные Римана-Лиувилля. Для нахождения явных решений дифференциальных уравнений дробного порядка, возникающих из определяющих 692 Применение линейных дробных аналогов реологических моделей . . . соотношений, используется в основном преобразование Лапласа. В работе [1] показано, что они могут быть найдены собственно методами дробного анализа с использованием интегральных операторов с началом в нуле или в другой точке, соответствующей, например, конечному моменту времени. Проблема идентификации параметров дробных реологических моделей и необходимость их определения из экспериментальных данных обсуждалась уже авторами ранних работ по наследственно упругим моделям (см., например, [8]). В монографии Ю. Н. Работнова [20, с. 86] приведен метод прямого определения параметров дробно-экспоненциального ядра ползучести, предложенный ранее авторами работы [21]. Однако указанный метод основан на сравнении образов преобразований Лапласа явного (аналитического) решения задачи с образом полиномиальной аппроксимации экспериментальных точек и минимизации суммы квадратов отклонений методом координатного спуска в пространстве параметров модели. В действительности проблема может быть решена непосредственно путем сопоставления экспериментальных и теоретических значений. В многочисленных публикациях последних десяти лет проблема идентификации параметров дробных моделей в основном решается на теоретическом уровне, например, методами спектрального анализа [22]. В работе [23] параметры модели определяются исходя из нескольких характерных точек, полученных в эксперименте, путем подстановки значений деформации в аналитические решения соответствующей задачи. 1. Описание эксперимента и идентификация параметров математических моделей. Рассмотрим дробные аналоги известных структурных моделей вязкоупругого тела: модели Ньютона, Максвелла, Фойхта, Кельвина и Зенера с дробным порядком α ∈ (0, 1). В работах [1, 24] были найдены явные решения задачи о ползучести при постоянном напряжении в случае нагружения и разгрузки для всех вышеупомянутых моделей. Рассмотрим применение этих математических моделей к конкретному эксперименту. В работе [2] исследовано поведение поливинилхлоридных трубок площадью поперечного сечения 1.2 мм2 при температуре 20 ℃ и постоянных напряжениях σ0 ∈ 4.655; 6.288; 8.738; 10.372; 12.005 МПа, действующих на образцы в течение 8 часов, после чего производилась полная разгрузка образцов (σ = 0 при t > 8 часов). При каждом напряжении испытывалось от 3 до 5 образцов, а затем результаты испытаний осреднялись. Анализ экспериментальных кривых показал, что данный материал обладает следующими свойствами: имеет мгновенную упругую деформацию при приложении и снятии нагрузки, а под действием постоянной нагрузки проявляет свойство вязкоупругости. Простейшей математической моделью, описывающей подобное поведение деформируемых тел, является дробный аналог модели Максвелла α ηD0t ε(t) = η α D σ(t) + σ(t), E0 0t (1) α ϕ = (d/dt)n I n-α ϕ - левосторонняя дробная производная Римана- где D0t 0t Лиувилля порядка α > 0 функции ϕ(t) (t > 0, n = [α] + 1, [α] - целая часть β β числа α ∈ R), I0t ϕ = I0+ ϕ (t) - интеграл Римана-Лиувилля порядка β > 0 [15], E0 и η - постоянные величины, которые при α = 1 совпадают с модулем упругого элемента и коэффициентом демпфирования соответственно в клас- 693 У н г а р о в а Л. Г. сической модели Максвелла [1]. Для напряжения, изменяющегося по закону σ(t) = σ0 [H(t) - H(t - t1 )], (2) где σ0 = const > 0, H(t) - функция Хевисайда, t1 - момент времени разгрузки, решение задачи ползучести для модели (1) имеет вид ε(t) = σ0 E0 1+ βtα β(t - t1 )α H(t) - 1 + H(t - t1 ) , Γ(α + 1) Γ(α + 1) (3) где β = E0 /η, Γ(z) - гамма-функция [25]. Рассмотрим методику идентификации параметров решения для дробного аналога модели Максвелла (3). Будем аппроксимировать данные эксперимента теоретическими значениями деформации на стадии нагружения по формуле σ0 tα σ0 + . (4) ε(t) = E0 ηΓ(α + 1) Идентификации подлежат порядок дробности α и величины E0 , η. Мгновенная упругая деформация выражается соотношением ε(0+) = σ0 /E0 . Тогда E0 = σ0 /˜ ε(0), где ε˜(0) = ε˜0 - первое значение деформации, взятое из массива экспериментальных данных. Для нахождения остальных неизвестных величин воспользуемся методом наименьших квадратов. Пусть ti - набор n моментов времени t, при которых было произведено измерение экспериментальных значений деформации ε˜i = ε˜(ti ) в каждый момент времени ti . Тогда равенство (4) можно записать в виде ε˜i - ε˜0 = σ0 tαi . ηΓ(α + 1) Разделив левую и правую части этого равенства на величину σ0 и прологарифмировав, получим следующее выражение: ln(˜ εi - ε˜0 ) - ln σ0 = α ln ti - (ln η + ln Γ(α + 1)). (5) Введем следующие обозначения: ui = ln(˜ εi - ε˜0 ) - ln σ0 , τi = ln ti , c = ln η + ln Γ(α + 1). Тогда (5) запишется в виде ui = ατi - c. Найдем параметры α и c, минимизирующие функцию n (ui - ατi + c)2 → min . f= i=1 694 (6) Применение линейных дробных аналогов реологических моделей . . . Используя необходимые условия экстремума ∂f = 0, ∂α ∂f =0 ∂c и решив систему уравнений метода наименьших квадратов, получим формулы для определения значений параметров α и c: α= c= n n n n i=1 ui i=1 τi i=1 ui τi - , n n 2 n i=1 τi - ( i=1 τi )2 n i=1 ui τi n n i=1 τi - n 2 i=1 τi - ( n n 2 i=1 ui i=1 τi . 2 n i=1 τi ) Возвращаясь к введенным выше обозначениям, находим коэффициент η: η = ec /Γ(α + 1). Найденные значения параметров позволяют по формуле (3) найти теоретические значения деформации в указанные выше моменты времени. Отметим, что если в решении задачи о ползучести в рамках дробного аналога Максвелла (3) отбросить мгновенную упругую деформацию, то получим решение соответствующей задачи для модели Скотт Блэра: ε(t) = σ0 [tα H(t) - (t - t1 )α H(t - t1 )]. ηΓ(α + 1) (7) При аппроксимации экспериментальных данных функцией (7) вычтем из тестовых экспериментальных значений деформации величину мгновенной упругой деформации, совпадающую со значением ε˜(0). Тогда дальнейший поиск параметров α и η аналогичен алгоритму для модели Максвелла, при этом в минимизирующей функции (6) значения (8) ui = ln ε˜i - ln σ0 , так как ε˜0 = 0. Таким образом, параметры α и η в модели Скотт Блэра окажутся такими же, как и в дробном аналоге Максвелла. Рассмотрим теперь дробный аналог модели Кельвина. В работе [1] для этой модели приведено решение в случае, когда σ(t) определяется формулой (2): ε(t) = σ0 E2 H(t) 1 + E2 E1 1 - Exp α, 1; - -H(t - t1 ) 1 + E2 E1 E1 ;t η - 1 - Exp α, 1; - E1 ; t - t1 η , (9) где E1 , E2 - коэффициенты упругих элементов, Exp(α, µ; λ; t) = tµ-1 Eα (λtα ; µ) 695 У н г а р о в а Л. Г. - обобщенная (двупараметрическая) дробная экспоненциальная функция [26], Eα (z; µ) - функция типа Миттаг-Леффлера [27]. Из (9) следует, что идентификации подлежат коэффициенты α, E1 , E2 , η. Вновь ограничимся случаем нагружения σ(t) = σ0 = const (0 < t t1 ). Деформация образца описывается выражением ε(t) = σ0 σ0 + E2 E1 1 - Exp α, 1; - E1 ;t η . (10) Для нахождения параметра, характеризующего мгновенную упругую деформацию, используем первую экспериментальную точку ε˜(0). Тогда E2 = σ0 /˜ ε(0). Известно, что дробный аналог модели Кельвина имеет асимптоту [1] lim ε(t) = σ0 /E, t→+∞ где E = E1 E2 /(E1 + E2 ), если σ(t) = σ0 = const, t ∈ (0, +∞). Подставив в это выражение последнее экспериментальное значение деформации ε˜n при нагружении, определим параметр E1 ≈ σ0 /˜ εn в грубом приближении. Следующие параметры α и η также в грубом приближении будут определяться с помощью метода наименьших квадратов. Для этого разложим дробную экспоненциальную функцию Exp(α, µ; λ; t) в ряд [26] ∞ Exp(α, 1; -β; t) = k=0 (-β)k tαk -βtα β 2 t2α = 1+ + + . . . , (11) Γ(αk + 1) Γ(α + 1) Γ(2α + 1) где β = -E1 /η, и ограничимся первыми двумя членами ряда (11). Тогда зависимость деформации ε(t) в случае нагрузки при 0 < t t1 примет вид ε(t) - ε(0) ≈ σ0 βtα σ0 tα · = . E1 Γ(α + 1) ηΓ(α + 1) (12) Очевидно, что равенство (12) совпадает с решением задачи ползучести для дробного аналога модели Максвелла (4), для которого уже были найдены необходимые значения параметров. Дальнейшие уточнения параметров α, η и E1 осуществляются методом покоординатного спуска (методом Хука-Дживса) [28]. Таким образом, аппроксимация дробным аналогом модели Максвелла является грубым приближением для аппроксимации дробной моделью Кельвина. Найденные параметры подставим в решение (9) и вычислим среднеквадратические оценки отклонений теоретических расчетных значений от экспериментальных значений деформаций на стадиях нагрузки и разгрузки по формуле n i=1 696 n [˜ ε(ti ) - ε(ti )]2 / ∆= [˜ ε(ti )]2 · 100%. i=1 (13) Применение линейных дробных аналогов реологических моделей . . . Заметим, что исключив мгновенную упругую деформацию из уравнения (9), получим решение задачи о ползучести в рамках дробного аналога модели Фойхта: ε(t) = σ0 E1 H(t) 1 - Exp α, 1; - ; t E1 η - -H(t - t1 ) 1 - Exp α, 1; - E1 ; t - t1 η . (14) Вычитая из экспериментальных данных первое значение деформации ε˜(0), соответствующее упругой деформации, можно аппроксимировать полученные экспериментальные тестовые кривые уравнением (14). Для этого следует воспользоваться идеей, изложенной выше при использовании в качестве аппроксимирующей модели дробного аналога модели Кельвина (9). Минимизирующая функция в этом случае совпадает с минимизирующей функцией для модели Скотт Блэра (7), в которой ui задается формулой (8). Рассмотрим далее дробный аналог модели Зенера, для которого решение задачи о ползучести имеет вид [1] ε(t) = σ0 E1 + E2 H(t) 1 + E1 E2 1 - Exp α, 1; - -H(t - t1 ) 1 + E1 E2 E0 ;t η 1 - Exp α, 1; - - E0 ; t - t1 η , (15) где E0 = E1 E2 /(E1 + E2 ). Ограничимся стадией нагружения: σ(t) = σ0 = const (0 < t t1 ). Деформация будет определяться равенством ε(t) = σ0 E1 σ0 + E1 + E2 E2 (E1 + E2 ) 1 - Exp α, 1; - E0 ;t η . (16) Выделим мгновенную упругую компоненту деформации, благодаря которой вычисляется величина E1 + E2 = σ0 /˜ ε(0), где ε(0) - экспериментальное значение упругой деформации. Устремляя в уравнении (16) t → +∞, получим lim ε(t) = σ0 /E2 ≈ ε˜n , t→+∞ откуда следует, что E2 ≈ σ0 /˜ εn . Тогда E1 ≈ σ0 /˜ ε(0) - E2 ≈ σ0 (1/˜ ε(0) - 1/˜ εn ). 697 У н г а р о в а Л. Г. Решение (16) по структуре аналогично (10). Вновь разложим дробную экспоненциальную функцию в ряд и ограничимся первыми двумя членами. Тогда на стадии нагрузки из (16) получим ε(t) - ε(0) = σ0 E12 tα . (E1 + E2 )2 ηΓ(α + 1) Минимизируя функцию (6), в которой ui = ln (ε(t) - ε(0))(E1 + E2 )2 , σ0 E12 находим коэффициенты α и η в грубом приближении. Далее, следуя методу Хука-Дживса, аппроксимируем экспериментальные данные уравнением (15) при любом t и находим отклонение теоретической зависимости ε(ti ) с найденными параметрами α, E1 , E2 , η от экспериментальных данных ε˜(ti ) по формуле (13). Нетрудно увидеть, что при отбрасывании начального упругого значения деформации решение задачи о ползучести для дробного аналога модели Зенера принимает вид, подобный решению для дробного аналога модели Фойхта (14): E0 E1 σ0 1 - Exp α, 1; - ; t . ε(t) = E2 (E1 + E2 ) η Последнее равенство, в свою очередь, при замене дробной экспоненты на первые два члена ее ряда принимает вид решения уравнения Скотт Блэра (7), в котором значения параметров уже известны. 2. Анализ и сравнение результатов аппроксимаций. В результате первичной обработки каждой кривой вязкоупругого деформирования при всех пяти напряжениях σ0 при нагрузке σ(t) = σ0 = const (0 t 8), σ(t) = 0 (8 < t 12) найдены значения параметров для всех рассмотренных моделей. В табл. 1 приведены значения параметров аппроксимации для дробного аналога модели Кельвина. В колонке ∆ с использованием нормы (13) приведены отклонения теоретической зависимости ε(ti ), рассчитанной по формуле (9) при каждом значении напряжения σ0 с найденными параметрами Таблица 1 Значения параметров аппроксимации (9) и погрешности для дробного аналога модели Кельвина [The values of the approximation parameters for Eq. (9) and measure of inaccuracy for fractional Kelvin model] 698 σ0 , MPa α E1 , MPa E2 , MPa η ∆, % 4.655 6.288 8.738 10.372 12.005 0.355 0.326 0.394 0.318 0.400 145.525 133.601 140.194 97.571 127.962 1163.800 898.329 1028.000 829.733 857.500 190.067 158.793 128.956 158.840 113.856 2.586 4.097 3.168 3.468 3.121 Применение линейных дробных аналогов реологических моделей . . . α, E1 , E2 , η, от экспериментальных данных. В качестве примера на рис. 1 штриховой линией представлена полученная по формуле (9) расчетная зависимость при σ0 = 4.655 MПa с параметрами из первой строки табл. 1. Погрешность аппроксимации конкретной реализации в данном случае равна 2.586 %. Анализ данных табл. 1 свидетельствует о том, что, вообще говоря, для дробного аналога модели Кельвина наблюдается существенный разброс параметров модели для реализаций при различных значениях напряжений. 0.03 Viscoelastic creep ε(t) Рис. 1. Экспериментальная (сплошная линия) и расчетная по модели (9) (штриховая линия) кривые вязкоупругого деформирования поливинилхлоридного пластиката при напряжении σ0 = 4.655 МПа с последующей разгрузкой [Figure 1. Experimental (solid line) and calculated by the Kelvin’s fractional model (9) (dashed line) viscoelastic creep curves of the flexible PVC at the stress σ0 = 4.655 MPa with subsequent unloading] 0.02 0.01 0 0 4 8 12 Time (hours) Аналогичный расчет проведен и для других рассмотренных в настоящей работе моделей при всех пяти уровнях напряжений. Для построения исследуемых моделей, «работоспособных» при напряжениях σ0 от 4.655 МПа до 12.005 МПа, выполнено усреднение параметров по пяти реализациям для каждой из них. В табл. 2 представлены полученные значения усредненных ¯1 , E ¯2 , η¯ для всех рассмотренных моделей. параметров α ¯, E В табл. 3 приведены погрешности отклонений расчетных значений вязко- Усредненные значения параметров для дробных моделей [The averaged parameter values for the fractional models] α ¯ ¯1 , MPa E ¯2 , MPa E η¯ 0.179 0.359 0.179 0.359 0.285 - 128.971 - 128.971 1078.692 - - 955.472 955.472 92.418 289.847 150.102 289.847 150.102 143.029 Fractional Models The The The The The Scott Blair’s fractional model (7) Voigt’s fractional model (14) Maxvell’s fractional model (3) Kelvin’s fractional model (9) Zener’s fractional model (15) Таблица 2 Таблица 3 Погрешности аппроксимаций (∆, %) для всех исследуемых моделей после осреднения параметров для пяти реализаций [The approximation errors ∆ (in percentages) for all investigated models after averaging parameters for the five realizations] Fractional Models The The The The The model model model model model (7) (14) (3) (9) (15) Stresses σ0 , MPa 4.655 6.288 8.738 10.372 12.005 19.139 18.580 19.792 19.325 19.018 12.538 8.892 9.671 6.606 7.394 10.192 4.841 8.578 3.964 4.937 13.237 11.062 12.609 11.153 11.438 15.745 12.273 14.429 11.843 12.527 699 У н г а р о в а Л. Г. Viscoelastic creep ε(t) Viscoelastic creep ε(t) 0.04 0.02 0 0 4 8 0.04 0.02 0 12 Time (hours) 0 4 a Viscoelastic creep ε(t) Viscoelastic creep ε(t) 0.02 0 4 8 12 Time (hours) c Viscoelastic creep ε(t) 8 12 0.08 0.04 4 8 Time (hours) 0 0 4 d 0.04 0 0.04 Time (hours) 0.08 0 12 b 0.06 0 8 Time (hours) 12 Рис. 2. Экспериментальные (сплошные линии) и расчетные по модели (9) с усредненными значениями параметров (штриховые линии) кривые вязкоупругого деформирования поливинилхлоридного пластиката при напряжениях σ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 Kelvin’s fractional model (9) 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] 700 Применение линейных дробных аналогов реологических моделей . . . Viscoelastic creep ε(t) Viscoelastic creep ε(t) 0.05 0.03 0.015 0 0 4 8 0.025 0 12 0 4 8 12 0.08 Viscoelastic creep ε(t) Viscoelastic creep ε(t) 12 b a 0.06 0.03 0 8 Time (hours) Time (hours) 0 4 8 12 0.04 0 0 4 Time (hours) Time (hours) d c Viscoelastic creep ε(t) 0.1 0.05 0 0 4 8 Time (hours) 12 Рис. 3. Экспериментальные (сплошные линии) и расчетные по модели (3) с усредненными значениями параметров (штриховые линии) кривые вязкоупругого деформирования поливинилхлоридного пластиката при напряжениях σ0 = 4.655 МПа (a), σ0 = 6.288 МПа (b), σ0 = 8.738 МПа (c), σ0 = 10.372 МПа (d) и σ0 = 12.005 МПа (e) с последующей разгрузкой e [Figure 3. Experimental (solid line) and calculated by the Maxvell’s fractional model (3) 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] 701 У н г а р о в а Л. Г. упругой деформации от экспериментальных для всех пяти моделей с усредненными параметрами, вычисленные по формуле (13) при каждом уровне напряжений. Средняя же погрешность по всем пяти реализациям для модели Скотт Блэра равна 14.170 %, для дробных аналогов модели Фойхта - 11.130 %, Максвелла - 13.016 %, Кельвина - 10.578 % и Зенера - 11.063 %. Таким образом, наилучшее приближение теоретической кривой к экспериментальным значениям дает дробный аналог модели Кельвина. В качестве иллюстрации на рис. 2 и 3 штриховой линией построены расчетные значения вязкоупругой деформации по моделям (9) и (3) соответственно с усредненными коэффициентами. Отметим, что в настоящей работе все вычисления и графические построения производились в пакете MATLAB, а для вычисления дробной экспоненциальной функции Exp(α, µ; λ; t) использован «Автоматизированный исследовательский комплекс «MitLef».1 3. Выводы. Сравнительный анализ погрешностей отклонений расчетных кривых от экспериментальных для всех моделей показал, что дробный аналог модели Кельвина наилучшим образом аппроксимирует экспериментальные данные, полученные при растяжении поливинилхлоридного стержня. Однако по результатам вычислений средней погрешности отклонений расчетных кривых с усредненными параметрами для пяти уровней напряжений по каж¯1 , дой модели видно, что дробная модель Кельвина с коэффициентами α ¯, E ¯ E2 , η¯ не столь существенно улучшает аппроксимацию экспериментальных данных в сравнении с другими моделями. Близкие результаты дают также дробный аналог модели Фойхта с тремя параметрами и Зенера с четырьмя параметрами. Декларация о финансовых и других взаимоотношениях. Исследование не имело спонсорской поддержки. Автор несет полную ответственность за предоставление окончательной версии рукописи в печать. Окончательная версия рукописи была одобрена автором. Автор не получал гонорар за статью. ORCID

About the authors

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