Диссертация (1137359), страница 9

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

Текст из файла (страница 9)

Область, в которой ищется решение задачи, не ограничена, т.е.θ0 ∈ [0, ∞).2. Искомая функция u∗1 не является ограниченной при θ0 → ∞.Для решения первой проблемы можно сделать замену координат σ = θ0 /(1 + θ0 ), переводящую полуось [0, ∞) в луч [0, 1). Однако,эту замену можно не делать, т.к. достаточно хорошая точность вычислений может быть получена простым путем ограничения области в которой ищется решение, например можно взять θ0 ∈ [0, 100],см. [73].Для решения второй проблемы нужно сделать заменуdefH =u∗1f 00 (0)− (θ + µ) √ .x0(1.113)Тогда уравнение (1.112) примет следующий вид:∂H ∂Hf 00 (0) 0f 00 (0) ∂µ+H + √ (θ + µ) + √H − limH−θ0 →∞∂t∂ξxx ∂ξ ∂Hf 00 (0)f 00 (0) ∂µ 0∂ 2Hθ =, (1.114)− √D − limD − 0 D+ √θ0 →∞∂θ∂θ02xx ∂ξгдеdefZθ0D =∂H 00dθ .∂ξ(1.115)0Применим технику осреднения.

Считая граничное условие H θ0 →∞e так, чтобы lim He = 0.неизвестным, мы выберем функцию H0θ →∞Тогда осредняя уравнение (1.114), учитывая что D ≡ 0, мыполучим уравнение на H:∂H∂ 2H−W =∂t∂θ02(1.116)73с граничными условиямиHθ0 =0= 0,∂H → 0,∂θ0 θ0 →∞(1.117)гдеZ2π e 0000 (0) ∂µe∂Hf1f(0)∂µ∂Hdef0W =D+ √θ0 =D+ √θ dξ.∂θ02π∂θ0x ∂ξx ∂ξ0(1.118)e будет иметь видТогда уравнение на осциллирующую часть H00ee ∂Hf 00 (0) ∂µ ef(0)∂H0e√H +H ++(θ + µ) + √H + H − limH−θ0 →∞∂t∂ξxx ∂ξe ∂Hf 00 (0) ∂µ 0∂ 2Hf 00 (0)√D − limD−D+θ=, (1.119)− √θ0 →∞∂θ0∂θ02xx ∂ξe будут иметь вид:а граничные условия на функцию Hf 00 (0)eH θ0 =0 = −µ √ ,xe 0H→ 0,θ →∞e = HeH.ξξ+2π(1.120)Начальное условие мы выбрали в следующем виде: 00− f√(0) µ(ξ)(1 − 0.2θ0 ), θ0 6 5,e = = 0. (1.121)xHHt=0t=00, θ0 > 5,Отметим, что для задачи (1.116)–(1.121) в [50] была построенанеявная схема, но на практике удобнее использовать явные схемы(см.

[42]), ввиду легкости ее распараллеливания (см. [14]) для вычислений на современных многопроцессорных ЭВМ (а также длявычислений с использованием графических ускорителей, например,технологии CUDA, см. [29; 30]). Построим для задачи (1.116)–(1.121)явную разностную схему.74Введем сетку ω с шагами hξ , hθ0 , ht на области0Ω = [0, 2π] × [0, θmax] × [0, Tmax ],в которой мы будем искать решение:nω = (i, j, k) : ξi = ihξ , θ0 j = jhθ0 , tk = kht ,i = 0, . . . , N ;j = 0, . . .

