Модифицированный метод граничных элементов для решения связных задач математической физики Текст научной статьи по специальности «Математика»

УДК 519.634

МОДИФИЦИРОВАННЫЙ МЕТОД ГРАНИЧНЫХ ЭЛЕМЕНТОВ ДЛЯ РЕШЕНИЯ СВЯЗНЫХ ЗАДАЧ МАТЕМАТИЧЕСКОЙ ФИЗИКИ

В. П. Федотов

Институт машиноведения УрО РАН,

Россия, 620049, Екатеринбург, ул. Комсомольская, 34.

E-mail: fedotov@imach.uran.ru

Предложен модифицированный метод граничных элементов для физико-математического моделирования связных процессов. Физическое моделирование основано на теореме Онзагера связи обобщённых сил и потоков, где предполагается нелинейность коэффициентов взаимности. Подход иллюстрируется на примере диаграммы деформирования. Математическое моделирование основано на модификации метода граничных элементов, основанной на исключении некорректных операций численного дифференцирования и интегрирования.

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

Введение. Развитие методов решения задач математической физики и высокопроизводительных компьютерных комплексов позволяет помимо предварительных расчётов ставить задачу мониторинга сложных систем в условиях реальной эксплуатации.

Поведение элементов сложных систем определяется их физико-механическим откликом на силовые, температурные, диффузионные, магнитные, динамические и другие воздействия. Это приводит к постановке связных, в основном нелинейных задач. Перекрестные эффекты определяются феноменологическими составляющими соответствующих систем уравнений, в которых необходимо учитывать характерные времена и представительные объемы для моделирования связных процессов. Они могут содержать не только искомые величины, но и их производные первого и второго порядков. Это требует больших вычислительных затрат и точности вычислений не только искомых величин, но и их производных минимум до второго порядка, что приводит при численном счете к некорректной задаче.

Для повышения скорости счета и точности вычислений предлагается модифицированный метод граничных элементов (МГЭ) для решения задач эллиптического, параболического и гиперболического типов как базовой основы для решения нелинейных и связных задач математической физики. Модификация базируется на аналитическом вычислении интегралов и производных от функций влияния задач указанных трёх типов и распараллеливании на всех этапах решения задач на уровне алгоритма. Для получения аналитических формул для точного вычисления интегралов от функций влияния и их производных вводится фиксированный «базовый» элемент. Интегрирование компонент функций влияния по произвольному участку границы сводится

Владимир Петрович Федотов (д.т.н., проф.), главный научный сотрудник, лаб. прикладной механики.

к интегрированию по фиксированному «базовому» элементу с последующим линейным преобразованием сдвига и поворота. Интегрирование для всех типов задач производится один раз, в результате чего получены компактные формулы, удобные для программирования.

1. Физическое моделирование. Для мониторинга сложных систем в реальных условиях необходимо знание текущих физико-механических характеристик. Например, для анализа прочности необходимы феноменологические уравнения связи напряжений деформаций, температуры, концентрации легирующих элементов и т.п. Но получение их экспериментальными методами в действующей конструкции трудоемко, а зачастую невозможно. Дальнейшее изложение основано на предположении о невозможности получения физико-механических характеристик и введения понятий устойчивости, разрушения, ресурса и т.п. в рамках замкнутой системы, поскольку они напрямую связаны с характером обмена с внешней средой. Феноменологические соотношения взаимности связных задач могут быть сформулированы в рамках открытых неравновесных систем на основе подхода Онзагера для связных физико-механических задач, согласно которому устанавливается связь между обобщёнными потоками и обобщёнными силами Xд.:

к — Т^тХт.

Как гипотезу, имеющую подтверждение на практике, примем, что коэффициенты взаимности Ь^т-, линейные по Онзагеру, могут зависеть от обобщённых сил = Ькт(Хт). Произведение обобщённых сил Xд. и обобщённых потоков 1к представляет собой скорость роста энтропии Б в изолированной системе [1]. Если система открытая, то полную скорость роста энтропии можно представить в виде двух слагаемых

>5 = ¿>е + = Хе1е + Хг1г, (1)

где Бе = Хе1е — изменение энтропии за счет обмена энергией с внешними системами, Si = Х^Т >0 — изменение энтропии за счет внутренних неравновесных процессов.

