Алгоритмы для определения коэффициента диффузии и пористости грунта

Размещено на http://www.allbest.ru/

УДК 532.546

АЛГОРИТМЫ ДЛЯ ОПРЕДЕЛЕНИЯ КОЭФФИЦИЕНТА ДИФФУЗИИ И ПОРИСТОСТИ ГРУНТА

Ч.Дж. Джаныбеков, М.У. Мурзакматов



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

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

а) С целью демонстрации работоспособности алгоритмов, приближенно вычисляющих распространение концентрации загрязнителей, значения сопряженной функции и уточняющих значения коэффициента конвективной диффузии и пористости грунта, обратимся к уравнению диффузии при Di =D(x, y, z):

c = f(x, y, z,t), (x, y, z) V, t (0, T) (1)

с начальным

c(x, y, z,0) = c0(x, y, z), (x, y, z)V (2)

и краевым

lc = (x, y, z, t), (x, y, z) , t(0,T) (3)

условиями, где

n1 = cos(n, x), n2 = cos(n, y), n3 = cos(n, z),

=( v x , v y , v z ), m = m(x, y, z, t), =(x, y, z, t),

=(n1, n2, n3)_ единичный вектор внешней нормали к поверхности .

В начально_краевой задаче (1) - (3) предполагается, что заданные коэффициенты и граница области =V удовлетворяют всем условиям, обеспечивающим существование единственного решения задачи и что функция D(x, y, z) принадлежит пространству Соболева (D(x, y, z)W2(0)).

Рассмотрим случай, когда D(x, y, z) - заданная функция в области V и коротко изложим применение МКЭ при численной реализации задачи (1) - (3) на ЭВМ. С этой целью область V разбиваем на тетраэдральные элементы с общим количеством узлов (вершин тетраэдров) n.

Согласно МКЭ решение краевой задачи ищется в виде

(4)

где ai (t)_искомые коэффициенты, Ni _ линейные базисы МКЭ.

С помощью обобщенного метода Галеркина из краевой задачи (1) - (3) имеем n уравнений

, j =1,2,…,n. (5)

Используя формулу Грина в равенствах (5), приходим к системе интегро_дифференциальных уравнений относительно ai(t):

j =1, 2,…, n. (6)

Интегрируя (6) по t в интервале (tk-1, tk)(при = tk - tk-1,

= T / m, k=1,2,…,m) и используя квадратурные формулы [2]

(7)

получаем систему алгебраических уравнений относительно

aik =ai(tk):

(8)

k=1,2,…,m; j=1,2,…,n,

где

v xk = v x(x, y, z, tk), v yk = v y(x, y, z, tk), v zk = v z(x, y, z, tk),

v nk = v n(x, y, z, tk)= v xk n1+ v yk n2+ v zk n3, k = (x, y, z, tk),

k=(x, y, z, tk), mk = m(x, y, z, tk). (9)

Здесь использованы результаты интегрирования

Учитывая, что в узловых точках (xi, yi, zi)

c(xi, yi, zi, t) cn(xi, yi, zi, t) =ai( t ),

из начального условия (2) имеем n условий (при t=0)

c(xi, yi, zi, 0)= ai( 0 ) =ai,0 , i=1,2,…,n, (10)

которые будут условиями Коши для системы уравнений (6).

Из системы линейных алгебраических уравнений (СЛАУ) (8) при известном векторе =(a10 ,a20 ,…,an0 ), определенном из (10), можем вычислить вектор =(a11 ,a21 ,…,an1 ), тем самым, согласно (4), приближенно построим поле концентрации c(x, y, z, t1). Далее, решая СЛАУ (8) на временном слое t=t2 с использованием , находим c(x, y, z, t2) и т.д. После m шагов получаем поле концентрации c(x, y, z, tm)= c(x, y, z,T).

Несколько слов о решении СЛАУ (8).

