Modeling of the study of viscoelastic deformation of elastic bodies

Full Text

1. Введение Механические свойства наследственно-упругих материалов, к которым можно отнести многие горные породы, пластмассы, металлы и сплавы, а также полимерные и композитные материалы на полимерной основе, при длительном нагружении в значительной степени зависят от времени. В условиях эксплуатации элементы конструкций из таких материалов подвергаются постоянным или меняющимся во времени режимам нагружения, в результате действия которых происходит перераспределение напряжений и деформаций, заметно изменяющее картину распределения напряжений и деформаций по сравнению с упругим случаем. Уточнение перераспределения напряжений и деформаций может быть полезным при оценке долговечности и несущей способности элементов конструкций. Одним из наиболее используемых при решении задач линейной теории вязкоупругости является метод, основанный на принципе В. Вольтерра, идея которого заключается в следующем: для решения задач вязкоупругости необходимо построить решение задачи линейной теории упругости и в окончательном результате заменить упругие постоянные операторами, расшифровав полученные комбинации операторов или функции от операторов. Непосредственное применение принципа Вольтерра к анализу напряженно деформированного состояния анизотропных сред является весьма затруднительным. Это связано с тем, что существует лишь незначительное число задач теории упругости, для которых известна явная зависимость от упругих постоянных, если использовать упругое решение в обычной форме. В случае, когда упругие решения содержат дробно-рациональные функции от упругих постоянных, решение задач вязкоупругости может быть легко получено при помощи алгебры резольвентных операторов Вольтерра, которая обладает следующими свойствами [1]: Свойство 1 (теорема умножения). Для любых регулярных точек x и y имеет место равенство P ∗ (x) P ∗ (y)= [P ∗ (x) - P ∗ (y)] / (x - y) , x ×= y. Свойство 2. Для любого резольвентного оператора P ∗ (β) имеет место формула обращения [1 - λP ∗ (β)]-1 =1 + λR∗ (β + λ) . Оператор R¯λ =1 + λR∗ (β + λ) называется обратным для оператора P¯λ =1 - λP ∗ (β) , а оператор R∗ (β + λ) - резольвентным для P ∗ (β) . Резольвентный оператор образован из функций того же класса, что и P ∗ (β) , сдвигом по параметру λ. В общем случае (например, для анизотропных сред) это решение может содержать иррациональные или трансцендентные функции от упругих постоянных. Это делает невозможным непосредственное использование алгебры резольвентных операторов к решению таких задач. Поэтому для решения упругих задач применяются приближенные методы, идея которых состоит в поиске решения, когда его зависимость от упругих постоянных устанавливается в явном виде. Степенные ряды, в которые раскладываются иррациональные и трансцендентные функции, как правило, сходятся медленно. Поэтому для получения достаточно точного решения нужно сохранять большое число членов ряда. Это приводит к громоздким вычислениям и потере точности. В задачах вязкоупругости для анизотропного тела, в отличие от задач теории упругости, вместо уравнений закона Гука используются уравнения состояния в одной из следующих форм: t ( e (t)= P (0) s (t)+ 0 d P (t - τ ) d (t - τ ) s (τ ) dτ, (1.1) t ( s (t)= R (0) e (t)+ 0 d R (t - τ ) d (t - τ ) e (τ ) dτ. (1.2) Здесь приняты обозначения e (t) = [e1, e2, e3, e4, e5, e6]T , s (t) = [s1, s2, s3, s4, s5, s6]T - векторы деформаций и напряжений в произвольный момент времени t; P (t) = Pij (t) , R (t) = Rij (t) (i, j = 1, 6) - матрицы функций ползучести и релаксации. В представлениях (1.1) и (1.2) учтены начальные условия. В начальный момент времени при t = 0 уравнения (1.1) и (1.2) переходят в уравнения закона Гука. Элементы матриц функций ползучести P или релаксации R определяются из эксперимента. Если представить их в виде произведений P (t) = {Pij (t)} = {aijpij (t)} , R (t) = {Rij (t)} = {Aijrij (t)} , то уравнения (1.1), (1.2) можно записать в виде ei = a¯ij sj, sj = A¯ij ei (i, j = 1, 6) (1.3) где ij a¯ij = aij 1+ p∗ , ij A¯ij = Aij 1 - r∗ . (1.4) Каждый из временных операторов (1.4) состоит из двух элементов: матриц упругих постоянij ных aij или Aij, а также интегральных операторов aijp∗ ij или Aijr∗ с ядрами наследственности dpij (t - τ ) или drij (t - τ ) , где d (t - τ ) d (t - τ ) t ( ∗ pijsj (t)= 0 t j ij j - dpij (t - τ ) s (τ ) dτ, r∗ e (t)= ( d (t - τ ) 0 j drij (t - τ ) e (τ ) dτ. d (t - τ ) Обращение систем (1.1), (1.2) или (1.3) представляет собой трудноразрешимую в аналитическом плане задачу. Экспериментально найденные значения ядер P или R задаются дискретным набором величин, соответствующим некоторым фиксированным значениям времени. По этим экспериментальным значениям различными методами строят аналитические аппроксимации ядер в специальной форме. Одной из таких форм является класс дробно-рациональных функций, построенных Ю. Н. Работновым [8] ∞ n n Э (-β; t - τ )= '\' (-1) β (t - τ )α+n(1+α) . (1.5) n=0 Γ [(n + 1) (1+ α)] Операторы с ядрами вида (1.5) обладают алгеброй операторов Вольтерра, и их резольвенты образованы из функций того же класса. Если принять dpij (t - τ ) d (t - τ ) ij α = λp Э p ij ij (-βp ; t - τ , то интегральный оператор aij представлений (1.4) примет вид a¯ij = aij 1+ λp Э∗ p p ( -β . ij αij ij Необходимо отметить, что аппроксимация экспериментально найденной кривой ползучести является источником дополнительных погрешностей. Так, ядро вида λЭα (-β; t - τ ) имеет произвол, содержащий всего три константы λ, α и β. Это обстоятельство не всегда позволяет построить функцию, хорошо описывающую экспериментально найденные величины на достаточно большом временном интервале. Результаты аппроксимации, приведенные в работах исследователей в этом направлении, показали, что погрешность в отдельных случаях достигала тридцати процентов. Поэтому точное описание вязкоупругих свойств реальных тел должно приводить к операторам с ядром более сложной формы. Для таких операторов построение резольвенты в аналитической форме наталкивается на непреодолимые трудности. В работах [1, 6] показано, что построение алгебры операторов Вольтерра не связано с каким-либо их специальным видом и может быть осуществлено для любых резольвентных операторов. Это обстоятельство делает возможным проводить алгебраическими методами реализацию решений граничных задач теории вязкоупругости во времени. Решение может быть осуществлено c использованием произвольных исходных операторов и выражено через значения воздействия на единицу этих операторов, заданных непосредственно таблицей экспериментальных данных. В данной работе предложен новый подход к решению таких задач. Предлагаемый подход не требует построения аналитического представления ядер ползучести и релаксации в специальной форме. Метод основан на численном определении резольвент интегральных уравнений состояния среды. 2. Обращение интегральных уравнений линейной теории вязкоупругости В опыте на ползучесть мгновенно прикладываются и поддерживаются постоянными напряжения s (t)= σ. Интегрирование уравнения (1.1) в этом случае дает e (t)= P (t) σ, (2.1) где σ = sk (t = 0)= [σ1, σ2, σ3, σ4, σ5, σ6]T - вектор-столбец напряжений констант. По формуле (2.1) определяются значения функции ползучести P (t) по величине измеренной деформации во времени. После интегрирования уравнения (1.2) по частям и подстановки в него выражения (2.1) можно найти соотношение, связывающее между собой функции ползучести P (t) и релаксации R (t) где I - единичная матрица. t ( I = R (t - τ ) 0 dP (τ ) dτ , (2.2) dτ Из уравнения (2.2) находится функция релаксации R (t) , соответствующая функции P (t) , найденной из эксперимента на ползучесть. В опытах на релаксацию мгновенно прикладываются и поддерживаются постоянными деформации. После интегрирования уравнения (1.2) получается s (t)= R (t) ε, (2.3) где ε = ek (t = 0)= [ε1, ε2, ε3, ε4, ε5, ε6]T - вектор-столбец деформаций констант. По формуле (2.3) определяются значения функции R (t) по величине измеренных напряжений во времени. После интегрирования уравнения (1.1) по частям и подстановки в него выражения (2.3) можно найти еще одно соотношение, связывающее между собой функции релаксации R (t) и ползучести P (t): t ( I = P (t - τ ) 0 dR (τ ) dτ . (2.4) dτ Из уравнения (2.4) находится функция ползучести P (t) , соответствующая функции R (t) , найденной из эксперимента на ползучесть. Решение уравнения (2.2) можно получить численным методом. В этом случае нет надобности использования аналитической аппроксимации экспериментальных значений функции P (t) . Пусть временной отрезок [0, t] разбит на N равных интервалов длиной h точками, соответствующими временам t0 = 0, t1 = h,... , ti = ih,... , tN = t, и пусть P0, P1,... , PN - экспериментально найденные значения функции ползучести P (t) в этих точках. Тогда соотношение (2.2) в точках сетки можно приближенно записать в виде n I = '\' Rn-i (Pi - Pi-1) (n = 0, 1,... , N ) , (2.5) i=0 где Pi = P (ti) , Rn-i = R (tn - ti) - значения функций в точке τ = ti. При этом значение Pi-1 = P (ti-1) при отрицательном значении индекса равно нулю. Из (2.5), с учетом изложенного, получается / n \ Rn = I - '\' Rn-i (Pi - Pi-1) i=1 A (n = 0, 1,... , N ) . (2.6) С помощью представлений (2.6) по экспериментально полученным значениям функции ползучести в точках сетки находятся функции релаксации. Решение уравнения (2.4) осуществляется аналогичным образом. В результате получается / n \ Pn = I - '\' Pn-i (Ri - Ri-1) i=1 a (n = 0, 1,... , N ) . (2.7) Из полученного соотношения по заданным функциям релаксации в точках сетки находятся функции ползучести. Таким образом, соотношения (2.6) и (2.7) полностью решают задачу обращения уравнений (1.1) или (1.2). Функции ползучести P (t) или релаксации R (t) строятся на основании табличных данных, которые определяются при проведении эксперимента. Таблица экспериментальных данных может содержать погрешности измерений, которые обусловлены различными причинами, в том числе связанными с длительностью интервалов замеров изменения деформаций или напряжений. Поэтому данные таблицы должны подлежать определенной математической обработке. Эта задача тесно связана с задачей сглаживания и восполнения, когда по заданным точкам измерения (tk, εk ) необходимо провести гладкую кривую w (t) при минимальной погрешности. При этом в качестве функции w (t) выбирают выражение, которое соответствует ожидаемому поведению рассматриваемого физического процесса. Установлено, что функции ползучести εk = P (tk ) должны быть монотонно возрастающими, а функции релаксации σk = R (tk) - монотонно убывающими функциями времени. Если уровень нагрузок меньше предела прочности материала, то развитие деформаций практически полностью прекращается при достижении определенного времени. В этом случае имеет место затухающий процесс, а представление для функции w (t) можно выбирать в виде ряда, содержащего затухающую экспоненциальную функцию / w(t)= m \ 1+ '\' aktαk e-βtα . (2.8) k=1 В представлении (2.8) подлежат определению величины α, β и ak, а процесс аппроксимации носит название нелинейное приближение. Пусть на отрезке [0, t] задана дискретная сетка значений временной переменной t: 0 = t0 < t1 < ... < tn = t. В узлах сетки даны экспериментально полученные значения функции ε [t]: ε(t0)= ε0, ε1, ... , εn = ε(t). Необходимо осуществить вычисление величин α, β и ak так, чтобы функция (2.8) сглаживала кривую, заданную таблицей, и, используя эту функцию, провести восполнение таблицы. Алгоритм сглаживания таблицы исходных экспериментальных данных основан на использовании дискретного метода наименьших квадратов. Для этого требуется, чтобы искомая сглаживающая функция минимизировала функционал n Φ (a1, a2, ... , am, α, β)= '\' [w (tp) - εp]2 . (2.9) p=0 Приравнивая производные по коэффициентам ak нулю, получаем систему линейных алгебраических уравнений относительно этих коэффициентов. Эта система решается при варьировании величин α и β для достижения выполнения условия (2.9). После нахождения постоянных α, β и ak определена функция, позволяющая осуществить гладкое восполнение сеточных значений в любой точке отрезка [0, t] . Матричные уравнения (2.1) и (2.3) можно записать в форме ei = Pij (t) sj или si = Rij (t) ej i, j = 1, 6 , (2.10) где Pij = aijpij (t) , Rij = Aijrij (t) . Для случая, когда из эксперимента на ползучесть определена матрица функций ползучести Pij (t) , необходимо взять представление ei = a¯ij sj = aij ( p∗ sj из (1.3) и сравнить его с соответствующим соотношением представления (2.10). 1+ ij Используя предлагаемую методику, получаем соотношение ∗ pij · 1= pij - 1. Таким образом, из экспериментальных данных для кривых ползучести, с учетом их восполнения, получены значения результатов воздействия на единицу интегральных операторов для записи уравнений состояния. Временная матрица в этом случае имеет вид ij Pn = a¯ij · 1= aij 1+ p∗ · 1 . Матрица A¯ij · 1= Rn находится по формуле (2.6). Для случая, когда из эксперимента на релаксацию определены элементы матрицы функций ре- 1 лаксации Rij (t) , необходимо взять представление si = A¯ij ej = Aij ( ij - r∗ ej из (1.3) и сравнить его с соответствующим соотношением представления (2.10). В результате получается ∗ rij · 1= rij - 1. Временная матрица уравнений закона Гука в момент времени t = tn определяются соотношениями ij Rn = A¯ij · 1= Aij 1 - r∗ · 1 . Матрица a¯ij · 1= Pn находится по формуле (2.7). Дальнейшее решение задач теории вязкоупругости в любой момент времени ничем не отличается от решения задач теории упругости. Возможности расшифровки произвольных операторов рассмотрены на примере для изотропного случая. Пусть из эксперимента получены опытные данные об изменениях коэффициента Пуассона во времени. Если не вводить функции Работнова, то кривые ползучести можно записать в виде ν ν¯ = ν (1 + p∗ · 1) = νpν (t) , (2.11) где pν (t) - функция ползучести, найденная в результате эксперимента. ν Отсюда находятся значения действия оператора p∗ на единицу: ∗ pν · 1= pν (t) - 1. (2.12) Используя предположение о постоянстве оператора объемного сжатия, согласно которому (1 - 2ν¯) /E¯ = (1 - 2ν) /E, находим ν E¯ · 1= E (1 - λE p∗ · 1) = E rE (t) , λE = 2ν/ (1 - 2ν) . (2.13) Здесь функция релаксации rE (t) для модуля Юнга определилась через функцию ползучести pν (t) для коэффициента Пуассона посредством формулы (2.12). При помощи свойства для резольвентных операторов [1 - λP ∗ (β)]-1 = 1 + λR∗ (β + λ) определяется оператор, обратный оператору (2.13): 1 1 1 ν E¯ · 1= E [1 + λE p∗ (λE ) · 1] = E pE (t) . (2.14) ν Результат воздействия оператора p∗ (λE ) на единицу представляет собой функцию ползучести pE (t) , значения которой выражаются через значения соответствующей функции релаксации rE (t) из соотношения (2.7): n pE (tn)=1 - '\' pE (tn-k) [rE (tk) - rE (tk-1)] (0, 1,... ,N ) . k=1 ν Значения p∗ (λE ) · 1 находится из соотношения (2.14): ∗ pν (λE ) · 1= [pE (t) - 1] /λE . (2.15) Учитывая свойство 1/G¯ = 2 (1 + ν¯) /E¯, можно найти 1 1 1 где ν G¯ · 1= G [1 + λG p∗ (λE ) · 1] = G pG (t) , λG = λE + ν/ (1 + ν) , ν G¯ · 1= G [1 - λG p∗ (βG) · 1] = G rG (t) , βG = λE - λG, (2.16) n rG (tn)= 1 - '\' rG (tn k=1 -k) [pG (tk ) - pG (tk-1 )] , (0, 1,... ,N ) , ∗ pν (βG) · 1= [1 - rG (t) ] /λG. ν Соотношения (2.11)-(2.16) связывают значения воздействия интегральных операторов p∗ (.. .) ν и r∗ (.. .) на единицу с функциями ползучести p...(t) и релаксации r... (t) , заданными таблично. Таким образом, решение граничных задач вязкоупругости сводится к решению задач теории упругости в произвольный момент времени. При этом уравнения закона Гука заменяются уравнениями состояния (1.1), (1.2), которые в результате интегрирования принимают вид (2.2) и (2.4), соответственно. 3. Численные исследования Проверка предложенного метода на достоверность осуществлялась исследованием напряженнодеформированного состояния пластины в условиях обобщенного плоского напряженного состояния. Решение упругой задачи для изотропной пластины, ослабленной круговым отверстием радиуса R, когда на бесконечности заданы усилия σ0 = p, σ0 = q, σ0 = 0, имеет вид [5] r r2 ( σ = 2α + b1 - β + 1 4a1 r2 - 2 3b3 r4 6 cos (2θ) , θ r2 ( σ = 2α - b1 + β - 3b3 r4 ( cos (2θ) , τrθ = β + 3b3 r4 - 2a1 r2 sin (2θ) , (3.1) r 1 г b1 ( a1 (κ +1) b3 l ur = 2G [α (κ - 1) - β cos (2θ)] + 2G - r + r - r3 cos (2θ) , r 1 ( a1 (κ - 1) b3 uθ = 2Gβ sin (2θ) - 2G + r r3 sin (2θ) . В представлениях (3.1) введены обозначения: α = (p + q) /4,β = (q -p)/2, κ =(3 - ν) / (1 + ν) , a1, b1, b3 - неизвестные постоянные; черточками сверху обозначены величины G и κ, которые при решении вязкоупругой задачи будут меняться во времени. В формулах (3.1) видна явная зависимость перемещений и напряжений от постоянных ν и G. Это обстоятельство позволяет построить точное решение вязкоупругой задачи методом Вольтерра и провести численные исследования. Рассмотрим два вида граничных условий. 1. Контур отверстия не подкреплен. Из граничных условий σr = τrθ = 0, когда r = R, коэффициенты функций (3.1) будут такими: a1 = -βR2, b1 = -2αR2, b3 = -βR4. (3.2) Из соотношений (3.1) и (3.2) видно, что в рассматриваемой задаче напряжения от времени не зависят, а изменяются только перемещения. 2. Контур отверстия жестко подкреплен. Коэффициенты функций (3.1) находятся из граничных условий ur = uθ = 0, когда r = R: G 2 G 2 G 4 a1 = GκβR , b1 = G αR (κ - 1) , b3 = GκβR . (3.3) Из соотношений (3.1) и (3.3) видно, что в рассматриваемой задаче от времени зависят и перемещения, и напряжения. В формулы для напряжений и перемещений (3.1) входят временные операторы. В первом случае это 1/G и κ/G, а во втором, с учетом значений коэффициентов (3.3), это операторы G, G/κ и 1/κ. При расшифровке их учитывалось, что соотношения между упругими операторами такие же, как и соотношения между упругими постоянными, а также считалось, что объемная деформация материала пластинки является упругой [8, 9]. Это позволяет воспользоваться равенствами 1 - 2ν - 1 2ν E = , G = - 3 ν , κ = . (3.4) E E 2 (1 + ν) 1+ ν В работе [9] коэффициент ν предложено заменить временным оператором α ν¯ = ν [1 + δ Э∗ (-β, t)] , (3.5) α где α, β и δ - реологические характеристики материала, а интегральный оператор Э∗ (-β, t) выражается через функцию Миттаг-Леффлера E1+α (-η) (см. [9]) Э∗ 1+α α (-β, t)= [1 - E1+α (-η)] /β, η = βt . (3.6) α Из равенств (3.4) с учетом (3.5) и свойств операторов Э∗ (-β, t) находятся указанные выше временные операторы. Исследования проводились для алюминиевой изотропной пластины со следующими упругими и реологическими параметрами [2]: E = 7.1 · 104 МПа, ν = 0.25, α = -0.5, δ = 0.00615c-0.5 и β = 0.00665c-0.5; а также для ортотропной пластины, изготовленной из композитного ма- 11 териала на основе эпоксидного связующего [7]: E0 = 23.0 · 103 МПа, λ1 = 0.0323 c-(1+α), 22 β1 = 0.157 c-(1+α), E0 12 = 16.0 · 103 МПа, λ2 = 0.1295 c-(1+α), β2 = 0.2745 c-(1+α), G0 = 3.08 · 103 МПа , λg = 0.0717 c-(1+α), βg = 0.0276 c-(1+α), ν12 = 0.11, ν21 = 0.0765, α = -0.846. Некоторые результаты для одного свободного или жестко подкрепленного отверстия в случае, когда пластинка растягивается на бесконечности усилиями интенсивности p вдоль оси Ox1, представлены ниже. Такие же результаты получились при решении этой задачи методом, изложенным в данной статье. При этом уравнения закона Гука следует выбирать в форме (2.1) для случая, когда контур Таб. 1 Перемещения / Movements t, час / t, hours 0 100 500 1000 5000 ur · G/(pR)(θ = 0) 1.2 3.8115 5.7501 6.6905 8.6486 ur · G/(pR)(θ = π/2) -0.4 -1.0529 -1.5375 -1.7726 -2.2622 uθ · G/(pR)(θ = π/4) 0.8 2.4322 3.6438 4.2316 5.4554 не подкреплен, и в форме (2.3) для подкрепленного контура. Матрицы Pn = P (tn) и Rn = R (tn) в случае, когда материал является изотропным, имеют вид ⎛ 1/E¯ -ν¯/E¯ 0 ⎞ ⎛ E¯/ 1 - ν¯2 ν¯E¯/ 1 - ν¯2 0 ⎞ Pn = ⎝ -ν¯/E¯ 1/E¯ 0 ⎠ , Rn = ⎝ ν¯E¯/ 1 - ν¯2 E¯/ 1 - ν¯2 0 ⎠ . (3.7) 0 0 1/G¯ 0 0 G¯ Элементы матрицы Rn в случае, когда материал изотропный, имеют достаточно простой вид. С целью тестирования элементы матрицы Rn вычислялись путем расшифровки и использования соотношений (2.6). Результаты вычислений сближаются, если шаг h сетки уменьшается. В случае, когда материал анизотропен, матрицу Rn следует вычислять, используя соотношения (2.6). При этом следует учитывать, что P0 = a, R0 = A. Исследование напряженно-деформированного состояния в любой момент времени ведется путем использования функций обобщенных комплексных переменных [4]. При проведении исследований для изотропных сред были использованы алгоритмы, разработанные для анизотропных материалов с учетом методики, предложенной в работе [3]. В таб. 1 приведены значения перемещений ur и uθ в характерных точках кругового отверстия изотропной пластинки в различные моменты времени. Как видно, стабилизация изменения перемещений происходит через достаточно большой промежуток времени. Напряжения в этом случае, как было указано ранее, от времени не зависят. Проведенные исследования для ортотропной пластинки показали, что напряжения, возникающие на контуре свободного отверстия, зависят от времени, а полученные численные результаты хорошо согласуются с данными, приведенными в работе [7]. Исследования напряженного состояния ортотропной пластинки с жестким круговым включением показали, что напряжения, возникающие на контуре спая, со временем уменьшаются за исключением небольшой зоны, находящейся вблизи возникновения максимальных напряжений σθ, где напряжения увеличиваются с течением времени. 4. Заключение Основные результаты представленной работы заключаются в следующем: 1. С применением аппарата алгебры резольвентных операторов и матричной алгебры разработан новый численно-аналитический метод обращения операторов определяющих уравнений деформирования анизотропных вязкоупругих сред без использования аналитических представлений ядер ползучести и релаксации, позволяющий распространить область возможного применения численно-аналитических подходов на новые классы моделей деформирования вязкоупругих тел и элементов конструкций. 2. Разработана численно-аналитическая методика исследования задач линейной вязкоупругости анизотропных сред, которая для произвольных моментов времени сводится к применению алгоритмов, подобных используемым при решении задач теории упругости анизотропного деформируемого тела без учета реологических эффектов, и позволяет изучить наследственные деформационные процессы в случаях конструкций, изготовленных из анизотропных материалов.

About the authors

R. N. Neskorodev

Donetsk State University

Author for correspondence.
Email: nromn72@mail.ru
Donetsk, Russia

A. V. Zyza

Donetsk State University

Email: z9125494@mail.ru
Donetsk, Russia

References

Supplementary files