С учётом (1) запишем уравнение сохранения энергии:

(¿Л

—гг = То в = ТоХе1е + ТоХгТ = 0, оЬ

где параметр То имеет смысл температуры. Здесь А(ук) называется функцией Ляпунова, или энергетической функцией. В замкнутой системе положительно определенная функция Л(ук) обладает таким свойством, что функция

является знакоотрицательной, а движение системы, соответствующей А(уь), устойчиво, тогда как фазовый переход от упругости к пластичности и разрушение связаны с потерей устойчивости. Поэтому для их описания необходимо рассматривать открытую систему.

Рассмотрим физическое моделирование на примере процесса деформирования. Будем рассматривать деформированный образец как открытую систему. В качестве внутренней обобщённой силы Х^ будет выступать деформация образца бе/То, в качестве внутреннего обобщённого потока ^— внутреннее напряжение а. В качестве обобщённого внешнего потока возьмем внешнее напряжение ар. Отметим, что отождествления внешнего ар и внутреннего а напряжений нет. В этом случае можно сразу записать условие сохранения энергии

(1к

— — = <5Ф = —арбе + аде = —ар5е + Ь(е)бе = 0.

Следуя Коши, будем считать, что в упругой области коэффициент Онзагера зависит от обобщенной силы следующим образом:

L(e) = Е{ 1 - ке + те2).

(2)

Поведение рассматриваемой системы (деформирование образца) полностью определяется потенциалом

Ф = f Е(1 - fry + Ш72)7^7 - аре =

J о

Ее2 Eke3

+

Eme

(3)

и описывается вариационным уравнением

5Ф = (Ее - Eke2 + Eme3 - ар) 5е = 0.

Исследования этой системы методами синергетики [3] окончательно приводят к уравнению в упругой области:

а = Е 1-

1 Е 3 а,

■£ +

\_Е* 27 al

)е = 0, к =

е2 ]е

1 Е 3 а„

m

\_Е* 27 <

Схематично отклик системы на внешнее воздействие представлен на рис. 1. Здесь Q — критическая точка. Если внешнее напряжение выше напряжения в этой точке — хрупкое разрушение. АВ — метастабильное состояние, где возможен фазовый переход второго рода с образованием диссипативных структур [3], ВС — неустойчивый участок, на котором возможен перескок на другую устойчивую ветвь, CD—«площадка текучести», ABC — «зуб текучести».

Для построения модели диссипативного процесса будем постулировать аддитивность процесса дислокационного упрочнения и разупрочнения от образования пространственных дефектов несплошностей [4], т.е. диаграмму а ~ е представим в виде

'в

^—Ч.С D