, M ;N = 2π/hξ ,M=k = 0, . . . , K;0θmax/hθ0 ,oK = Tmax /ht .ek = He tk , ξi , θ0 j , H kj = H tk , θ0 j ,Введем сеточные функции: Hi,jkk 0kk0Di,j = D t , ξi , θ j , Wj = W t , θ j , uki,j = u∗ tk , ξi , θ0 j , µi = µ(ξi ),kvi,j= v ∗ tk , ξi , θ0 j , где функции D и W были определены ранее,см. (1.115) и (1.118) соответственно.Введем следующие обозначения:defµ̂i =µi+1 − µi−1,2hξk defβi,j=kDi,j00def f (0)F = √ ,x+F θj0 µ̂i ,kk def e kαi,j= Hi,j + H j + F (θj0 + µi ),defχki,j =kβi,jek − HekHi,ji,j+12k:и для любой сеточной функции fi,jkfi,j+ kkf +fdef i,ji,j=,2kfi,j− kkf f−def i,ji,j=,2kгде интеграл Di,jвычисляется методом трапеций (см. [43]):kDi,j= hθ 0j−1dki,0 + dki,j X k+di,j 0 ,20j =1гдеdki,jek − HekHi+1,ji,j=hξ75Тогда явная разностная схема для задачи (1.119), (1.120), (1.121)будет иметь видek − Hekek − Hek HHi,ji−1,ji+1,ji,jk −+ αi,j+hξhξkk kkke i,j+ F µ̂i H+ H j − H M − F Di,j− Di,M−kkkk H j − H j−1 H j+1 − H jk +k −− βi,j+ βi,j+hθ 0hθ0e k − 2Hek + HekHi,j+1i,ji,j−1,+ ht2hθ0ke k+1 = He i,jH− hti,jk +αi,jkkkke i,0e i,Me 0,je N,jH= −µi F, H= 0, H=H,−F µ (1 − 0.2θ0 ), θ0 6 5,ijj0eHi,j =0, θ0 j > 5.для задачи (1.116), (1.117), (1.121) будет иметь видkk+1Hj=kkHjkkH j+1 − H j + H j−1+ ht Wjk + ht,h2θ0H 0 = 0,kkH M = HM−1 ,0H j = 0,где интеграл Wjk вычисляется методом трапеций (см.

[43]): kN−1kXχ+χh0,jξN,jWjk =+χki,j .2π2i=1Легко показать, что построенная разностная схема удовлетворяет принципу максимума (см. [42]), и, следовательно, устойчива.Искомая функция uki,j находится по следующей формуле:kke i,juki,j = H+ H j + F (θ0 j + µi ),76kа функция vi,jнаходится из интеграла (2.52) с помощью численногоинтегрирования методом трапеций:vi,j = hθ0M−1kkXζi,0+ ζi,Mk+ζi,j ,2j=1гдеkζi,juki,j+1 − uki,juki+1,j − uki,j+ µ̂i.=−hξhθ 0Отметим, что результаты моделирования, полученные по построенной явной разностной схеме, совпадают с результатами, приведенными в работе [50].

В силу этого, не претендуя на оригинальность, приведем их здесь. Заметим, что построенная разностная схема будет использоваться при моделировании задачи о течении в трубе с малыми периодическими неровностями на стенке, которая будетисследоваться в Главе 3.Также введем функцию eps(t) — критерий стационарности:kk+1keps(tk+1 , tk ) = max kuk+1i,j − ui,j k, kvi,j − vi,j k ,kk| Соответственно, если с некоторого моk = max |gi,jгде норма kgi,j(i,j)∈ω∗мента времени t eps(t) = 0 ∀ t > t∗ , то наше решение вышло настационар.Будем считать, что неровность задается функцией видаµ(ξ) = A cos(ξ),A = const,и зафиксируем x = 1 (x — параметр — расстояние от кромки пластины). Из результатов [50] видно, что зависимость течения от амплитуды неровностей имеет несколько характерных типов.

При A 6 1течение ламинарное (см. рис. 1.7, 1.8), а при A = 1.8 в течении возникает периодический вихревой процесс — на левой стенке «ямки»возникает вихрь (см. рис. 1.10–1.13), который двигаясь по потокуразрушается у правой стенки «ямки», и затем на левой стенке воз-77никает новый вихрь и т.д. Из этого следует, что существует некоторая критическая амплитуда A∗ , такая, что при A < A∗ течениеламинарное, а при A > A∗ течение становится вихревым.Полученные нами результаты подтверждают результаты из [50]и уточняют их: существуют две критические амплитуды A∗1 и A∗2 >A∗1 , такие, что при A < A∗1 течение ламинарное и с некоторого момента времени — стационарное (например, при A = 0.8 см.

рис. 1.7–1.9), при A∗1 < A < A∗2 в течении образуется вихри, но начинаяс некоторого момента времени течение становится ламинарным истационарным (например, при A = 1.8 см. рис. 1.10–1.14), а приA > A∗2 в потоке также наблюдается периодический процесс образования вихрей, но начиная с некоторого (большого) времени течениестановится стационарным, но не ламинарным, т.е.

в «ямке» возникает стационарный вихрь (например, при A = 2.3 см. рис. 1.15).Отметим, что использование явной схемы вкупе с использованием многопроцессорных ЭВМ позволяет моделировать задачу напромежутках времени, существенно больших, чем при использовании неявной схемы, даже не смотря на необходимость использованияочень маленьких шагов по времени.Заметим, что на приведенных рисунках 1.10–1.15 отчетливовидно явление отрыва пограничного слоя и его повторного прилипания к обтекаемой поверхности.78Рисунок 1.7 – Поле скоростей при µ = 0.8 cos ξ, t = 0Рисунок 1.8 – Поле скоростей при µ = 0.8 cos ξ, t = 10Рисунок 1.9 – График функции eps(t) при µ = 0.8 cos ξ79Рисунок 1.10 – Поле скоростей при µ = 1.8 cos ξ, t = 2Рисунок 1.11 – Поле скоростей при µ = 1.8 cos ξ, t = 380Рисунок 1.12 – Поле скоростей при µ = 1.8 cos ξ, t = 5Рисунок 1.13 – Поле скоростей при µ = 1.8 cos ξ, t = 9Рисунок 1.14 – График функции eps(t) при µ = 1.8 cos ξ81а) t = 2.25б) t = 6.7Рисунок 1.15 – Поле скоростей при µ = 2.3 cos ξ824.2.Алгоритм численного решения уравнения типаРэлеяВ этом параграфе мы приводим алгоритм численного решениязадачи (1.13), (1.14), для удобства перепишем ее здесь, подставивявный вид u†0 (см.

(1.16)):√ Zξ000√fτ/x∂= 0,ε1/3 ∆ ve2II dξ + f 0 τ / x ∆ev2II − ve2II∂txve2II τ =0 = lim ve2∗ , lim ve2II = 0, ve2II ξ = ve2II ξ+2π .τ →∞θ→∞(1.122)(1.123)В качестве начального условия возьмем2ve2II t=0 = v∞ e−τ /10 ,(1.124)гдеdefv∞ (t, x, ξ) = limve2∗ .0(1.125)θ →∞Отметим, что интегралмулеZξZξfedξ =Rξfedξ можно вычислить по следующей фор-1fe(ξ ) dξ −2π000Z2π Zξ 000fe(ξ ) dξ000dξ 0 .(1.126)0Сделаем замены:defg =Zξve2II dξ,defG = ∆g.(1.127)С учетом замен (1.127) уравнение (1.122) примет вид√ 000√fτ/x∂G∂G+ ε−1/3 f 0 τ / x= ε−1/3 ve2II.∂t∂ξx(1.128)83Тогда алгоритм нахождения ve2II следующий.1.

