CC BY
56
УДК 533.12 Дата подачи статьи: 31.03.2014
МОДЕЛИРОВАНИЕ ПРОЦЕССА ДЕФОРМИРОВАНИЯ ПРИ НАЛИЧИИ ФАЗОВОГО ПЕРЕХОДА ОТ ПОТЕНЦИАЛЬНОГО К ДИССИПАТИВНОМУ ТЕЧЕНИЮ
(Работа выполнена при частичной поддержке Региональной целевой программы развития вычислительных, телекоммуникационных и информационных ресурсов УрО РАН, проект № РЦП-13-П4)
В.П. Федотов, д.т.н., профессор, главный научный сотрудник; А.В. Горшков, к.ф.-м.н., доцент, научный сотрудник (Институт машиноведения Уральского отделения РАН, ул. Комсомольская, 34, г. Екатеринбург, 620049, Россия, fedotov@mach.uran. ru, A.V. Gorshkov@urfu. ru)
Перекрестные эффекты в сложных системах определяются феноменологическими составляющими соответствующих систем уравнений, в которых на основе подхода Онзагера устанавливается связь между обобщенными потоками и обобщенными силами. Для деформационных задач эта связь, как правило, устанавливается линейным постулатом Гука. Однако на его основе невозможно описать фазовый переход второго рода от упругости к диссипа-тивному пластическому течению, а непрерывные физические величины, такие как температура, концентрация примеси и т.п., связать с условным пределом текучести, поскольку последняя величина является экспериментально-точечной.
Как гипотезу, имеющую подтверждение на практике, примем, что коэффициенты взаимности Lkm, постоянные по Онзагеру, могут быть нелинейными функциями. Для деформационных задач примем постулат Коши, согласно которому в упругой области коэффициент Онзагера есть квадратичная функция деформации. В этом случае решение задачи упругости сводится к решению нелинейных уравнений. Применение метода конечных элементов требует решения систем большой размерности и, как следствие, больших затрат времени. Для сокращения времени счета используется метод граничных элементов в сочетании с применением там, где возможно, аналитических вычислений. В статье рассмотрен пример применения алгоритма к решению плоской задачи нелинейной упругости для неоднородной области.
Ключевые слова: нелинейная упругость, граничные элементы, аналитические вычисления.
MODELING OF THE DEFORMATION PROCESS IN THE PRESENCE OF A PHASE TRANSITION FROM A POTENTIAL TO DISSIPATIVE STREAM
(The work is done with partial support from the Regional special-purpose program of computational, telecommunication and information resources development of the Ural Division of the Russian Academy of Sciences, project № РЦП-13-П4) Fedotov V.P., Dr.Sc. (Engineering), Professor, Chief Researcher;
Gorshkov A. V., Ph.D. (Physics and Mathematics), Associate Professor, Research Associate (Institute ofEngineering Science of the Ural Branch of the Russian Academy of Sciences, Komsomolskaya St. 34, Ekaterinburg, 620049, Russian Federation, fedotov@imach.uran.ru, A.V.Gorshkov@urfu.ru)
Received 31.03.2014
Abstract. The paper considers the algorithm for solving nonlinear elasticity issue.
Cross-effect in complex systems is determined by the phenomenological components of relevant systems of equations. In these systems the relation between the generalized flows and generalized forces is established based on the Onsager's approach. For the deformation problems, this relation usually is established by Hook's linear postulate. However, this postulate can not be a basis for describing the second-kind phase transition from elasticity to dissipative plastic flow. And relating continuous physical quantities (temperature, admixture concentration etc.) can not offset yield strength since the latter is an experimental pointwise quantity.
As a practically confirmed hypothesis, it is assumed that the reciprocity factors Lkm, constant by Onsager, can be nonlinear functions. For deformation problems we use the Cauchy postulate. It states that the Onsager factor in the elastic region is a quadratic function of strain. In this case the elasticity problem solution is reduced to the solution of nonlinear equation. The application of the finite elements method requires a solution of large-dimension systems and, consequently, it is extremely time-consuming. To reduce the computation time, the boundary element method is used in combination with analytical computation if possible. The article considers the example of the algorithm application to solve a nonlinear plane problem of elasticity for an inhomogeneous region.
Keywords: nonlinear elasticity, boundary element, analytic calculation.
Поведение элементов сложных систем определяется их физико-механическим откликом на силовые, температурные, диффузионные и другие воздействия. Это приводит к постановке связных, в основном нелинейных задач. Перекрестные эффекты определяются феноменологическими составляющими соответствующих систем уравнений, в которых на основе подхода Онзагера для связных физико-механических задач устанавливается связь [1] между обобщенными потоками Jk и обобщенными силами Х^. Jk=LkmXm.
Для деформационных задач эта связь, как правило, устанавливается линейным постулатом Гу-ка. Однако на его основе невозможно описать фазовый переход второго рода от упругости к диссипативному пластическому течению [2], а непрерывные физические величины, такие как температура, концентрация примеси и т.п., связать с условным пределом текучести, поскольку последняя величина является экспериментально-точечной. Как гипотезу, имеющую подтверждение на практике, примем, что коэффициенты взаимности Lkm, постоянные по Онзагеру, могут нелинейно зависеть от обобщенных сил Lhn=Lhn(Xm). Для деформационных задач примем постулат Коши, в соответствии с которым в упругой области зависимость коэффициента Онзагера от деформации в одномерном случае будет следующей.
L(e)=£• (1^е+те2). (1)
Тогда поведение рассматриваемой системы (деформирование образца) полностью определяется потенциалом
S
Y = J E (l - k y+ my2) y d у-ст p e =
Ee ' 2
Eke3 Eme4 - + ■
-ст , e
3 4
и описывается вариационным уравнением
ST = (Ее - Eke1 + Eme3 - стp ) Se = 0 .
Исследования этой системы методами синергетики [3] окончательно приводят к уравнению в упругой области:
ст = E
1 E
1---е+ —
1 E2
Л
3 а„
27 ст„
8= 0.
k=i-E.
3 а„
1 E2
(2)
27 < '
Схематично отклик системы на внешнее воздействие показан на рисунке 1. Здесь к.т. - критическая точка. Если внешнее напряжение выше напряжения в этой точке, происходит хрупкое разрушение. АВ - метастабильное состояние, в котором возможен фазовый переход второго рода с образованием диссипативных структур [2]; ВС -неустойчивый участок, на котором возможен перескок на другую устойчивую ветвь; CD - площадка текучести; АВС - зуб текучести.
Такой подход приводит к решению задач нелинейной упругости. Использование метода конечных элементов требует решения систем большой размерности и увеличенного времени счета из-за трудоемкости итерационного процесса. Возможны несколько путей повышения скорости и точности численного решения. В алгоритме применяются два из них: метод граничных элементов,
0 0,002 0,004 0,006 0,008 0,01 0,012 0,014
Рис. 1. Схема отклика системы на внешнее воздействие. Варианты кривой о-е
Fig. 1. System response on external input
500
400
300
200
100
0
2
e
2
позволяющий уменьшить размерность матрицы разрешающей системы, и максимальное использование аналитических вычислений [4-6].
При таком подходе наиболее затратны по времени два этапа вычислений: построение матрицы разрешающей системы и вычисление интегралов от нелинейных слагаемых в процессе итерационного решения нелинейной задачи.
Интегралы по граничным элементам для построения разрешающей матрицы, в том числе сингулярные, вычислены аналитически. Это позволяет существенно сократить время и повысить точность вычисления матрицы, так как элементы матрицы вычисляются по конечным формулам.
Для вычисления интегралов от нелинейностей используется смешанная схема. Область поиска решения разбивается на конечные элементы, и вычисляются интегралы по элементам. Наибольший вклад в результат дают элементы, содержащие особую точку. Подынтегральная функция на таких элементах терпит разрыв второго рода, поэтому интегралы вычисляются аналитически, что обеспечивает более высокую скорость и точность вычислений по сравнению с обычным численным интегрированием. Интегралы по элементам, не содержащим особую точку, вычисляются по обычным формулам численного интегрирования, так как на этих элементах подынтегральные функции непрерывны, а их вклад в общий результат мал (принцип Сен-Венана) и быстро убывает по мере удаления от особой точки.
Следует отметить, что используемый алгоритм позволяет проводить распараллеливание процесса вычислений.
Построение граничного интегрального уравнения
Рассматривается задача нелинейной упругости в области V с границей Г для однородного изотропного материала. Предполагается, что в упругой области коэффициент Онзагера Ь(е) определяется согласно постулату Коши (1). В многомерном случае Деу) для изотропного тела зависит от инварианта тензора деформаций 12, и его можно представить в виде Ц12)=Е(1-Р(12)), Р(12)=И2--м122, где k и т - положительные постоянные, (2).
Связь тензора напряжений Сту и тензора деформаций примет вид
2С Гов5,, , (3)
= (1 - Р (2 ) ) у
+ и
1 - 2и "
Здесь G - модуль сдвига; и - коэффициент Пуассона; Еу - компоненты тензора деформаций; е - объемная деформация; 5у - символ Кронекера.
В плоском случае 12 = ^(еп -е22 )2 + е^2 .
Коротко можно записать СТу=(1-Р(12))Еукте1<т,
где Е']кт=ь^и8 * 8 кт+° (§л 81т8 *) - тензор
упругих постоянных.
На части границы Г заданы перемещения ы1 (х) = й (х), а на части Г2 заданы силы р1 (х) = р (х). Чертой отмечены заданные функции.
Как известно, метод граничных элементов [7] сводит решение нелинейной упругой задачи для тела V с границей Г к решению интегрального уравнения:
ик (I) = | Р (х)и1 гх - |М, (х)Р* (!> гх +
Г1 Г2
+\ F(12 (г))Ет1гп1 (х)е^ (%, x)dVx - (4)
V
- | М,. (х)Р*к x)dгх + | Р (х)ил гх .
Г1 Г2
Здесь х=(х1, х2) - координаты точки интегрирования; |2) - координаты точки наблюдения; и,(х) - /-я компонента перемещения точки тела; р(х) - /-я компонента усилий; и* (х) - компоненты перемещений; р* (2,, х) - компоненты усилий фундаментального решения. Индекс х у символов <ГХ и «Ух показывает, что интегрирование проводится по переменным х.
Построение дискретного уравнения плоской задачи
Согласно методу граничных элементов [7], разобьем границу области Г на N граничных элементов Г„, в качестве которых в данном случае взяты отрезки прямых. Элементы и узлы нумеруются последовательно, против часовой стрелки. При этом область получается слева, а внешняя часть -справа. Каждому элементу ставятся в соответствие два узла: начало элемента, точка |1л, и его конец, точка |2", в направлении обхода нумерации. Индекс показывает, что эти точки относятся к элементу с номером п.
На элементе перемещения аппроксимируются линейными функциями, усилия - кусочно-постоянными: и" М = и" ф" + и?" ф? р" О) = р" .
Здесь р" - значение /-й компоненты усилия на элементе с номером п; ик - значение /-й компоненты перемещения в узле k элемента с номером п; фП (я) - базисные функции п элемента, Ф1 = 1 - s / Ьп, ф2 (я) = ^ / Ьп; Ь„ - длина элемента. Функции подобраны так, что в одном из узлов она равна 1, в другом 0. Такая система функций называется нормализованной. Вне элемента п на границе функции ф"к (я) обращаются в ноль. Заме-
тим, что функции фП (я) определены только на границе и здесь 5 не вектор, а скалярный параметр, меняющийся вдоль границы. Для вычисления интеграла от нелинейных слагаемых
1Р (I (в)^ в., (х)е; х)с1Ух
v
разобьем область V на М треугольных конечных элементов Vq, q=l, 2, ...,М. Возможны четыре варианта взаимного положения точки и элемента V,! области V.
1. Точка лежит вне элемента. Подынтегральная функция не имеет особенностей, и интеграл по элементу вычисляется по формуле (в плоском случае)
I Р(А (ет ))Ет1 £т1 £**,■ & =
V,
3
= Е Р (12 (£ т ))Ет 1 £ т1 & хк )£/6,
к =1
где - площадь элемента; хк - координаты k узла
£
элемента; ет1 - значения компонент тензора на-
к
пряжений в узле х .
2. Точка совпадает с вершиной элемента Vn. Интеграл становится сингулярным со слабой особенностью. В этом случае интеграл вычисляется аналитически.
3. Точка лежит на границе.
4. Точка внутри элемента.
Случаи 3 и 4 сводятся к случаю 2 разрезанием элемента на вспомогательные.
Запишем (4) для узлов '%„, £=1, 2, ..., N, с учетом непрерывности перемещений и= и1"+1, у=1,
2, где - перемещение £ узла элемента п в направлении j:
" ( 1
< = Е {Ррир (I я) - < РР (I „)} +
p=1
M
+Т J F(12 (в))Eijlmвы (*Кк (Iя, x)dVx, (5)
q= 4 n = 1, 2,..., N.
Здесь Щ = щ £я), Uj (|) = J (I, x) dYx,
r„
P (I) = J p (I, x) ф2 (x) dгx +
r„
+ J p (I, x) фГ (x) dгx .
(6)
В данном случае и далее первый индекс в обозначении перемещений будем опускать. Интегралы (5) по граничным элементам вычислены аналитически. Система (4) решается методом последовательных приближений. Первое приближение [7] ищем, отбросив нелинейные слагаемые:
и,и0 = Z {р 0Up (¡U) - РР )},
p=1
N . .
< = Z {pn U я) - я)}+
p=i
M
+ Z J F(12 (el ))Eijlm80m (x)S; (I„ , X)dV%
q='v
N ( 1
«Г = I {ртир (I я) - «р рг (I я)}+
р=1
м
+Е| р (12 (в;-1 (х)в; (|я, х)^.
Я=1г
ч
Здесь ъг1т(х), Е(12 (е„)) - значения компонент тензора деформаций и, соответственно, функции ^2(в)), вычисленные по г-приближению решения. Верхние индексы показывают номер узла (п) и номер приближения (г). Решение линейной задачи описано в [6]. Блок-схема алгоритма приведена на рисунке 2.
^ Вход ^
Формирование граничных элементов и краевых условий
Построение сетки области
Построение матрицы разрешающей системы и вектора правых частей
Решение линейной системы
Построение поля тензора деформаций
Вычисление интегралов от нелинейных слагаемых
Оценка погрешности
Нет
Коррекция вектора правых частей
Выход
Рис. 2. Блок-схема программы Fig. 2. The algorithm block diagram
г
Пример применения алгоритма
Для иллюстрации работы алгоритма рассматривается задача о растяжении неоднородной изотропной полосы с круговым включением (рис. 3). Длина полосы ¿=100 м, ширина й=10 м, радиус отверстия г=1 м.
Область включения обозначим У(2), соответственно обозначим границы областей: Г(1) - внешняя граница области У(Г), Г(2) - внутренняя граница области и граница включения, - упругие коэффициенты внешней области, Е^ - упругие коэффициенты включения.
Верхний индекс в скобках указывает, к какой области относится параметр или переменная.
На общей границе области V1 и включения выполняются условия непрерывности перемещения: и(1)(х) = и(2)(х), х е Г(2), /-1, 2, и условия равновесия усилий: р(1) (х) = -р(2) (х), х е Г(2), /=1, 2.
Использование изложенного выше подхода приводит к системе нелинейных уравнений:
,,л(1)
N1 / Л
= S {pn 0)p (I я) - <(1) P(1) (|я)} +
р=1
+S {pn(2)ир(Х) (Iя) - up(2)Pjp(1) (|я)} + р=1
M1
+E J F(12 (6(1) ))Epme® (x)e*p (Iя, x)dVx,
q=1 v
uf2 = £ {p(2)Uf(2) и) - uj(2)Pj(2) и)} +
(7)
+Е \ Р(/2 (в(2) »Е^в® (*)в™ (Iя, ,
Я=1У
ч
которая решается методом последовательных приближений.
На рисунке 4 приведены графики нормального ст„ и касательного стх напряжений в сечении х1=0 на отрезке 0<х2<5 относительно растягивающего усилия. Модули сдвига заданы относительно предела текучести стр=70 МПа. Вычисления проводились при значениях относительного модуля сдвига 0(1)=500, 0(2)=250, коэффициент Пуассона и=0,3.
Для сравнения были проведены расчеты при значениях относительного модуля сдвига 0(1)=500,
0(2)=0, что соответствует отверстию, и при 0(1)=500, 0(2)=500, что соответствует однородной полосе. Полученные решения совпадают с теоретическими [8].
На рисунке 4 приведен график касательных напряжений в сечении а-а.
^——_
X
) 0,5 1,5 2 2,5 3
а)
an
б)
Рис. 4. Графики напряжений: а) касательное аТ, б) нормальное а„
Fig. 4. Voltage diagram: a) shearing aT, b) normal an
Таким образом, в данной работе предложена и обоснована физически нелинейная модель упругого тела, позволяющая осуществить корректный фазовый переход (по Пригожину) от потенциального упругого к диссипативному пластическому течению. Построен алгоритм реализаций таких задач, основанный на методе граничных элементов, для чего созданы соответствующие системы граничных интегральных уравнений. В качестве иллюстрации работы алгоритма рассмотрена задача о растяжении неоднородной изотропной полосы с круговым включением.
Литература
1. Пригожин И. Введение в термодинамику необратимых процессов. М.: Иностранная литература, 1960. 127 с.
2. Дьярмати И. Неравновесная термодинамика. М.: Мир, 1974. 301 с.
3. Томпсон Дж.М. Неустойчивости и катастрофы в науке и технике. М.: Мир, 1985. 256 с.
0.12
0
x
0
0
0.5
1.5
2.5
4. Федотов В.П., Нефедова О.А. Применение модифицированного метода граничных элементов для решения параболических задач // Вестн. Сам. гос. техн. ун-та: Сер. Физ.-мат. науки. 2011. Вып. 4 (25). С. 93-101.
5. Федотов В.П., Спевак Л.Ф. Аналитическое интегрирование функций влияния для решения задач упругости и теории потенциала методом граничных элементов // Математическое моделирование. 2007. Т. 19. № 2. С. 87-104.
6. Федотов В.П., Горшков А.В. Численно-аналитический метод решения задач упругости с особенностями // Вестн. Сам. гос. техн. ун-та: Сер. Физ.-мат. науки. 2005. Т. 38. C. 29-34.
7. Бреббия К., Теллес Ж., Вроубел Л. Методы граничных элементов. М.: Мир, 1987. 524 с.
8. Седов Л.И. Механика сплошных сред. В 2 т. М.: Наука, 1983.
References
1. Prigogin I. Introdaction to termodynamics of irreversible processes. Springer Publ., Illinois, USA, 1955, 160 p. (Russ. ed.: Moscow, Inostrannaya literatura Publ., 1960, 127 p.).
2. Dyarmati I. Neravnovesnaya termodinamika [Non-equilibrium thermodynamics]. Moscow, Mir Publ., 1974, 301 p.
3. Thompson J.M.T. Instabilities and catastrophes in science and ingineering. New York, Wiley Publ., 1982, 226 p. (Russ. ed.: Moscow, Mir Publ., 1985, 256 p.).
4. Fedotov V.P., Nefedova O.A. Application of modified method of boundaty elements to solve parabolic problems. Vestn. Sam. gos. tekhn. un-ta: ser. Fiz.-mat. nauki [The bulletin of Samara State Tech. Univ. Physics and Mathematics Science Series]. 2011, vol. 25, pp. 93-101.
5. Fedotov V.P., Spevak L.F. Analytical integration of influence functions to solve the elastic problems and potential theory using the boundary elements method. Matematicheskoe modeliro-vanie [Mathematical modeling]. 2007, vol. 19, no. 2, pp. 87-104.
6. Fedotov V.P., Gorshkov A.V. Numerical-analytical method of solving elastic problem with singularities. Vestn. Sam. gos. tekhn. un-ta: ser. Fiz.-mat. nauki [The bulletin of Samara State Tech. Univ. Physics and Mathematics Science Series]. 2005, vol. 38, pp. 29-34.
7. Brebbiya C.A., Telles J.C.F., Wrobel L.C. Boundary element techniques. Berlin, Springer - Verlaq Publ., Heidelberg New York, Tokyo 1984. (Russ. ed.: Moscow, Mir Publ., 1987, 524 p.).
8. Sedov L.I. Mekhanika sploshnykh sred [Mechanics of continua medium]. Moscow, Nauka Publ., 1983, vol. 1-2.