[А

о

Рис. 1. Упругий участок при различных значениях к и m

а = ат +

а\—ат ,- <72 — ат , ^

е — ет — '

л/1 в — t т

{ев ~ ет)2

(е-етУ

(4)

Здесь <7i и (т2 —значения напряжений, соответствующих деформации £в при отсутствии в первом случае — разупрочнения, во втором случае — упрочнения.

Второе слагаемое правой части отражает известный из физики металлов механизм дислокационного упрочнения, для которого справедлива параболическая связь напряжений и плотности дислокаций р, которая, в свою очередь, линейно связана со сдвиговыми деформациями р ~ е, третье слагаемое является решением линейного кинетического уравнения накопления повре-ждённости.

Для определения а\ и <Т2 используются обычно два очевидных условия:

da

a(sB) = aB] -

О,

а, МПа

500

=£в

400

300

200

100

0

8

Ю е,<

Рис. 2. Диаграмма деформирования стали Х70: 1 — упругость; 2 — упрочнение; 3 — эксперимент

из которых следует <ti = (4<тв — ат)/3, <72 = + 2ат)/Ъ.

Рассмотрим моделирование диаграммы а ~ е стали Х70 при следующих исходных данных: внешнее напряжение ар = 590 МПа, предел текучести ат = = 430 МПа, временное сопротивление а в = 636,75 МПа, деформация временного сопротивления S в = 11,455 %. На рис. 2 показано сравнение модели и эксперимента.

Следует отметить, что при необходимости можно использовать более сложные кинетические уравнения, и это приведет к другому виду последнего слагаемого (4), однако, как показано в [2], приведенная модель (4) даёт достаточное для расчётов приближение к экспериментальным данным для многих металлов при различных температурах и концентраций примесей.

2. Деформационный ресурс и разрушение. Согласно предлагаемой модели поведение образца (материальной точки) определяется внешним воздействием и внутренним откликом материала на каждой стадии. Под разрушением материала в точке будем понимать критическую ситуацию, когда энергия внешнего воздействия сравняется с внутренним потенциальным и диссипа-тивным энергетическим откликом материала. Текущий ресурс в каждый момент деформации будет определяться разностью энергий внешнего воздействия и внутренней реакции материала.

Энергия внешнего воздействия определяется формулой

Wp = [ ар(т)£(т)с1,т. J о

На первой стадии упругий потенциал определён формулой (3), тогда

Её

Фе = [ Е(1 - + Ш72)7^7 =

J о

Ек£3 Етё + —:—•

Для простоты изложения опустим неустойчивое состояние системы с образованием зуба и площадки текучести.

На второй стадии действуют два диссипативных механизма, которые определяются формулой (4), тогда энергетический отклик вычисляется по формуле

а1~аТ .(7_£г)1/2_ <Г2-<ТТ .—— Л 21

Ф«Ш =

^dis

(ев ~ етУ2

(7 - ьт) dq.

-(ев-ет)1/2

На последней стадии разупрочнения, если оно есть, поведение системы определяется третьим слагаемым функции (4), а энергия — интегралом

Фвой =

<72 — <7Т

(ев - £т)2

(7 - ьт)2 dr/

^di,

0,5

0

-0,5

1 \ **

I

0

Текущий остаточный определяется разностью

ресурс системы

Фгев = Wp - Фе - Ф<ш - Ф;

soft. •

(5)

10 £,%

Рис. 3. Деформационный ресурс: 1 — внешняя энергия; 2 — внутренний отклик; 3 — остаточный ресурс

Уравнение Фге8 = 0, или в приведенном виде Ф = Фгез/Щ) = 0, позволяет определить значение деформации efr, при которой произойдет разрушение. Отметим: разрушение может произойти на любой из рассматриваемых стадий, тогда соответствующие слагаемые в (5) отсутствуют. На рис. 3 представлен приведенный ресурс деформирования, соответствующий свойствам материала по модели (4).

Следует отметить, что после а в диаграмма приходит в неустойчивое состояние и для вычисления распределенных характеристик необходимо применение аппарата статистической механики, однако для вычисления интегральных характеристик, входящих в (5), детерминированные зависимости (4) могут быть достаточны для первого приближения.

3. Математическое моделирование. Для решения связных физико-механических задач предлагается модифицированный метод граничных элементов [4,5]. Построение общего алгоритма для таких задач не представляется возможным, поэтому здесь рассматриваются алгоритмы решения базовых задач математической физики, конкретизация алгоритма для каждой реальной задачи требует отдельного рассмотрения.

Эллиптические задачи. Решение линейной упругой задачи определяется интегральным соотношением

Ui(0 = / [u*ij(£,x)f*(x) - f*j(C,x)uj(x)] d,S(x)+

+ [u*ij(C,x)fj(x) - f*j(C,x)u*(x)]dS(x), (6)

J Su

где «*•(£, ж), ж) — тензоры влияния; и^(х), — поверхностные перемещения и напряжения, звездочкой обозначены известные величины. Разбивая границу области на элементы и устремляя внутреннюю точку £ на границу, получим из (6) систему 2р уравнений относительно Щ и

1и-р) = У.(/Г'" / и-Ахур>,х)с1,${х) -и?' I ЦЛхур>,х)с1,${х) | +

2

/< , u*j(x{p)>x)dS(x) ~ иТ /( , f*j{x{p\x)dS(x)

N+M , ч

+ Е f? / (ж(р)'х) dsw -п?* / f*i (xiP)>ж)dSм ■

„_дг I -I \ J J stä /

Здесь звёздочками обозначены известные величины.

Модификация заключается в аналитическом интегрировании тензоров влияния по базовому элементу (рис. 4) с последующим аффинным преобразованием. Влияние любого отрезка границы АВ на точку £ будет эквивалентно влиянию базового отрезка ОС на точку

Х2

Ol l XI

Рис. 4. Базовый элемент

Например,

[ и*Д, x)dS(x) = --U6 (In (Й + Й) - In ((6 - Lf + Щ) , Joe А

Ai = -(8тгМ1-^))"1-

Остальные значения интегралов приведены в [5].

Предлагаемый подход позволяет построить эффективную итерационную процедуру для учёта нелинейности процесса, характерного для связных задач. Физическая модель удовлетворяет требованию фазового перехода второго рода при переходе от упругости к пластичности, что математически исключает угловую точку, производная в которой отсутствует, а задача определения границы раздела упругой и пластической областей является некорректной. Формулировка задачи в обобщенном смысле, т.е. через граничные интегральные уравнения, позволяет учитывать нелинейность только через изменение свободного вектора разрешающей системы, не меняя её матрицы, что существенно экономит время расчёта. Нелинейные свойства можно учесть так:

1 -2v~

2ߣij + -, _ 7./Л/ - = veij ~ <7%,

где согласно (2)

2 E(-kI£ + mI2) , _

ij~ 2(1 +v) l ' 1 2// U U

— начальные напряжения на каждой стадии итерации, 1е — интенсивность деформаций. После стандартных преобразований [5] задача сведётся к системе граничных интегральных уравнений. Входящие в систему интегралы вычисляются аналитически. В итерационном процессе интегралы по поверхности уже вычислены, изменяются только интегралы по области от известных функций.

Гиперболические задачи. В гиперболических задачах помимо граничных условий возникают начальные условия в области для искомой функции и ее производной по времени, что требует вычисления соответствующих интегралов. Применяя к ним второе тождество Грина и проводя преобразования, сокращающие «неберущиеся» интегралы, перейдём к граничному интегральному уравнению

_ , Г*р Г ди(х,т) О -т)Н(фр -т) -г) , ,

*> = -1 /г ВТ ^Л.-гГ-^^-

- Г [у„(х,г) Н^ГТ)~Г) а«Г(»)*г+

Jto JГ

При интегрировании выражений необходимо различать различные случаи расположения волны влияния от точки £ относительно базового элемента. Скорость распространения волны конечна, а значит, влияние точки £ на базовый элемент проявится только после того, как волна, возникшая в точке £ в момент времени т, достигнет базового элемента. Будем различать следующие пять случаев:

1) волна еще не достигла базового элемента; множество таких точек будем обозначать как

Мх = Ц € П : фр-т) -г(£,х) <0, хе [а, Ь}},

