Диссертация (1145332), страница 20
Текст из файла (страница 20)
К уравнению (4.41)присоединим следующие граничные условия:Θ(0) = Θ1 ,Θ(τ ) = Θ2 .(4.42)Для нахождения решения системы (4.39)-(4.42) применим метод простойитерации с параметром. В соответствии с которым, начальное приближениедля температуры Θ(z) (например, линейное приближение, которое соответствует нулевому значению правой части уравнения (4.41)) обозначим черезΘh0i (z). Затем, подставляя Θh0i (z) в (4.39) вместо функции Θ(z), находимрешение задачи (4.39)-(4.40) и обозначим его через f h1i (z, ν).
Далее находимрешение задачи (4.41)-(4.42), соответствующее заданной функции f h1i (z, ν), иe h1i (z). Выберем малый положительный параметр α и пообозначим его как Θe h1i (z)+(1−α)Θh0i (z) в качестве следующего приближенияложим Θh1i (z) = αΘфункции Θ(z). Затем, подставляя Θh1i (z) вместо функции Θ(z) в уравнение(4.39), находим следующее приближение f h2i (z, ν) и т.д. Таким образом, наe hji (z) для нахождения слеj-ом шаге мы используем функции Θhj−1i (z) и Θдующего приближения функции Θ(z) по формуле:e hji (z) + (1 − α)Θhj−1i (z).Θhji (z) = αΘ(4.43)Основная сложность в численной реализации итерационного алгоритма заключается в решении краевой задачи для уравнения переноса теплового излучения. Для нахождения решения задачи (4.39), (4.40) будем использоватьрекурсивный алгоритм, основанный на методе Монте-Карло.Предварительно сформулируем некоторые ограничения. Пусть функцияΘ(z) является неотрицательной и Θ(z) ∈ Cb (0, τ ), где Cb (X) есть Банахово пространство функций непрерывных и ограниченных на X с нормой142kϕkCb (X) = sup |ϕ(x)|.
Также, пусть p(ν, ν 0 ) ∈ Cb ([−1, 0) ∪ (0, 1]) иx∈X121Zp(ν, ν 0 )dν 0 = 1.−1Заметим, что оператор B : Cb (Γ+ ) → Cb (Γ− ) является линейным, ограниченным и неотрицательным, kBk < 1.Для дальнейших рассуждений обозначимX = (0, τ ) × [−1, 0) ∪ (0, 1] ,(0, ν ∈ (0, 1],ξ(ν) =τ , ν ∈ [−1, 0).Введем также следующие интегральные выражения:1(Aϕ)(z, ν) =νZzz − z0ϕ(z 0 , ν)dz 0 ,exp −νξ(ν)λ(Sϕ)(z, ν) =2Z1p(ν, ν 0 )ϕ(z, ν 0 )dν 0 ,−1z − ξ(ν)(T f )(z, ν) = (Bf )(ξ(ν), ν) exp −+ (ASf )(z, ν).ν(4.44)Отметим, что согласно [152] при выполнении неравенств: kBk < 1 и λ < 1,существует единственное решение задачи (4.39)-(4.40), которое может бытьпредставлено в виде следующего ряда Неймана, сходящегося в норме Cb (X):f (z, ν) =∞X(T k f0 )(z, ν),(4.45)k=0z − ξ(ν)f0 (z, ν) = exp −h(ξ(ν)) + (1 − λ)(AΘ4 )(z, ν).νПомимо задачи (4.39)-(4.42) определения температурного профиля, рассмотрим также задачу улучшения теплоотдачи от границ слоя.
Для ее фор-143мулировки введем следующую величину, характеризующую отдачу тепла отнагретых стенок слоя:ZτΘ(z)dz.J=(4.46)0Очевидно, что большее значение J соответствует лучшей отдачи тепла отстенок слоя, то есть лучшему их охлаждению.
Величина J зависит от коэффициентов задачи (4.39)-(4.42), в частности тех, которые характеризуютматериал внешней части слоя: коэффициентов зеркального и диффузногоотражения. Пусть ρs1 = ρs2 = ρs и ρd1 = ρd2 = ρd . Сформулируем следующуюзадачу определения оптимальных отражающих характеристик границ слоя,обеспечивающих лучшую теплоотдачу:Требуется определить значения коэффициентов ρs , ρd , 0 ≤ ρs + ρd < 1,обеспечивающих максимальное значение функции (4.46), то естьJ(ρs , ρd ) → max,8.2.0 ≤ ρs + ρd < 1.(4.47)Рекурсивный алгоритм решения краевой задачи для моделирадиационно-кондуктивного теплообменаРассмотрим итерационный алгоритм (4.43). Ограничимся случаем изотропного рассеяния. Для вычисления решения задачи (4.39)-(4.40), соответствующего заданной функции Θ(z), предложим рекурсивные алгоритмы, основанные на методе Монте-Карло.
В соответствии с [152] существует единственное решение задачи (4.39)-(4.40), которое представляется в виде рядаНеймана (4.45).Применим алгоритм метода Монте-Карло для вычисления конечной суммыfN (z, ν) =NX(T n f0 )(z, ν).(4.48)n=0Для реализации метода перепишем (4.48) в виде следующего рекуррентного144соотношения:fn (z, ν) = (T fn−1 )(z, ν) + f0 (z, ν),n = 1, 2, ..., N.Рассмотрим структуру оператора T (см. равенство (4.44)). Он состоит издвух слагаемых: первое описывает эффекты отражения, второе – эффектырассеяния.
Используя простые преобразования, запишем второе слагаемое вправой части (4.44) в следующей форме:z−ξλ1 − exp −×I(z, ν) =2νZ zZ 1exp (−(z − z 0 )/ν)p(ν, ν 0 )f (z 0 , ν 0 )dν 0 dz 0 , (4.49)−1 ν(1 − exp (−(z − ξ)/ν))ξгде ξ = ξ(ν). В соответствии с методом Монте-Карло мы аппроксимируеминтеграл в этом выражении как математическое ожидание от случайных последовательностей, которые определяются через случайные величины z 0 и ν 0 ,распределенные на интервалах (ξ, z) и (−1, 1) с плотностямиexp(−(z − z 0 )/ν)иν(1 − exp(−(z − ξ)/ν))1p(ν, ν 0 ).2(4.50)Следовательно, интеграл (4.49) аппроксимируется следующей суммойλI(z, ν) =M XMz − ξ(ν)1 − exp −f (zk , νk ).νk=1Здесь zk , νk , k = 1, 2, ..., M , есть независимые реализации случайных величинz 0 и ν 0 , распределенных на интервалах (ξ, z) и (−1, 1) с плотностями (4.50).Следовательно, мы можем следующим образом аппроксимировать функцииfn (z, ν), n = 1, 2, ..., N :M1 Xfn (z, ν) ≈ f n (z, ν) =sn,k (z, ν),Mk=1f 0 (z, ν) = f0 (z, ν),(4.51)145z − ξ(ν)sn,k (z, ν) = Bf n−1 (ξ(ν), ν) exp −+νz − ξ(ν)f n−1 (zk , νk ) + f0 (z, ν).+λ 1 − exp −ν(4.52)Таким образом, конечная сумма (4.48) может быть вычислена на основе рекуррентных соотношений (4.51), (4.52).Для конструирования устойчивого варианта итерационной процедуры запишем аналитическое представление решения задачи (4.41)-(4.42) при заданной функции f .
После интегрирования равенства (4.41) получим1Θ(z) =2NcZ zZ1f (ζ, ν)νdνdζ + C1 z + C2 .0(4.53)−1Постоянные C1 и C2 определяются из граничных условий (4.42). Предполоhj−1iжим, что приближение температуры Θ(z) уже найдено на (j − 1) -омшаге итерационной процедуры (4.43). На основе (4.43) и (4.53) мы получаемследующую аппроксимацию для температуры на j-ом шаге:hjihj−1iΘ (z) = (1 − α)Θ(z) + αzM NcMX!hjiνk f N (zk , νk ) + C1 z + C2 , (4.54)k=1где случайные величины zk и νk являются равномерно распределенными наhjiсоответствующих интервалах (0, z) и (−1, 1). Вычисление f N (z, ν) осуществляем на основе формул (4.51) и (4.52), полагая, что температура равнаhj−1iΘ(z). Постоянная C1 также вычисляется с помощью метода Монте-Карло.hjiТаким образом, приближение Θ (z) функции Θhji (z) вычисляется с помощью формул (4.51), (4.52) и (4.54).Для вычисления второго слагаемого в правой части равенства (4.52) генерируются точки (zk , νk ), k = 1, ..., M , на соответствующих интервалах(ξ(ν), z) и (−1, 1) с плотностями (4.50).
Генерация точек zk осуществляетсяпо следующей формуле:zk = z + ν ln (1 − αk (1 − exp(−(z − ξ)/ν))) ,146здесь αk является независимой реализацией равномерно распределенной на(0, 1) случайной величины. Генерация угловой переменной νk определяетсяконкретным видом фазовой функции p(ν, ν 0 ).При вычислении первого слагаемого в правой части равенства (4.52) используется угловая переменная, соответствующая зеркальному отражению:−ν, и также генерируется угловая переменная, соответствующая диффузно√му отражению: −sgn(ν) αk .Альтернативный подход в реализации метода Монте-Карло основывается на другом аналитическом представлении решения задачи (4.39)-(4.42).
Всоответствии с [125] из уравнения (4.39) следует равенствоZ01f (z, ν)νdν−11= 2(1 − λ) Θ4 (z) −2Z1f (z, ν)dν .−1После подстановки этого представления в (4.41) и последующего интегрирования мы получим1−λΘ(z) =NcZ zZ00ζZ1 14Θ (x) −f (x, ν)dν dxdζ + C1 z + C2 .2 −1(4.55)Постоянные C1 и C2 определяются из граничных условий (4.42) следующимобразом:C1 = τ −11−λΘ2 − Θ1 −NcZ0τZ0ζ1Θ4 (x) −2Z1f (x, ν)dν dxdζ ,−1C2 = Θ1 .На основе (4.43) и (4.55), используя метод Монте-Карло, мы получаем следующую аппроксимацию для температуры на j-ом шаге итерационной процедуры:hjihj−1iΘ (z) = (1 − α)Θα(z)+!MX4z(1 − λ)hjihj−1ixkΘ(zk ) − f N (zk , νk ) + C1 z + C2 ,M Nc(4.56)k=1где xk , zk , νk есть независимые реализации случайных величин на соответ-147hjiствующих интервалах (0, z), (0, xk ) и (−1, 1).
Вычисление f N (z, ν) осуществляется на основе формул (4.51), (4.52), предполагая, что температура равhj−1iна Θ(z). Постоянная C1 также вычисляется с помощью метода МонтеhjiКарло. В итоге, приближение Θ (z) функции Θhji (z) вычисляется на основеформул (4.51), (4.52) и (4.56).Таким образом, построены два вида рекуррентных соотношений, основанных на методе Монте-Карло.
Предложенные подходы позволяют избежатьнеустойчивости, вызванной операцией дифференцирования в правой частиравенства (4.41).8.3.Диффузионное приближение моделирадиационно-кондуктивного теплообменаРассмотрим подход, основанный на использовании диффузионного приближения (P1 приближения) для уравнения переноса излучения (4.39) [103].С этой целью представим функцию f (z, ν) двумя первыми слагаемыми ееразложения в ряд Фурье по полиномам Лежандра:f (z, ν) ' Φ0 (z) + νΦ1 (z).(4.57)Это дает следующую аппроксимацию уравнения (4.39):−Φ000 (z) + 3(1 − λ)Φ0 (z) = 3(1 − λ)Θ4 (z),(4.58)Как уже отмечалось в §7, существуют различные подходы для вывода граничных условий для диффузионного приближения (4.58). Остановимся наследующих двух.
В первом, для построения граничных условий мы подставляем выражение (4.57) в граничные условия (4.40) вместо функции f (z, ν)и интегрируем (4.40) по всем "входящим" в слой направлениям ν. Это даетследующие условия:411 + ρs1 + ρd1 Φ00 (0) = ε1 Θ41 ,ε1 Φ0 (0) −23(4.59)14814 dsε2 Φ0 (τ ) +1 + ρ2 + ρ2 Φ00 (τ ) = ε2 Θ42 .23(4.60)Во втором случае, будем использовать граничные условия Маршака [142],которые дают21 + ρs1 + ρd1 Φ00 (0) = ε1 Θ41 ,(4.61)321 + ρs2 + ρd2 Φ00 (τ ) = ε2 Θ42 .ε2 Φ0 (τ ) +(4.62)3Заметим, что использование весовой функции (4.1), применяемой для оптиε1 Φ0 (0) −чески толстых слоев, дает следующие граничные условия:16 d17s1 + ρ1 + ρ1 Φ00 (0) = ε1 Θ41 ,ε1 Φ0 (0) −241717161 + ρs2 + ρd2 Φ00 (τ ) = ε2 Θ42 ,ε2 Φ0 (τ ) +2417которые близки к граничным условиям Маршака (4.61), (4.62) и дают в сравнении с ними визуально неразличимые результаты.В новых обозначениях перепишем уравнение (4.41) в следующем виде:Θ00 (z) = −σΦ000 (z),σ=1.3Nc(4.63)Заметим, что функция Φ0 (z) интерпретируется как функция f (z, ν), усредненная по всем направлениям ν.Из уравнения (4.63) мы получаемΘ(z) = −σΦ0 (z) + C1 z + C2 .(4.64)Константы C1 , C2 определяются из граничных условий (4.42):C1 = τ −1 (Θ2 − Θ1 + σ(Φ0 (τ ) − Φ0 (0))),C2 = Θ1 + σΦ0 (0).Таким образом, задача (4.39)-(4.42) сводится к системе уравнений (4.58)-149(4.60) и (4.64) (или (4.58), (4.61), (4.62) и (4.64)).В следующем разделе представим результаты численных экспериментов,реализующих предложенные подходы.8.4.Численное моделирование радиационно-кондуктивноготеплообменнаЧисленные эксперименты выполнены для двух расчетных примеров, представленных в [159, 160], где использовался метод простой итерации [159] иитерационный метод Ньютона [160] в комбинации с устойчивой версией PNприближения.