При =0.5 СЛАУ совпадает с системой, получаемой из (6) с использованием схемы Кранка_Никольсона (после интегрирования по t в интервале (tk-1,tk)). Поскольку эта система получается путем дискретизации с использованием МКЭ, ее матрица является разреженной. Способ решения системы (8) зависит от знака и величины элементов матрицы, образованных из слагаемых Q ij - q k ij +kij + v n k ij, точнее, от величины qkij . Если модуль вектора скорости фильтрации =(vx, vy, vz) будет значительной, то дифференциальное уравнение (1) по характеру будет приближаться к гиперболическому уравнению и соответственно характер аппроксимирующей ее системы (8) также будет изменяться. Такая ситуация требует более внимательного отношения к выбору метода решения этой системы. конвективный диффузия пористость грунт

б) Приближенное решение сопряженной ретроспективной задачи относительно c*(x, y, z, t) проводится с помощью МКЭ по алгоритму, изложенному в предыдущем пункте. При этом следует отметить, что разбиение области V на элементы и количество узлов сетки остается без изменения, как и в пункте а).

в) На вычислении поля функции c*=c*+c* при t(0,T) следует остановиться более подробно. Функцию с*, при известной функции c*(x,y,z,t), можно вычислить из интегро_дифференциального уравнения (4.10) или (4.13)*), используя соответствующие «начальные» условия. Не занимаясь непосредственным решением уравнения (4.10), предложим более экономичный алгоритм.

С этой целью возвращаемся к возмущенной краевой задаче (4.1), (4.2), откуда относительно возмущенной функции с* получим краевую задачу

c* = p, (x, y, z)V, t(0,T), (11)

l c* =*, (x, y, z) , t(0,T), (12)

с условием в конце отрезка времени

c*(x, y, z, t)|t=T =(x, y, z), (x, y, z) . (13)

Здесь операторы и l определяются согласно (2.5), при этом x1=x, x2=y, x3=z.

Приближенное решение задачи (11) - (13) получим МКЭ. Ее решение ищем в виде

(14)

где ai*(t)_неизвестные коэффициенты, Ni _линейные базисы конечных элементов .

Исходя из обобщенного принципа Галеркина относительно невязок уравнения и краевого условия, получим

(15)

j=1,2,…,n,

откуда на основании формулы Грина имеем

(16)

При t=T , согласно (13), конечными условиями для полученной системы обыкновенных дифференциальных уравнений относительно искомых функций ai*( t ) (i=1,2,…,n), будут

(17)

(Все формулы, имеющие другую нумерацию, взяты из работы [1]).

Приведение системы (16) к системе алгебраических уравнений и реализация ее на ЭВМ проводится аналогично системе (6).

г) Остановимся на проблеме уточнения значений коэффициента диффузии D(x, y, z). С этой целью возвращаемся к возмущенной задаче относительно коэффициента диффузии (5.1), (5.2), которая распадается на две краевые задачи, а именно, на краевую задачу относительно D:



*с* = p, (x, y, z)V, (18)

l*c*=*, (x, y, z) (19)

и на возмущенную задачу относительно D:

*с* = p, (x, y, z)V, (20) l*c*=*, (x, y, z) (21)

при любом t(0,T).

Здесь возмущенные операторы * и l* определены согласно (5.3) - (5.4). Поскольку функция с* определяется как решение задачи (4.1), (4.2) (или (18), (19)) совместно с конечным условием при t=T , то при известной функции с* и заданных функциях р и * из (20), (21) можно вычислить D.

Решение возмущенной задачи (20), (21) построим с помощью обобщенного метода Галеркина с применением МКЭ, для чего приближенное решение ее представим в виде

(22)

и из взвешенных невязок

j=1,2,…,n,

используя формулу Грина, приходим к системе n уравнений относительно искомых коэффициентов Di =D(xi, yi, zi):



или (23)

где j=1,2,…,n .

Из СЛАУ (23) можно вычислить Di , и тем самым, согласно (22), получим приближенное решение

а затем и уточненное значение коэффициента диффузии:

D(1)(x, y, z)=D(0)(x, y, z)+D(x, y, z)=(Di+Di)Ni. (24)

д) Обратимся теперь к задаче уточнения коэффициента пористости грунта при предположении, что m = m(x, y, z, t). Уравнение (5.8), при наличии приближенного значения величины m (т.е. ее начального приближения) в области V, дает возможность уточнить поле пористости грунта. Исходя из возмущенных уравнений (5.11) и (5.12), составляется уравнение, уточняющее поле функции m(x, y, z, t) в области V при t(0,T).

Итак, из (5.11)и (5.12) получаем цепочку уравнений:

краевую задачу относительно функции с*



с*=p, (x, y, z)V, t(0,T), (25)

l*с*=, (x, y, z), t(0,T) (26)

и возмущенное уравнение относительно т

(27)

Поскольку функция с* определяется как решение ретроспективной краевой задачи (25), (26) при дополнительном условии относительно с* при t=T, то из уравнения (27) при известных с* и можно определить m(x, y, z). С этой целью эту функцию представим в виде

(28)

и из уравнения (27) на основании метода Галеркина имеем

(29)

или, после подстановки (28)

(30)

Интегрируя (30) по t в интервале (tk-1, tk) и используя квадратурные формулы (7), приходим к СЛАУ относительно mik :



(31)

где j=1, 2,..,n; k =1,2,…,m; =tk - tk-1=T/m ,

mik=mi(tk), ck*=c*(x, y, z, tk), =(x, y, z, tk),

Таким образом, если известны (начальные) приближенные значения пористости грунта m(0)(x, y, z, t), полученные путем измерения, то, используя математическую модель задачи, согласно (31), построим поле вариаций m(x, y, z ,t), после чего получим возможность уточнить значения коэффициента пористости, т.е. получим следующее приближение

m(1)(x, y, z, t) = m(0)(x, y, z, t)+m(x, y, z, t).

Из вышеизложенного видно, что алгоритмы, полученные из теории возмущений, являются адаптирующимися в том смысле, что в процессе итераций значения коэффициентов D(x, y, z) и m(x, y, z, t) будут постепенно приближаться к реальной ситуации, описывая ее более точно. Ясно, что это возможно только при адекватной математической модели.

е) Совокупность вышеописанных алгоритмов дает возможность построить цельную вычислительную технологию, определяющую интересующие нас параметры.

Предполагается, что с помощью измерения или наблюдения нам известны поля коэффициентов диффузии D(0) и пористости грунта m(0), заданные с определенной погрешностью. Заметим, что на практике никакая из функций D(0)(x, y, z) и m(0)(x, y, z, t) не задается, а вместо них задаются их точечные значения в небольшом дискретном множестве точек. Образование из этих значений функций D(0) и m(0), принадлежащих к определенному классу, является задачей вычислителя.

Вернемся к рассмотрению задачи (1) - (3). Если, кроме искомой концентрации c(x, y, z, t), будут неизвестными еще D(x, y, z) и (или) m(x,y,z,t), то рассматриваемая начально_краевая задача будет существенно нелинейной и не доопределенной, а математически _ некорректно поставленной. Таким образом, кроме некорректности поставленной задачи возникает дополнительное затруднение, связанное с нелинейностью математической модели. Из_за этого, прежде чем образовать итерационную процедуру, мы должны быть уверены, что из многих ветвей решений найдем именно то решение, которое идентифицирует изучаемый нами объект, т.е. из множества решений мы должны выбрать единственно верное.

Для простоты изложения допустим, что функция m(x, y, z, t) задана, а коэффициент D(x, y, z) следует идентифицировать и пусть с помощью измерений или наблюдений с некоторой точностью получены значения функции в некотором множестве точек, т.е. имеем

D(xi, yi, zi)=Diэ , i=1,2,…,m , (32)

где т - натуральное число (т 1020). Считаем, что опорные точки (xi, yi, zi) разбросаны по возможности равномерно по всей области V. Дополнительные условия (32) называются внутренними условиями.

