Одномерная обратная задача диффузии-конвекции Текст научной статьи по специальности «Математика»

Раздел V. МЕХАНИКА И МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ

Э.А. Гребещенко

ОДНОМЕРНАЯ ОБРАТНАЯ ЗАДАЧА ДИФФУЗИИ-КОНВЕКЦИИ

Аннотация. Рассмотрены и разобраны аналитическое и численные решения уравнений диффузии и конвекции с последующей реализацией в среде пакета Mathcad.

Ключевые слова: обратная задача, аналитическое решение, численное решение, распределение, диффузия, конвекция.

E.A. Grebeshenko

ONE-DIMENSIONAL INVERSE PROBLEM OF DIFFUSION-CONVECTION

Abstract. Considered and dealt with analytical and numerical solution of the equations of diffusion and convection with the subsequent implementation in the environment of Mathcad.

Keywords: inverse problem, analytical solution, numerical solution, distribution, diffusion, convection.

Аналитическое решение уравнения диффузии

Рассмотрим первую начально-краевую задачу для неоднородного уравнения теплопроводности вида:

м; = ku'x +/,0 < x < 1,0 < t < т. (1)

u (x, 0) = u0( x), u (0, t) = 0, u (l, t) = 0. (2)

В отношении функций u (x,t) и /(x,t) будем предполагать их представимость в виде рядов Фурье по тригонометрическому базису:

i(x, t) = £ Cm' (t )sin(®mr), / (x, t) = £ Cf (t) sin(amx),

m-1 m-1

2 / 2 / где ю = ж //,С^ = /(x)sin(юmx)dx,С^ = — /(x)sin(юmx)dx.

Наряду с задачей (1) - (2), которую будем называть «точной» непрерывной задачей, рассмотрим первую начально-краевую задачу с «усеченной» правой частью

/("'( х, t) = £ С»!' (t )81п(ютх),

m-1

и «усеченными» начальными условиями

N-1

м(") (х, t) = £ С"' (0)sin(юmx)

т-1

для задачи вида

("N')' = k ("N')" + /(л'' ,0 < х < 1,0 < г < Т.

"(N' (х, 0) = "0и) (х),"N' (0, г) = "N' (/, г) = 0. Необходимо отметить, что такая постановка возникает, например, в случае табличного способа задания функции/ и и0 ряды будут ограничены N - 1 гармониками т.к. для восстановления непрерывной функции применяется интерполяционный тригонометрический полином, где N — количество дискретных значений функции. Функции и(м) и ^ при подстановке в (1) приводят к соотношению

С N-1 V С N-1 V N-1

I £ С" sin(юmx) I = к I £ С»") sin(юmx) I + £ С{/> sin(юmx).

V m-1 )г V m-1 ) ж п-1

Меняя очередность операции дифференцирования и суммирования ряда, получим:

N-1 ' N-1 N-1

YJC:) (t)) sin(fflmx) = £kC(m} (-ю2m2 sin(®mx)) + £Cf sin(amx).

m-1 t m-1 n-1

Принимая во внимание линейную независимость функций sin(wmx), получим:

(CI" )(t) )t = -ka>2 mC> + C

Решение уравнения примет вид:

cm )(t)=

c(u) -.

C

(f) Л

ka2

C

e

(f) m

ka2 m2

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