где г(£,ж) — расстояние между точками £ и х\

2) волна достигла некоторых внутренних точек базового элемента, но не достигла его границ; множество таких точек будем обозначать в виде

М2 = € П : фр -г)- г(£, а) <о(~]фР -т) - г(£, Ъ) < 0

р| Бх € (0, Ъ) : фр - т) - г(£, х) > о| ;

3) волна достигла только правой границы базового элемента; множество таких точек будем обозначать как

М3 = {С € П : фр-т)-г(£,а) < 0 р| фр - т) - г(£, Ъ) > о} ;

4) волна достигла только левой границы базового элемента; множество таких точек будем обозначать как

М4 = {С € П : c(tF - т) - r(£,a) > of|c(iF - г) - г(£,6) < о} ;

5) волна полностью достигла базового элемента; множество таких точек будем обозначать как

М5 = {{ € П : c(tF - т) - r(i,x) > 0,х € [a, Ь}} .

Аналитические формулы для интегралов по всем положениям точки влияния £ приведены в [4,5].

Для параболических задач построение интегральных уравнений подобно предыдущему, отличие заключается только в аналитическом интегрировании соответствующих функций влияния.

4. Заключение. Решение связных задач математической физики требует рассмотрения двух проблем: 1) разработки физических моделей нелинейных функциональных соотношений связности между обобщёнными потоками и силами, а также физическими параметрами процессов, которые учитывают характерные времена и объемы материальных точек; 2) разработки математических методов и алгоритмов, позволяющих решать громоздкие системы дифференциальных уравнений в частных производных достаточно быстро и обеспечивающих высокую точность вычислений не только искомых функций, но и их первых и вторых производных. В работе предложен подход, базирующийся на неравновесной термодинамике открытых систем и позволивший не только рассматривать связные процессы деформирования, диффузии, теплопроводности, но и отслеживать состояния неустойчивости систем с определением фазовых переходов, деформационного ресурса и разрушения. Модифицированный метод граничных элементов за счёт предварительного аналитического счёта позволил существенно сократить вычислительное время и за счёт точного вычисления производных и интегралов от функций влияния дают возможность исключить некорректные операции численного дифференцирования и интегрирования. Непрерывные в области решения дают возможность исключить сеточное и конечно-элементное дробление самой области, что вызывает трудности для связных процессов с различными характерными представительными объемами. Следует отметить, что предлагаемый подход пока нельзя назвать универсальным, но с его помощью в единой методологии можно рассматривать отдельные классы связных задач.