Найти функцию G из уравнения (1.128) c условиями периодичности: Gξ = Gξ+2π .2. Найти функцию g из уравнения∆g = G.(1.129)3. Найти функцию ve2II по формулеve2II =∂g.∂ξ(1.130)Как и в предыдущем параграфе, ограничим область, в которойищется решение до τmax = 100. Введем сетку ω2 с шагами hξ , hτ ,ht на области Ω2 = [0, 2π] × [0, τmax ] × [0, Tmax ], в которой мы будемискать решение:nω2 = (i, j, k) : ξi = ihξ , τj = jhτ , tk = kht ,i = 0, . . . , N ;j = 0, . . . , M ;N = 2π/hξ ,k = 0, .

. . , K;M = τmax /hτ ,oK = Tmax /ht .kВведем сеточные функции: vi,j= ve2II tk , ξi , τj , Gki,j = G tk , ξi , τj .Сначала найдем численно (например, с помощью метода Рунге–Кутта, см. [43]) функцию Блазиуса (см. (4) на стр. 13):2f 000 + f · f 00 = 0,f (0) = f 0 (0) = 0,f (∞) = 1.Введем обозначения:(1)fj√ def 0= f τj / x ,(3)fj√ 000fxτ/jdef=.xТеперь перейдем к приведенному алгоритму. Уравнение (1.128)будем решать численно с помощью разностной схемы (см. [42]), ко-84торая довольно тривиальна:Gk+1i,j=Gki,j−ε−1/3k(1) Gi,jht f j− Gki−1,j(3) k+ ε−1/3 ht fj vi,j.hξПункты 2 и 3 приведенного выше алгоритма удобно решать,разложив функции в ряды Фурье по переменной ξ:g=Xk6=0ĝk (t, τ )eikξ ,G=XĜk (t, τ )eikξ ,ve2II =Xv̂k (t, τ )eikξk6=0k6=0Тогда уравнение (1.129) примет следующий вид:∂ 2 ĝk− k 2 ĝk = Ĝk ,2∂τk 6= 0,а из граничных условий (1.123) легко получитьv̂∞,kĝk τ =0 =,ikĝk τ →∞ → 0,где v̂∞,k — коэффициенты разложения в ряд Фурье функции v∞ ,k 6= 0.

Эту задачу мы будем решать с помощью неявной разностнойсхемы методом прогонки (см. [42]).Коэффициенты v̂k будем находить по формуле (1.130), котораяпримет следующий вид:v̂k = ikĝk .Мы не будем приводить численную реализацию данного алгоритмав силу его тривиальности.Перейдем к результатам численного решения. В качестве функции, описывающую неровность мы взялиµ(ξ) = 2.3 cos(ξ).Как было показано в предыдущем параграфе, при таком выборефункции в тонком слое возникают вихри, см. рис. 1.15.85Заметим, что в переменных (x, y) при малых ε классическийпограничный слой Прандтля сосредоточен в малой окрестности осиx, см. рис.

1.16. Поэтому мы приводим картину течения в переменных (τ, ξ), см. рис. 1.17, и, для наглядности, в переменных (τ, x),см. рис. 1.18. В последнем случае мы растягиваем период по переменной ξ в окрестности выбранных точек на оси x (напомним, чтов «главном» скорости потока зависят от x как от параметра).На рис. 1.18, 1.17 мы изображаем главные слагаемые построенного нами решения, которое имеет вид√ 2/3eII,u(τ, ξ) = f 0 τ / x + ε1/3 u1 +O εv(τ, ξ) = ε2/3 ve2II + O(ε),где функция ueII1 определяется формулой (см. (1.21))√ 00τ/fx√ueII1 =µxНа рис.

Характеристики

Список файлов диссертации