N —1 ((

и ( X, t ) = £

m-1

C(и) —

m,0

W

C ( f ) m

ka2m2

\

c (f) m

ka2 m2

sin(amx).

Численное решение задачи конвекции

Рассмотрим уравнение переноса

dq dq — + и — = 0 , dt дх

где t е[0,Г], хе[0,L], q(0,х) = q0 (х), q(t,0) = q(t,L) = 0, и = const.

Введем равномерную расчетную сетку a = ah xaz, где ah = {х..|х( = ih, i = 0,1,..., N, Nh = Lj,

a = {ф = 0,1,..j, т = tn+i — tn = const.

Для численного решения поставленной задачи можно использовать следующие конечно-разностные схемы [1]

- схема «против потока» или при и > 0 левый уголок [1]

(3)

n+1 n n n

q — q + иq — qi—1 = 0

(4)

- центральная разностная схема

n+1 n n n

q — q + и^ — qi—1 = 0; т 2h

(5)

- схема «кабаре» [2]

nn

n+1 n n n—1

qn — q,_ + qi—1 — qi—1 + uqi — qi—1 = 0.

(6)

2т 2т h

Для решения задачи конвекции можно использовать схему, полученную в результате линейной комбинации центральной разностной схемы (5) и схемы «кабаре» (6)

п+1 п п п-1 п . л п г п

д -д + д,-1 -+ иям + Ч-4-1 = 0. (7)

т 2т 4h

Найдем численные решения на основе схем (4) - (7) модельной задачи с начальными условиями д0 (х) = Н(20 - х), где к(х) - функция Хэвисайда. На рис. 1 приведены решения модельной задачи

(1 - численное решение, 2 - точное). 1.5Г

0.5-

-0.5L

Л Л

V

50

100 150

2 ...2

ka m t

m

2 2

ka m t

e

h

z

1.5

0.5

-0.5

2

1 ЛУГаял

50

100

150

1 5г

0.5"

-0.5L

1 Л, 2

50

100

15

Рис. 1. Решения модельной задачи на основе схем (2) - (5)

а) левая разностная схема; б) центральная разностная схема; в) схема «кабаре»; г) предложенная

разностная схема

Из рисунка1 видно, что левая разностная схема аппроксимирует «скачок» недостаточно точно, центральная разностная схема неустойчива, решение, полученной на основе схемы «кабаре» имеет осцилляции (энтропийные возмущения) [3]. Предложенная разностная схема дает наиболее точное решение модельной задачи. Левая разностная схема имеет порядок аппроксимации О (т + к), центральная разностная схема аппроксимирует непрерывную модель с порядком

О (т + к2), схема «кабаре» - с порядком О (т2 + к2).

Далее будут рассмотрены численное решение задачи диффузии на основе явных схем и на основе схем с весами (неявных) в плане реализации данных задачи в среде пакета Mathcad.

Численное решение задачи диффузии на основе явных схем Будем рассматривать первую краевую задачу для уравнения теплопроводности с постоянными коэффициентами:

дТ д 2Т

— = к — + Пх,0,0 < х < 1х ,0 < х < I, (8)

дt дх

где к - коэффициент температуропроводности, со следующими начальным и граничными условиями:

Т(х,0) = щ(х) , (9)

Т(0,,) = Ц1(,), Т(1х,,) = Ц2«. (10)

Определим равномерную сетку юкт с шагом к по пространственной перменной и шагом т по временной переменной. Для сеточной функции Т(х,,) введем обозначение Т"=Т(х-,,п), где - = 0..^х (кх^ = 1х), п = 0...N1 N к, = I,). Правую часть заменим приближенно сеточной функцией ф". Явная разностная схема для уравнения (8) будет выглядеть следующим образом:

Т. -Тп Т" - 2Тп + Т."

■ = к-!+---^ + с

к

к

(11)

(12)

Уравнение (11) решается по слоям, соответствующим моментам времени. Если решение найдено на слое п, то решение на слое п+1 вычисляется по явной формуле. Схема (11) устойчива только при условии

7 = кк < 1 7 К < 2'

Рассматривается одномерная задача диффузии

5Т , д2Т ~дг

с начальным условием:

[1, х е [40,60] [0, х ^ [40,60] и граничными условиями:

Т(0,,) = 0, Т(1х,,) = 0.

Шаг по пространственной переменной будет равен кх = 1, шаг по времени к = 0,1. Ниже приведен фрагмент программной реализации решения уравнения диффузии на основе явной схемы:

— = к — + / (х,,), к = 0,0 < х < 100, 0 <, < 100

дх2

( х >={;

Рис. 2. Программный блок, предназначенный для решения задачи диффузии

на основе явной схемы

На рис. 3 приведено начальное распределение поля температур (слева) и решение задачи диффузии при временных интервалах равных 20,50,100 (справа).

Рис. 3. Решения задачи диффузии на основе явной схемы а) начальные условия; б) график решения задачи на основе явной схемы

Численное решение задачи диффузии на основе схем с весами

Чисто неявной схемой для уравнения (8) называется схема вида

Т^+1 т

грП+1 г\ грП+1 грП+1 !- = кТ+1 - 2Т + Т-1 + «

К К

(13)

Данная схема также решается послойно; и на каждом (п+1)-ом слое приводит к трехдиагональной системе с количеством неизвестных (К - 1). Схема (13) имеет первый порядок аппроксимации по и второй порядок по кх.

п

Рис. 4. Программный блок, предназначенный для решения задачи диффузии на основе а) неявной схемы; б) схемы Кранка- Николсон

Запишем для уравнения (8) разностную схему, представляя вторую производную функции Т по координате х в виде двух слагаемых

5 2Т = 1 а^Т 1 5 2Т Эх2 " 2 Эх2 2 Эх2

и аппроксимируя первое слагаемое на п-ом шаге по времени, а второе - на (п + 1)-ом шаге по времени:

Т«+1 тп

К

грП+1 г\грП+1 . гр П+1 грп г\грп . грп

= кТП1 - 2Т-- + Т-1 + кТП - 2Тп + Т-1 + ,

2К2

2К2

(14)

Данная разностная схема называется разностной схемой Кранка-Николсона ,названа в честь создавших её авторов. Разложение второй производной функции Т по координате х на две равноценные составляющие, одна из которых аппроксимируется на п-ом шаге по времени, а другая - на (п + 1)-ом шаге по времени, указывает на то, что аппроксимацию этой производной в целом следует рассматривать на шаге по времени (п + 1/2). В то же время конечная разность, аппроксимирующая производную функции и по времени, по отношению к точке (п + 1/2) является центральной конечной разностью, имеющей, как известно, второй порядок аппроксимации. Следовательно, разностная схема Кранка-Николсона аппроксимирует уравнение (8) со вторым порядком и по времени, и по координате.

На рис. 5 приведено решение задачи диффузии на основе неявной и схемы Кранка- Николсона при этом временной интервал равен 100 с.

О 20 40 60 80 100

Рис. 5. Сопоставление решения задачи диффузии на основе неявной и схемы Кранка- Николсона

Численное решение обратной задачи диффузии

Будем рассматривать некорректную задачу диффузии с обратным временем получаемую из соответствующей прямой задачи заменой t на Ч (т.е. переходом к обратному времени)

5Т __k &T

dt ~ dx2 ' Т(х,0) = u0(x) , T(0,t) = ^(t), T(4,t) = ^2(t>.

(15)

(16) (17)

Для приближенного решения некорректной задачи (15) - (17), будем использовать следующее уравнение

дТ , 52T 54T — _ _k—г _а—-dt дх2 дх4

(18)

Данный подход был предложен в работе [1]. Уравнение (18) может быть записано в следующем виде:

дТ д2 ( д 2Т

— _—т\ _kT _а—-дt дх2 I дх2

(19)

Для удобства последующей дискретизации уравнение (5) запишем в виде следующей сис-

темы уравнений:

дТ д2 s д 2Т

— _ —г , s _ _kl _ а —-

дt дх2 дх2

(20)

Дискретные аналоги уравнений (20) с запишутся в следующем виде:

ТП+1 _ ТП , s^ _ 2sn + s\Г, , „ Тм _ 2Т + Т_1

h

_ k-

2h2

-±± s __кТ _а-

h2

Ниже приведен фрагмент программной реализации решения уравнения обратной диффузии на основе явной схемы:

Рис. 6. Программный блок, предназначенный для решения обратной задачи диффузии

В качестве начального распределение поля температур было взято решение задачи диффузии при временных интервалах равных 20,50,100 (справа), приведенное на рисунке 3,ниже представлено распределение на основе восстановленного времени.

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

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

1. Головизнин, В.М., Самарский, А.А. Разностная аппроксимация конвективного переноса с пространственным расщеплением временной производной // Математическое моделирование. 1998. T. 10. № 1. - С. 86-100.

2. Ладонкина, М.Е., Неклюдова, О.А., Тишкин, В.Ф. Использование разрывного метода Галеркина при решении задач газовой динамики // Матем. Моделирование. 2014. № 26:1. - С. 17-32

3. Марчук, Г.И., Математическое моделирование в проблеме окружающей среды. М.: Наука, 1982.- 320 с.

4. Самарский, А.А., Вабищевич, П.Н. Численные методы решения задач конвекции-диффузии. - М.: Эдиториал УРСС, 1999.- 261 с.

5. Сухинов, А.И., Белова, Ю.В. Математическая модель трансформации форм фосфора, азота и кремния в движущейся турбулентной водной среде в задачах динамики планктонных популяций // Инженерный вестник Дона, 2015, Т. 37, № 3, - 50 с.

6. Сухинов, А. И., Чистяков, А. Е., Якобовский, М. В. Точность численного решения уравнения диффузии-конвекции на основе разностных схем второго и четвертого порядков погрешности аппроксимации// Вестн. ЮУрГУ. Сер. Выч. ма-тем. информ., 2016, Т. 5, № 1. - С. 47-62

7. Лапин, Д.В., Чистяков, А.Е., Сухинов, А.А. Численное решение прямых и обратных задач диффузии-конвекции на многопроцессорных системах для прогноза и ретроспективного анализа водных экосистем // Известия ЮФУ. Технические науки. 2014. № 12 (161). - С. 230-242.

А.А.Илюхин, А.С. Гондаревская, Л.А.Николаева

ВЫРОЖДЕНИЕ СИСТЕМЫ ИНТЕГРАЛОВ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ ДВИЖЕНИЯ ГИРОСТАТА

Аннотация. Преобразования уравнений движения в зависимости от цели преобразования в той или иной форме используют [4,5,6,7,8] общие интегралы уравнений движений. При этом независимость системы интегралов не подвергается сомнению. Однако, в работах [1,2,3] указана ситуация, при которой такое утверждение просто неверно. В данной работе выписан явный вид уравнений особого многообразия, на котором система интегралов вырождена.

Ключевые слова: уравнения движения, общие интегралы, особые решения.

A.A. Ilyukhin, A.S. Gondarevskaya, L.A.Nikolaeva

THE DEGENERATION OF THE SYSTEM OF INTEGRALS OF THE DIFFERENTIAL EQUATIONS OF MOTION OF A GYROSTAT

Abstract. Transformation of equations of motion depending on the goal of transformation in one form or another is used [4,5,6,7,8], the General integrals of the equations of motion. The independence of the system of integrals is not questioned. However, in [1,2,3] indicated the situation in which this statement is simply wrong. In this work written out the explicit form of the equations of special manifolds, in which the system of integrals is singular.

Keywords: the equations of motion, general integrals, particular solutions.

Основные уравнения. Запишем уравнение движения гиростата и интегралы этой системы дифференциальных уравнений в векторной форме:

d(М + Х) + ^ __ (М2 + я2)(м3 + Х = (e2r3 -еъу2) •р , dt

d(М2 + ¿2) + («3 _ai) • (М_3 + Я3)(А_1 + Х) = (еу -еууъ) • P , (1)

dt

d (М 3 + Я3) dt

+ (а1 - а2) • (М1 + Х)(М2 + Х2 = (еу2 - е2у) • P .

^ = «3 • (М3 + ¿3) У _ «2 • (М2 + ¿2) У , dt

^ = а1 • (М1 + Х1)• у3 -а3 • (М3 + Х3)• у1, (2)

dt

dy3 = «2 • (М 2 + Х 2) у - «1 • (М1 + Х1) • У2 . dt

У12 + У22 + У32 = 1 , (М1 + Х1) • у1 + (М2 + Х2) • у2 + (М3 + Х3) • у3 = k , (3)

1 ' ' - \2 , „ /ЪЛ ,1 \2 , „ /ЪЛ , - \2

—ахх • (М1 + Х1)2 + • (М2 + Х2)2 + • (М3 + Х3)2 | + +«23 • (М2 + Х2)• (М3 + Х3) +

+ а31 • (М3 + Хз)• (М1 + Я:) + а12 • (М1 + Я:)• (М2 + Х2)~(е1у1 + е2у2 + е3у3)• Р = 2h . Здесь М + Х - вектор кинетического момента гиростата, Х - гиростатический момент, ^ -вектор угловой скорости тела - носителя, у - единичный вектор вертикали, т - единичный вектор