Работа выполнена при финансовой поддержке Программ Президиума РАН № 15 (проект № 12—П—1 — 1031 «Математическое моделирование нелинейных связных физико-механических процессов на основе модифицированного метода граничных элементов») и № 18 (проект № 12—П—1 — 1028 «Алгоритмы и математическое обеспечение для вычислительных систем сверхвысокой производительности»)

БИБЛИОГРАФИЧЕСКИЙ СПИСОК

1. I. Gyarmati, Non-equilibrium Thermodynamics. Field Theory and Variational Principles.

Berlin: Springer, 1970. xi+184 pp.; русск. пер.: И. Дьярмати, Неравновесная термодинамика. Теория поля и вариационные принципы. М.: Мир, 1974. 301 с.

2. J. М. Т. Tompson, Instabilities and catastrophes in science and engineering. Chichester: Wiley, 1981. xvi+242 pp.; русск. пер.: Дж. M. Томпсон, Неустойчивости и катастрофы в науке и технике. М.: Мир, 1985. 256 с.

3. I. Prigogine, Introduction to thermodynamics of irreversible processes, 3d ed. New York: Interscience Publishers, 1967. xv+147 pp.; русск. пер.: И. Пригожим, Введение в термодинамику необратимых процессов. М.: Иностр. лит., 1960. 127 с.

4. В. П. Федотов, А. А. Контеев, "Модифицированный метод граничных элементов для задач о колебаниях плоских мембран"// Тр. ИММ УрО РАН, 2009. Т. 15, №2. С. 211 221; англ. пер.: V. P. Fedotov, A. A. Konteev, "Modified boundary element method for problems on oscillations of flat membranes" // Proc. Inst. Math. Mech., 2009. Vol. 267, no. 1. Pp. S78-S89.

5. В. П. Федотов, Л. Ф. Спевак, Решение связных диффузионно-деформационных задач на основе алгоритмов параллельного действия. Екатеринбург: УрО РАН, 2007. 172 с. [V. P. Fedotov, L. F. Spevak, Solution of Coupled Diffusive-Deformation Problems on the Basis of Parallel Algorithms. Ekaterinburg: UrO RAN, 2007. 172 pp.]

Поступила в редакцию 02/XI/2012; в окончательном варианте — 24/XII/2012.

MSC: 65М38, 74S15; 35Q74

MODIFIED BOUNDARY ELEMENT METHOD FOR THE SOLUTION OF CONNECTED PROBLEMS OF MATHEMATICAL PHYSICS

V. P. Fedotov

Institute of Engineering Science, Ural Branch of RAS,

34, Komsomolskaya St., Ekaterinburg, 620049, Russia.

E-mail: fedotov®imach.uran.ru

The modified boundary element method for physico-mathematical modeling of the multifactorial processes is offered for discussion. Physical modeling based on the Onsager's theorem about the relationship between generalized forces and fluxes where we assume the coefficients of reciprocity are nonlinear. The approach is illustrated by the example of the strain diagram. Mathematical modeling is based on a modification of the BEM where all incorrect procedures of the numerical differentiation and integration were replaced by preliminary analytical calculations.

Key words: phase transition, open systems, stresses, strains, deformation resource, integral equation, influence function, elliptic and hyperbolic problems.

Original article submitted 02/XI/2012; revision submitted 24/XII/2012.

Vladimir P. Fedotov (Dr. Sei. (Techn)), Chief Researcher, Lab. of Applied Mechanics.