CC BY
22
УЧЕНЫЕ ЗАПИСКИ КАЗАНСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА
Том 150, кн. 3
Физико-математические пауки
2008
УДК 532.546
РЕШЕНИЕ МЕТОДОМ ИТЕРАЦИЙ ЗАДАЧИ НАСЫЩЕНИЯ ДВУСЛОЙНОГО ПОРИСТОГО МАТЕРИАЛА
М.М. Алимов
Аннотация
Рассматривается плоская задача определения стационарной конфигурации межфаз-пой границы насыщения двуслойного пористого материала жидкостью. Для случая большого контраста толщин слоев и их пропицаемостей она формулируется как задача со свободной границей в толстом малопропицаемом слое со специфическим условием па границе слоев. Для ее численного анализа предлагается итерационный метод, аналогичный методу Леви Чивиты в теории струй идеальной жидкости. С его помощью в некотором диапазоне изменения параметров процесса с удовлетворительной точностью удалось получить числешго-апалитическое решение задачи. Проведено сравнение результатов с известными асимптотическими оценками.
Ключевые слова: задачи со свободными границами, многофазные среды, насыщение жидкостью пористых материалов, технология УАКГМ.
Введение
Пусть имеется горизонтальный слоистый пласт, кровля и подошва которого непроницаемы. Пласт состоит из двух слоев: нижнего толстого и верхнего тонкого. Слои имеют толщину Н±, пористость т± и проницаем ость к±, где все величины постоянны, верхний индекс плюс относится к тонкому, а минус к толстому слою. Между слоями возможны перетоки. Имеется межфазная граница некоторая поверхность, соединяющая кровлю и подошву пласта и разбивающая весь пласт на две зоны: справа от границы пласт заполнен воздухом, слева полностью насыщен водой (вязкой жидкостью). Воздух справа на бесконечности соединяется с атмосферой, и весь пласт первоначально поддерживается при постоянном атмосферном давлении. Затем слева из бесконечности по всей мощности пласта с постоянным расходом начинает подаваться вода, которая вытесняет воздух в атмосферу. По истечению некоторого промежутка времени устанавливается стационарная конфигурация границы между фазами и скорость ее продвижения. Основная задача найти эту стационарную конфигурацию межфазной границы и величину скорости ее продвижения в зависимости от параметров пласта.
Задача интересна с точки зрения как возможных геофизических приложений (для анализа процессов совместного движения воды и воздуха в слоистых природных пластах), так и возможных технологических приложений. В частности, она возникает в производстве композитных материалов по технологии вакуумного вливания смолы в закрытые формы (УАИТМ). когда к основному образцу-слою пористого и малопроницаемого материала с целыо существенного ускорения процесса его насыщения полимером добавляется тонкий «распределительный» слой высокопроницаемого материала [1. 2].
1. Постановка задачи
Основные предположения, касающиеся характера процесса насыщения пористого материала, состоят в том, что вязкая жидкость несжимаема, полностью насыщает пористую матрицу (нет двухфазных зон), движение ее подчиняется закону Дарен, силы тяжести отсутствуют. Соответственно, при условии установившейся формы межфазной границы в системе отчета, связанной с этой границей, движение жидкости будет установившимся и потенциальным.
В описываемом процессе существенны два безразмерных комплекса
к+Ь+ ш~Ь~
Первый характеризует соотношение между пропускными свойствами слоев, второй между емкостными свойствами слоев. В технологических приложениях наиболее интересен случай сильного контраста толщин слоев (5 С 1) и их проницаемостей (е С 1). Как показано в работе [3], в этом случае задача определения стационарной формы межфазной границы насыщения двуслойного пористого материала сводится к краевой задаче для одного, а именно для толстого малопроницаемого слоя. В безразмерной постановке в терминах потенциала у скорости течения она формулируется следующим образом
ЛБ : ду/ду = 0; (2)
АС- {&Р
дх2 ду ' \дх' ду
(1,0); (з)
БС : <¿> = 0, ^г-=ивт/3, (4)
дп
где О г - область физической плоскости, отвечающая смоченной части толстого малопроницаемого слоя. Ее граница состоит из трех участков: ЛБ - непроницаемой подошвы пласта, ЛС - границы между слоями и Б С - межфазной границы
Л
ется скорость течения (см. вторую часть условия (3)). Через в обозначен угол наклона к горизонту касательной к контуру БС, через и - безразмерная скорость продвижения контура БС, которая полностью выражается через параметры е и 5
1 + е 15
Граничное условие (3) характеризует приток жидкости из верхнего тонкого высокопроницаемого слоя и после интегрирования вдоль границы может быть приведено к виду [3]
АС: || = 1+е[1-</>], фА = 1, Фс = и, (5)
где ф - функция тока, гармонически сопряженная у.
Граничное условие на свободной границе (4). как обычно, для задач типа Хеле-Шоу состоит из двух частей. Первая часть динамическое условие, которое означает. что действием капиллярных сил на межфазной поверхности иренебрегается. Вторая часть кинематическое условие материальности свободной границы. Условие (4) также можно переписать в двух эквивалентных формулировках [4 6]
DC : у = 0, ф = u(y +1); (6)
ду и\2 ( ду\2 и2
DC: * = {-à-2) =Т- (7)
В результате сравнения (5) и (7) можно найти точное значение вектора скорости жидкости v = Vy в точке C
vx,c = 1 + е(1 - и) « S, —l'y,с = vx,c («< - vx>c) « \[Щ~£
и, соответственно, угол наклона этого вектора к оси у
/Зс = \1>х,с/г'у,с \ ~
Отметим, что краевая задача (1) (4) существенно нелинейна. Нелинейность постановки в физической плоскости обусловлена наличием свободной границы БО. Результативный в некоторых случаях прием ухода от такой нелинейности с помощью переформулировки задачи во вспомогательной плоскости здесь не так эффективен, поскольку при этом становится нелинейным граничное условие на участке АО.
2. Схема построения численно-аналитического решения задачи
Поскольку течение потенциально и имеется свободная граница, естественно воспользоваться общей идеологией метода годографа скорости из теории струй идеальной жидкости. В первую очередь надо провести параметризацию решения ввести вспомогательную плоскость комплексного переменного, скажем, £, в результате чего краевая задача сведется к отысканию вида двух аналитических функций комплексного переменного, например:
г = г(С), Ш = Ш (С), (8)
где г(() - функция, конформно отображающая вспомогательную комплексную плоскость £ на физическую комплексную плоскость г, а Ш(() - функция, конформно отображающая вспомогательную комплексную плоскость £ на плоскость комплексного потенциала течения Ш = у + гф.
Как правило, в теории струй в качестве вспомогательной плоскости выбирается область канонического вида так, чтобы по крайней мере одна из функций параметризованного решения определялась сразу в замкнутом виде, а вторая представлялась в виде, удобном для численного анализа (если она не определяется сразу).
С (ш)
А
1
ш
А
-о
Б (0)
Рис. 2. Вид плоскости комплексного потенциала течения Ш
В данном случае, к сожалению, ни одна из функций (8) сразу не определяется и надо ориентироваться на численную процедуру поиска решения. При этом, выбирая в качестве вспомогательной плоскости £ область канонического вида, только для одной из функций параметризованного решения, скажем Ш(£), удается добиться представления, удобного для численного анализа. Для другой функции естественно ввести другую вспомогательную плоскость комплексного переменного, скажем п, чтобы функция г(п) представлялась в виде, удобном для численного анализа. Между двумя вспомогательными плоскостями С ш п канонического вида можно установить взаимнооднозначное соответствие п = п(С) > которое будет выражаться замкнутой аналитической формулой. В результате, помимо параметрического вида решения (8), можно использовать и следующий вид
Соответственно, из четырех функций (8), (9) достаточно определить численно две функции Ш) и г(п), предварительно представив их, например, в виде функциональных рядов. В частности, при подходящем выборе вспомогательных плоскостей £ и п функции Ш(£) и г(-ц) в определяющих уравнениях можно представить в виде рядов Фурье. Тогда появляется возможность использования быстрого преобразования Фурье (далее БПФ).
В то же время качественный анализ задачи (1) (4) показывает, что функции Ш(С) и г(п) обязательно будут иметь дробно-степенные особенности в точке С (есть также логарифмические особенности в точке Л, но они не создают больших трудностей). Действительно, области Ог и О щ представляют собой треугольники ЛСБ с одной криволинейной и неизвестной по форме стороной: у области Ог это контур СБ (см. рис. 1), у области О щ это граница ЛС (см. рис. 2). В области Ощ известно положение всех угловых точек Л, С, Б, в то время, как в области Ог координата хп неизвестна и подлежит определению в ходе решения задачи. Все углы в областях О2 и О щ известны, в частности, угол при точке С в обеих областях составляет рс = п/2.
Как известно, наличие дробно-степенных особенностей приводит к плохой сходимости рядов Фурье н делает практически невозможным определение функций г(п) и Ш(£) непосредственно в виде таких рядов. Поэтому предлагается, аналогично методу Лови Чивиты в теории струй идеальной жидкости [7], аналитически в замкнутом виде выделять главную часть решения, содержащую нужные особенно-С
г = г(п), Ш = Ш (п).
(9)
г(п) = го(п) + гх(п), Ш (С) = Ш„(0 + Ш^С),
А (г)
СИ)
Б (0)
С (1)
А(0)
Б (1)
Рис. 3. Вид вспомогательных плоскостей: о) переменного С, б) переменного п
где и ^о(^) — функции известного вида, являющиеся главными частями
решения; а г^ц) и Wl(Z) — гладкие (в смысле отсутствия особенности в точке С) добавки неизвестного вида. Для их численного определения из краевой задачи (1) (4) будут выведены два нелинейных граничных интегро-дифференциальных уравнения.
В то же время есть одно существенное отличие от упомянутого метода Лови Чивиты. В данном случае невозможно путем локального анализа полностью выявить характер дробно-степенных особенностей функции г(п) и W(С) в точке С
тель при особенности, иначе ее момент, остается неопределенным (логарифмические особенности в точке А таких проблем не создают, поэтому на них внимание не акцентируется). Следовательно, можно ожидать, что получить решение методом итераций только функций г1(п), W1(Z) при предварительно построенных функциях г0(п), W0(Z) не удастся. Последние также необходимо перестраивать по ходу итераций.
Действительно, пробные расчеты показали, что при фиксированном виде функций го(п), Wо(Z) и наложении всех необходимых условий на искомые функции г1(п), W1(C) в их разложениях Фурье уже с первых итераций проявляется старшая гармоника возмущений, что неизбежно приводит к расходимости итерационного процесса.
В то же время успешной оказалась другая стратегия отказаться от некоторых необходимых условий на функции г^ц) и Wl(Z), а появившуюся свободу использовать для сглаживания искомых функций на каждой итерации, фактически для принудительного погашения старшей гармоники. Недостающие необходимые условия накладываются уже на полные функции г(п) и W(С) и удовлетворяются за счет перестраивания главных частей решения г0(п) и W0(Z). Таким образом, предлагается организовать итерационный процесс, каждая итерация которого состоит из двух шагов: па первом шаге определять гладкие добавки (п), Wlj) (С), а па втором - перестраивать главные части решения г^ (п), W д*1 (С), где верхний индекс (]) здесь и далее обозначает номер итерации.
В качестве вспомогательных плоскостей комплексного параметрического переменного С и п выберем области канонического вида — четверти единичного круга с соответствием точек, представленным на рис. 3, о, б.
Взаимнооднозначное соответствие между этими плоскостями устанавливает дробно-линейное конформное отображение [8]
3. Параметризация решения
(С + О2'
ъ
которое позволяет найти в замкнутом внде и производную отображения. Поясним, чем объясняется такой выбор вспомогательных плоскостей. В плоскости Ш неизвестен участок границы СА, и для определения функции Ш (£) этому участку границы в плоскости £ удобно сопоставить дугу единичной окружности. На двух остальных участках границы выполняются «хорошие» граничные условия
АВ : ф = 0; ВС : у = 0. (И)
Выполнение таких же граничных условий естественно потребовать для обеих аддитивных частей функции Ш(£): функций Шо(£) и Шх(^). Поэтому гладкую в П ^ функцию ) можно представить в виде степенного ряда по нечетным
степеням £
к
Ш1(С ) = ^С2й+1с2к+1, (12)
к=0
коэффициенты которого С2к+1 вещественны в силу граничных условий (11). Со-
СА
Фурье
кк У1И = ^ С2к+1 вш(2к + 1)а, ф1 (а) = ^ С2к+1 еов(2к + 1)а. (13)
к=0 к=0
Здесь символом обозначаются граничные значения соответствующих функций переменной £ на участке СА, где С = ега ■
Аналогичным образом, в плоскости г неизвестен участок границы ВС, и для определения функции г(п) ему в плоскостп п удобно сопоставить дугу единичной окружности. На двух остальных участках границы выполняются «хорошие» граничные условия
С А : у = 0; АВ : у = -1. (14)
Выполнение таких же граничных условий (с точностью до значения константы в правой части условия на АВ) естественно потребовать у обеих аддитивных частей функции г(п): функций г0(п) и г1 (п) • Поэтому функцию г1(п) можно представить в внде
к
= С-1 (1шУ - + (15)
к=0
Здесь коэффициенты ряда С2к, наряду с множителем С_1, вещественны в силу граничных условий (14). Под логарифмом понимается главная ветвь. Первое сла-
у1
А АВ АС
ВС
в ряд Фурье
к к
¿1(7) = ]ГС2йсо8(2А;7), у1(7)=С1-1(7-^)+^С12^т(2А7). (16) к=0 к=1
Здесь символом " Л" обозначаются граничные значения соответствующих функций переменной п на участке ВС, где п = еГу •
В результате такой параметризации определение функций Ш(£) и г(п) сводится к определению дискретного набора коэффициентов Ск, к = -1, 0,1,..., К, при условии, что функции Ш0(£) и г0(п) будут известны.
4. Вывод граничных уравнений и организация первого шага итераций
Условимся и далее для любых функций использовать символы "А" над знаком функции для обозначения граничных значений этих функций на участках границы CA и DC соответственно.
Граничное уравнение для функции W(Ç) на участке границы CA можно получить из граничного условия (5). Действительно, подставляя следующее граничное соотношение
СА- = — = /(—
\dx J y dX \da J / \da
CA
W(Ç):
dx
CA: *L
da
с условиями на концах интервала a G /2]
1 + £ — £-ф
А (17)
da
C : W(a) = iu; A : W(a) = —то + i. (18)
=0
n/2
При этом надо учитывать, что \¥(а) = ф(а) + гф(а) не просто функция вещественного переменного а, а обозначение относящихся к участку С А граничных значений функции W(^) комплексного переменного С € О Поэтому к задаче (17), (18) надо также добавить граничные условия для функции W(() на остальных участках границы области О ^ :
Ш )\лп =0, Ж)1 по =0 (19)
и оценку локального поведения функции W(С) в окрестности точки А:
С^г: IV (С) ~ -1п (С - «) • (20)
п
Отметим, что локальное поведение (20) функции W(О согласуется с граничными условиями (18), (19).
В свою очередь, граничное уравнение для функции г(п) на участке границы ВС очевидным образом можно получить из второго граничного условия в (6):
ВС : у = п-1ф - I. (21)
СА
ные г. Чтобы избежать при расчетах процедуры численного дифференцирования г как источника высокочастотной погрешности, граничное уравнение (21) на участке ВС
ВС п
более подходящий вид граничного уравнения для функции г(п) на участке грани-ВС
ВС: # = (22)
и условия на концах интервала 7 € [0,п/2]
в : у(1)\1=0 = -1; С : ¿(7) |7=П/2 =0. (23)
Здесь также необходимо учитывать, что ¿(7) = х(у) + гу(7) не просто функция вещественного переменного 7, а обозначение относящихся к участку ВС граничных
значений функции г(п) комплексного переменного п € П Поэтому к задаче (22), (23) надо добавить граничные условия для функции г(п) на остальных участках границы области П п :
у(п)1сл = 0; у(п)и = -1 (24)
А
?7 ~ 0 : 1п?;. (25)
п
Отметим, что локальное поведение (25) функции г(п) согласуется с граничными условиями (23), (24).
Таким образом, краевая задача (1) (4) сведена к совокупности двух граничных уравнений (17), (22) относительно функций Ш(£) и г(п), которые надо решить совместно с граничными условиями (18) (20), (23) (25).
В соответствии с принятой в п. 2 схемой решения перепишем граничные уравнения (17), (22) в виде итерационных уравнений для определения гладких добавок Ш0 (С), г1^)(п) по известному с предыдущей итерации в иду функций Ш0° 1)(0>
гГ^п), Шр'-^С), ^'-^(п):
¿а ¿а
е-1 + 1 - ф0Г-1) - ф0'-1)
¿х 0 ¿х 1
¿а ¿а
,
¿7 (¿7 ¿7 ^ ¿7 у
Ввиду принятого для ШХС) представления (13) искомая функция в левой части уравнения (26) иредставима рядом Фурье
л-О')/ \ к
_ \ - Г /о/. ,
¿а
[- (2к + 1)С(2°к)+^ ес8(2к + 1)а. (28)
к=0
Для текущей j-й итерации все функции в правой части уравнения (26) известны. Применяя к ней БПФ по почетным косинусам, можно найти коэффициенты Фурье-разложения левой его части (28), а значит, коэффициенты С^)+1 разложения (13), функцию Ш1(О) (С) и ее производную ¿Ш1О)у/'всюду в области П^.
Далее, используя связь (10) между £ и п> можно найти производную ¿Ш0^¿п всюду в области Пп и, в частности, ее значения на границе ВС.
Аналогично, ввиду принятого для ^(п) представления (16) искомая функция в левой части уравнения (27) продставима рядом Фурье
^ = + (29)
' к=1
Для текущей j-й итерации все функции в правой части уравнения (27) известны. Применяя к ной БПФ по четным косинусам, можно найти коэффициенты Фурье-разложения левой его части (29), а значит, коэффициенты С_1, С0 разложения (16), функцию г0)(п) и ее производную ¿гО'У ¿п всюду в области Пп.
Далее, используя связь (10) между £ и п> можно найти производную всюду в области П^ и, в частности, ее значения на границе СА.
Единственным препятствием к применению БПФ является наличие в правых частях уравнений (26). (27) неопределенностей типа «бесконечность минус бесконечность» на концах интервалов а € [0, п/2], 7 € [0, п/2], отвечающих точке С: а = ^ 7 = п/2. Эти бесконечности порождаются дробно-степенными особенностями главных частей решения. В уравнениях (26). (27) они с необходимостью взаимно погашаются, но при этом дают неопределенную константу. Заметим, что логарифмические особенности, присущие этим функциям в точке А: а = п/2, формально также порождают неопределенность типа «бесконечность минус бесконечность» в правой части уравнения (26). Однако при их взаимном погашении никакой ноопро-
А
а = п/2
Что касается точки В: 7 = 0, то она просто является точкой регулярности для всех функций, фигурирующих в правой части уравнения (27).
С
лонность предлагается использовать для подавления старшей гармоники, а именно: для задачи Фурье-анализа функции ду1/да потребовать выполнения условия
С2К+1 = 0, (30)
а для задачи Фурье-анализа функции ду1 /д^ условия
с 2К = 0. (31)
Каждое из этих условий линейно по отношению к искомым функциям дуч/да и ду1/д7. Поэтому условие (30) эквивалентно заданию определенного значения производной ду1/да в точке С, а условие (31), в свою очередь, эквивалентно заданию определенного значения производной ду1/д^ в точке С. Соответственно, обе задачи Фурье-анализа замкнуты.
Таким образом, первый шаг з'-й птераций завершен. На выходе получим значения коэффициентов | С ^, к = -1,0,1,..., к\ для второго шага 3 -й итераций.
5. Построение главных частей решения и организация второго шага итераций
Для главных частей решения ¿0(п) ^ Wо(Z) в качестве начальных приближений для итерационного процесса выберем некоторые базовые функции г00(п), W00(Z)
z00(ri) = z{0j)(ri) . , Wco(C) = (Z)
U),
=0
j=o
так, чтобы в плоскостях z00, Woo им отвечай области, аналогичные Q,z, Q w, отличающиеся от них только определенным фиксированным положением точки D в плоскости zoo и, соответственно, точки C в плоскости Woo • Имеющийся произвол естественно использовать с целыо добиться наиболее простого замкнутого вида функций zoo(n), Woo(Z)•
В качестве области Qz00 возьмем часть потенциального течения от цепочки источников, расположенных в точках
zoo = 2in; n = 0, ±1, ±2, ...
к стоку на бесконечности (см. затененную область на рис. 4). Мощность источников
2п
C
A
A
'■ W z
^ (0)
D ( ^d - i)
Рис. 4. Вид вспомогательной плоскости z00
_D (id) C_D (0)
W ,
C A (ip2) C (-Щ A
Рис. 5. Вид вспомогательных плоскостей: о) переменного i, б) переменного 0
Комплексный потенциал £ такого течения представляет собой полуполосу, изображенную на рис. 5, а, где d = вс + п/2. Соответственно, конформное отображение О200 ^ О I реализует функция [9]
. . / nzoo\ . . 4г ,
t {z00j = ln (^tg -j— J z0o(t) = — arctg (
(32)
Отсюда, в частности, следует, что
4
хоо,п = (гег,3с)} .
В свою очередь, области О г и О п канонического вида связаны между собой конформным отображением
п + г
¿d.
(33)
В результате комбинация отображений (32), (33) дает вид функции гоо(п)- При этом в многозначных функциях 1п и ак^ необходимо выделять однозначные ветви в соответствии с видом областей О200, О 4, О п.
Аналогично, в качестве области О^00 возьмем часть потенциального течения от цепочки источников, расположенных в точках
- 1 • « I--Фоо,с + —) ; п = 0, ±1, ±2, ...
Woo = i ( Фоо,с + J ; г ур^
к стоку на бесконечности (см. затененную область на рис. 6). Мощность источников несущественна и для определенности также выбирается равной 2п.
Комплексный потенциал в такого течения представляет собой полосу, изображенную на рис. 5, б. Соответственно, конформное отображение О^00 ^ О д реализует функция [9]
в (Woo) = ln [cos (¿всWoo) - cos (вс^oo,c)] - ln [1 - cos (в^оо,с)],
t
e
А
А
кРРс)
С (Уоо,С )
А
А
чЬ
о (0)
Рис. 6. Вид вспомогательной плоскости Woo
а обратное отображение функция
Т'Уоо {°) = агесов {е9 [1 - соэ (/Зсфоо,с)] + соэ (/Зсфоо,с)} • (34)
В свою очередь, области О д и О ^ канонического вида связаны между собой конформным отображением
В результате комбинация отображений (34), (35) дает вид функции Шею (О-При этом в многозначных функциях 1п и агееов необходимо выделять однозначные ветви в соответствии с видом областей О^00, О д, О
Отметим, что вид области О^00 и, соответственно, функции Ш00(С) зависит от параметра ф00,с. Наиболее простой вид функции получается в случае, когда Фоо,о = п/(2вс)■ Тем не менее, целесообразно оставить возможность изменения параметра ф00,с. Это объясняется тем, что при фиксированном значении ф00,с в области малых £ и 6, а именно при £ < 0.^ щи 6 < 0.3, возникают проблемы с достижением необходимой точности численного решения граничного уравнения для функции Ш(С). Подбирая определенным образом величину параметра ф00,с, удается достичь удовлетворительной точности решения для более широкой области £6
Организацию второго шага итераций естественно связать с удовлетворением тех граничных условий на полные функции г(п), Ш(С), которые не были удовлетворены за счет первого шага. В частности, после завершения первого шага итерации будет получено определенное значение С—1. В соответствии с представлением (15) функции г^ (п) это означает, что па гр аницах С А АВ будут выполняться усло-
Чтобы функция г(п), составленная го двух аддитивных частей г0(п) и на этих же границах удовлетворяла условиям (14), на ее главную часть необходимо наложить условия
(35)
ВИЯ
СМ: ^=0, ВС: у™ =
С А: уЫ= 0, ВС: у^ = - 1,
(36)
при этом базовая функция г00(п) (см. рис. 4) удовлетворяет условиям
СА : У00 =0, ВС : У00 = -I-
(37)
В соответствии с соотношениями (36), (37) для главной части 20 (п) функции г(п) можно организовать второй шаг j-YL итерации как конформное отображение
(7)
области П200 та область П^0 • Последняя то виду аналогична П200 и отличается только величиной скачка в точке А при переходе с границы АД на СА. Поскольку величина этого скачка единственный характерный размер области, организовать отображение П200 ^ П7) можно как простое умножение базов ой функции £00(п) на соответствующий множитель
(„)= 1 - гоо(п)- (38)
Аналогичным образом, после завершения первого шага j-й итерации будут получены определенные значения коэффициентов , к = 0,1,...,К.В соответствии с представлением (12) функции Ш7 (£) это означает, что в точке С она принимает определенное значение
С : Щ7 = , ^ - Е7
к=0
Чтобы функция Ш(£), составленная из двух аддитивных частей Шо(С) и (С), удовлетворяла уел овию Шс = ¿и (см. условие (5) и рис. 2), на ее главную часть необходимо наложить условия
С : Ш07) = ¿7, ^С - и - ЕС^^, (39)
к=0
при этом базовая функция Ш00(С) (см. рис. 6) удовлетворяет условиям
С : Ш00 = г^00,с, (40)
где "000,с _ вспомогательный параметр, имеющий определенное значение (оно выбирается один раз, перед началом всех итераций).
В соответствии с соотношениями (39), (40) для главной части Ш0(£) функции Ш(С) можно организовать второй шаг j-й итерации как конформное отображение области 0^00 та область П^ • Последняя то виду аналогична 0^00 и отличается только величиной отрезка СД. Надо сказать, что организовать это отображение аналогично (38), как простое умножение базовой функции Ш00(С) на некоторый множитель, невозможно. В отличие от области П200, у области П^00 не один характерный размер, а два - это расстояние СД = ^00,с и скачок на 1 в точке А при переходе с границы АД на СА, и величина скачка должна оставаться той же самой (7)
для области П^0 • Поэтому приходится строить более сложное конформное отоб-(7)
ражение П^00 ^ П^, причем здесь возможны различные варианты. В частности, можно предложить такой вариант отображения:
ш07) (С) =
1+ А
А 00 + Ш00>(С) ]
Ш00(С), А00 > ^020,с, (41)
А00
значение которого выбирается один раз, перед началом всех итераций, аналогично параметру ^00,с.
При условии А 00 > ^00 с отображение (41) конформно всюду в области П^00. Более того, сохраняются не только все углы на границах области, но и сам вид
области - меняется только величина участка границы СБ. Остается только подобрать итерационный параметр ^^ так, чтобы отображение (41) давало нужное значение 'Ф^с (см- условие (39))
= ( Аоо ~ Фоо,с \ 'Фоо,с
K
-Е - Фоо,о
k=0
(42)
Выражения (38), (41), (42) представляют собой второй шаг итерации. Необходимым условием сходимости итераций будет, очевидно, существование пределов
lim
И}
= •
K
K
^ Еc2jk+i\= Y.c2k+i.
k=0
k=0
Соответственно, за критерий выхода из итераций принималось условие
max
1 -
c
(j-i)
K
K
-i
c
(j)
i -Е cj^E C
2k+1
< 10
-6
k=0
k=0
6. Анализ результатов расчетов и сравнение с имеющимися асимптотическими оценками
Для использования БПФ задавалось равномерное разбиение интервалов а £ £ [0, п/2] и y € [0, п/2] с число уз лов N = 2s + 1, где s выбиралось от 10 до 14. По окончанию итераций находилась численная невязка r(a) и h(Y) уравнений (17), (21) _ _ + н
' dtp \ / f dx\
CA : r(a)
da
da J
- £ 1 - ф<
1,
DC : h(Y)
1 / *
ф - У
1,
(43)
(44)
где звездочкой помечены функции, полученные в результате описанного выше итерационного процесса. Уточним, что невязки гп вычислялись в узлах сетки ап, п = 1,..., N, а невязки Нп - в узлах се тки п = 1,..., N. Как правило, все они оказывались меньше 10-4.
В качестве параметра фоо,с выбиралось значение, близкое к 0.5, скажем, в интервале [0.4,0.55], в качестве параметра А 00 - значение из интервала [2, 5]. Выбор различных значений фоо,с, А оо из указанных интервалов практически не влиял па скорость сходимости итераций и величину невязок тп, Нп.
Вместе с тем у описанного метода решения есть один неустранимый недостаток, который существенно ограничивает эффективность метода в области малых е и 6, несмотря та сходимость итераций и малость невязок тп, Нп. В результате конформного отображения П п ^ П 2 равномерное разбиение дуги единичной окружности Б С : 7 £ [0, п/2] в плоскос ти п узлами 7п приводит к неравномерному разбиению участка границы БС в физической плоскости г: в окрестности С
го отображения П ^ ^ П ^ равномерное разбиение дуги единичной окружности С А : а £ [0, п/2] в плоскос ти £ узлами ап приводит к неравномерному разбиению границы С А в плоскости Ж: сетка будет сильно разреженной в окрестности точки СА
сетки в плоскости Ж по мере того, как а приближается к п/2.
Рис. 7. Конфигурации контура ПО, полученные численно (сплошные линии) и с помощью асимптотического разложения [1](штрих-пунктирные линии), для е = 0.3 и 6 = 0.3 (о), 5 = 0.1 (б)
Рис. 8. Конфигурации контура ПО, полученные численно (сплошные линии) и с помощью асимптотического разложения [3](штрих-пунктирные линии), для е = 0.1 и 6 = 0.3 (о), 6 = 0.1 (б)
Как следует из формулы (6), чем меньше е и 5, тем угол вс меньше, и эффект разрежения сетки в окрестности точки С в плоскостях 2 и Ш будет выражен сильнее. В частности, при е = 0.1 и 5 = 0.1 и числе узлов сетки N = 210 + 1 величина интервала сетки в окрестности точки С в плоскостях 2 и Ш достигает нескольких единиц, и увеличение числа узлов, скажем, до N = 214 + 1 практически ничего не меняет. Соответственно, большие повязки между узлами сетки такого большого интервала приводят к значительной погрешности определения конфигурации межфазной границы ДС в целом, причем невязки уравнения (21) - непосредственно, а невязки уравнения (17) опосредованно.
Это диктует необходимость контролировать невязки (43), (44) не только в узлах сетки 7„ и ап, то также та последнем интервале разбиения по 7, то есть па [7«-1,7^]) и, кроме того, та первом и последнем интервалах разбиения по а, то есть на [а1?а2] и [а^-1,а«]. Обозначим эти невязки через Н(7), г 1,2(а), г N-1,« (а) соответственно.
В результате такого контроля было выявлено, что только в области, где е > 0.3 и 5 > 0.3, при любых ^оо,с> А 00 из указанных выше диапазонов невязки
Табл. 1
(5 e
0.3 0.1
0.3 0.38, 0.35 0.36, 0.36
0.1 0.49, 0.45 0.48, 0.46
hn —i,n(y) > r 1,2(a), rn—i,n(а) те превышают 0.03. С уменьшением как £, так и 6 (вместе или то отдельности), невязки h n —i,N (y) растут слабо. В то же время невязки r i,2(a), r n—i,n (а) начинают расти быстро и достигают 10 ^ 20% уже при £ = 0.1 и/ил и 6 = 0.1. Подбор определенных з начений ф00,0, А 00 из указанных выше интервалов позволяет значительно уменьшить эти невязки. В результате удается добиться невязок h n-i,n (y) , r i,2(a), r n—i,n (а) не более 0.03, вплоть до £ = 0.1 и/ил и 6 = 0.1.
Соответственно, в диапазоне 0.1 ^0.3 значений £ и 6 можно с удовлетворительной точностью рассчитать конфигурации межфазной границы DC и сравнить их с найденным в [3] главным членом асимптотического разложения решения. Отметим, что основной результат [3], состоит в том, что в крупномасштабном приближении контур DC можно считать прямой линией, наклоненной к горизонту под углом ¡5°с = arctg (1/а/и — 1) ~ е1/2 . При этом оставался открытым вопрос, какова точность такой асимптотической оценки? Ответ очевидным образом связан с выясне-
£
разложение.
DC
могцыо описанного выше итерационного метода решения задачи (1) (4) (сплошные линии) и полученные с помощью асимптотического разложения [3] решения той же задачи (штрих-пунктирные линии). Рис. 7 отвечает £ = 0.3 и 6 = 0.3 (о), 6 = 0.1 (б). Рис. 8 отвечает £ = 0.1 и 6 = 0.3 (а), 6 = 0.1 (б). Конфигурации совмещались D
наименьших ошибок разложения.
В табл. 1 для различных £ и 6 в диапазоне 0.1 ^ 0.3 приведены относительные ошибки асимптотической оценки величины xo - xd (первое число) и величины радиуса кривизны R d свободной границы в точке D (второе число). Как видно, эти ошибки практически совпадают и по порядку величины составляют O(£i/2) . В результате можно сделать вывод, что асимптотическое разложение [3] организовано по степеням £i/2 и относительная ошибка определения конфигурации контура DC по формулам главного члена этого разложения составляет O(£i/2).
Работа выполнена при поддержке Российского фонда фундаментальных исследований (проект 08-01-00548).
Summary
М.М. Alimov. Iterative Solution of the Problem of Liquid Impregnation into Laminated Porous Material.
The problem of liquid impregnation into laminated porous material is considered. For case with high contrast of layers' thickness and permeability, this problem is formulated as specific free boundary problem. Numerical solution is obtained by the iterative method similarly to the Levi Civit.a method in liydrodynamic of an ideal liquid. A comparison of asymptotic analysis results and numerical results is presented.
Key words: free boundary problems, multiphase mediums, liquid impregnation into porous material, VARTM technology.
Литература
1. Hsiao К, Т., Mathur R., Gillespie J. W. Jr., Fink B.K., Advani S.G. A closed form solution for flow during the vacuum assisted resin transfer molding process // J. Manuf. Sci. Eng. 2000. V. 122, No 3. P. 463 475.
2. Mathur R., Heider D., Hoffmann C., Gillespie J.W. Jr., Advani S.G., Fink B.K. Flow front measurements and model validation in the vacuum assisted resin transfer molding process // Polym. Compos. 2001. V. 22, No 4. P. 477 490.
3. Alimov M.M., Kornev K.G. Impregnation of liquids into a laminated porous material with a high permeability contrast // Pliys. Fluids. 2007. V. 19, No 10. P. 102108-1 102108-11.
4. Muskat M. The Flow of Homogeneous Fluids through Porous Media. Ann Arbor, Michigan: I.W. Edwards, Inc., 1946. = Маскет M. Течение однородных жидкостей в пористой среде. М., Л.: Гостоптехиздат, 1949. 628 с.
5. Полубарииова-Кочииа П.Я. Теория движения грунтовых вод. М.: Наука, 1977. 664 с.
6. Saffman P.G., Taylor G.I. The penetration of a fluid into a porous medium or Hele-Sliaw cell containing a more viscous liquid // Proc. Roy. Soc. London. Ser. A. 1958. V. 245, No 1242. P. 312 329.
7. Гуревич М.И. Теория струй идеальной жидкости. М.: Наука, 1979. 536 с.
8. Лаврентьев М.А., Шабат Б.В. Методы теории функции комплексного переменного. М.: Наука, 1973. 736 с.
9. Koppenfels W., Stalman F. Praxis der Konformen Abbildung. Berlin: Springer, 1959. = Коппеифельс В., Штальмаи Ц. Практика конформного отображения. М.: ИЛ, 1963. 406 с.
Поступила в редакцию 26.02.08
Алимов Марс Мясумович кандидат физико-математических паук, ведущий научный сотрудник НИИ математики и механики им. Н.Г. Чеботарева Казанского государственного университета.
Е-шаП: Mars.AlimovQksu.ru