Допустим, что функция D(x,y,z) принадлежит пространству Соболева W2(0), а функцию концентрации ищем в пространстве W2(1). Область V разбивается на тетраэдральные элементы с таким расчетом, чтобы опорные точки (xi,yi,zi), в которых заданы точечные значения функции D(x, y, z), совпадали с узлами сетки. Количество узлов равно n.

Начальное приближение коэффициента D образуем с помощью линейных базисов конечных элементов таким образом, чтобы в опорных точках было D(0)(xi, yi, zi)=Diэ, а в остальных узлах принимает значения из интервала (minDiэ, max Diэ), т.е.



Dn(0)(x, y, z) =Di(0) Ni(x, y, z), (33)

где Di(0)=D(0)(xi, yi, zi), Ni - линейные базисы конечных элементов.

Понятно, что построенная вышеописанным способом функция является интегрируемой в области V.

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

Шаг 1. Принимая за начальное приближение искомого коэффициента значения Dn(0)(x, y, z) из (33) и решая начально_краевую задачу (1) - (3), находим функцию сn(1)(x, y, z,t) - первое приближение искомой функции.

Шаг 2. Решив ретроспективную задачу (2.11) - (2.14) при том же коэффициенте Dn(0) , находим сn*(1)(x,y,z,t)- первое приближение сопряженной функции.

Шаг 3. Решив СЛАУ, вытекающую из (16), находим вариацию с*n в виде (14) и образуем поле сопряженной функции

c(1)* = cn*(1) + cn*. (34)

Шаг 4. Используя найденные значения cn(1) и сn(1)* , вычислим Dn из СЛАУ (23). Вариация коэффициента D ищется в виде (22). После нахождения Dn из (24) образуется коэффициент Dn(1) .

Шаг 5. Подставляя найденные значения коэффициента Dn(1) в начально_краевую задачу (1) - (3), решая ее с помощью обобщенного метода Галеркина с применением МКЭ и выполнив первый шаг, находим cn(2). При значении Dn(1), переходя к шагу 2, находим cn(2)*, и далее, переходя к шагу 3, вычислим cn* , а затем на шаге 4 находим Dn(2) и т.д. Во всех узлах проверяется условие



(35)

где >0 _ заданное малое число.

Теперь рассмотрим случай, когда функцию D(x, y, z) можно считать заданной, а пористость грунта m(x, y, z, t) нуждается в уточнении. Считаем, что в данном случае заданы внутренние условия типа (32) относительно функции m(x, y, z, t). Процесс идентификации этой функции проводится в следующей последовательности.

Шаг 1. При известных функциях c0(x, y, z), D(x, y, z), f(x, y, z, tk), (x,y,z,tk), (x, y, z, tk) и заданной m(0)(xi, yi, zi, tk) (i=1,2,…,n) решаем задачу (1) - (3) и определяем первое приближение концентрации cn(1).

Шаг 2. При тех же значениях пористости mn(0)(tk) из ретроспективной задачи (2.11) - (2.14) находим первое приближение сопряженной функции cn*(1).

Шаг 3. Решая ретроспективную задачу, т.е. систему (16), вычислим значения вариации сопряженной функции сn* и по формуле (34) получаем поле сопряженной функции cn(1)*.

Шаг 4. Используя значения функций cn(1)* и , из системы (31) находим поле вариации m(xi, yi, zi, tk) и образуем первое приближение коэффициента пористости:

m(1) = m(0) +m при t = tk.

Шаг 5. После проверки условия

(35)

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

После каждой итерации проверяется выполнение условия типа (35) относительно функции с(x, y, z, t), а также условия, связанного с данными наблюдений за концентрацией загрязнителей.



Литература

1. Джаныбеков Ч. Дж., Мурзакматов М.У. Идентификация параметров в задаче о переносе загрязнителей в пористой среде // В настоящем номере этого журнала.

2. Самарский А.А. Теория разностных схем.- М.: Наука, 1977.-656с.

Размещено на Allbest.ru

...