Диссертация (Математическое моделирование теплофизического эксперимента на основе численных методов расщепления и идентификации), страница 11

Обозначив эти∂y ⎝ ∂y ⎠∂ ⎛ ∂ ⎞⎜λ ⎟ ,∂x ⎝ ∂x ⎠операторы Λ x , Λ y , запишем (2.9) в более простом виде:cV ijql ijU τ +1 − U τ= Λ xU ij + Λ yU ij + 2 .Δtπ r0(2.10)Схема (2.9) имеет порядок аппроксимации O(Δt + h 2 ) на точном решении T ( xi , y j ,τΔt ) ,где обозначено h 2 = (Δxi−)2 + (Δxi+) 2 + (Δy −j )2 + (Δy +j ) 2 , и верна для регулярных узлов < i, j > , новблизи границ контакта областей, где имеем разрывы 1-го рода теплофизических игеометрических характеристик, должны выполняться условия сопряжения (2.3)–(2.4),поэтому схема должна быть соответствующим образом модифицирована (см.

п. 2.5). Длярасчета температурного поля согласно схеме (2.9) необходимо также выполнение краевогоусловия (2.5), что возможно реализовать естественным образом с помощью принятия заграничный контур так называемой «неотражающей границы» (см. п. 2.5). Начальное условие(2.6), очевидно, записывается для временного слоя t 0 = 0 в виде U ij0 = T0 = Const . ЗначенияcV , ql в (2.9) вычисляются в зависимости от текущей координаты расчетной области:cV ij = cV ( xi , y j , tτ ) , ql ij = ql ( xi , y j , tτ ) , причём ql ij = 0 , если ( xi , y j ) ∉ Q0 , и ql ij = ql , если( xi , y j ) ∈ Q0 ,азначениякоэффициентовλтеплопроводностивычисляютсявсоответствующих «полуузлах», например λ i −1/ 2, j = λ ( xi − hxi −1 / 2, y j , tτ ) .2.2.4. Метод и алгоритм вычисленийРассмотримразностнуюсхему(2.9),котораяможетбытьявной,неявной,комбинированной, и реализована разными методами.

Если рис. 2.2. соответствует текущемувременному слою tτ , то разностная схема (2.9) будет явной, и все элементы формулы вправой части равенства, а также cV ij записываются с верхними индексами τ , а еслиследующему временному слою tτ +1 , то – неявной, и записываются индексы τ + 1 . Прииспользовании явной схемы, неизвестная величина U ijτ +1 явно выражается через известныеU ijτ , U iτ−1, j , U iτ+1, j , U iτ, j −1 , U iτ, j +1 и вычисляется для каждого узла < i, j > отдельно. Прииспользовании неявной схемы, разностная задача представляет собой систему линейныхалгебраических уравнений (СЛАУ) относительно неизвестных величин U ijτ +1 , U iτ−+1,1 j , U iτ++1,1 j ,U iτ, +j −11 , U iτ, +j +11 , i = 0, M , j = 0, N :5111⎛ λiτ++1/2,⎞λiτ−+1/2,2jjτ +1τ +1− ++UU⎜⎟⎟ ++−ijij1,1,Δxi + Δxi− ⎜⎝ Δxi+Δxi−⎠11⎡⎛ λiτ, +j +11/2 λiτ, +j −11/2 ⎞ cVτ +ij1 ⎤ τ +1⎛ λiτ++1/2,⎞λiτ−+1/2,22jj+++⎢ +⎥ Uij⎜⎟+⎟− ⎜+Δxi− ⎟⎠ Δy +j + Δy −j ⎜⎝ Δy +jΔy −j ⎟⎠ Δ t ⎥⎦⎢⎣ Δxi + Δxi ⎜⎝ Δxi⎛ λiτ, +j +11/2 τ +1 λiτ, +j −11/2 τ +1 ⎞ cVτ +ij1 τ qlτij+12U i , j +1 +U i , j −1 ⎟ =U −,− +⎜⎟ Δ t ij π r02Δy −jΔy j + Δy −j ⎜⎝ Δy +j⎠(2.11)решив которую, сразу получим значения U ijτ +1 во всех узлах < i, j > временного слоя tτ +1 .Схема (2.11) представляет собой разностный аналог уравнений (2.1), (2.2), отвечающийнеявному стандартному 6-точечному шаблону [71, Гл.9], изображенному на рис.

2.3а. ВСЛАУ (2.11) имеется ( M + 1) × ( N + 1) неизвестных, обозначенных как элементы двумерногомассива, тогда как в стандартной для СЛАУ индексации неизвестных используетсяодномернаяиндексация.Поэтомунеобходимопроизвестисквознуюнумерациюнеизвестных, например в виде U ijτ +1 = uk , k = 1 + i + j ⋅ M , то есть вдоль горизонталей.

Тогдаматрица данной СЛАУ будет иметь ленточный характер, с шириной ленты 2M − 1 иупорядоченной структурой ненулевых элементов внутри ленты. Очевидно, что сувеличением числа узлов в расчетной области, ширина ленты будет расти, что создаёттрудности при решении данной СЛАУ, так как число машинных операций равно KL2 , гдеK = ( M + 1) × ( N + 1) – количество узловых точек, L = 2M − 1 – ширина ленты [57, Гл.

3].Уменьшить число машинных операций позволяют экономичные методы, в которых числоопераций пропорционально K . Это, например, группа методов расщепления, [57, Гл. 3],[71, Гл.9], [151, Гл. 7], [163], [179] представленная методами: расщепления (на уровне задач),полного расщепления В.Ф. Формалёва–О.А. Тюкина, переменных направлений – МПНПисмена-Речфорда (Peaceman D.W., Rachford H.H.) и с экстраполяцией – МПНЭВ.Ф.Формалёва, дробных шагов – МДШ Н.Н.

Яненко – локально-одномерная схема.i, j +1слой τ+1слой τ+1i −1, j i, j hyi, j −1слой τi, j +1i +1, jΔti, ji, j −1слой τ+1/2hxi, jслой τhyΔt/2i, ji −1, jhxi +1, jΔt/2i, j(а)(б)Рис. 2.3. Шаблоны разностных схем./ – узлы с известным (текущим)/неизвестным значением сеточной функции;а – неявный 6-точечный; б – неявный 7-точечный с промежуточным слоем;52Разностная схема (2.10) в явном или неявном вариантах имеет свои преимущества инедостатки.

Явная схема достаточно проста: в отличие от неявной схемы, при том жепорядке аппроксимации [57, Гл. 3], она не требует решения на каждом шаге СЛАУ большойразмерности. Главный недостаток явной схемы – её условная устойчивость, то естьустойчивость при ограничении на величину шага по времени Δ t ≤ cV / (hx−2 + hy−2 ) / λ . Этотнедостаток становится особенно существенным для решаемой задачи (п.

2.2.2). Согласнофизической постановке (разд. 1.3), область Q0 очень мала (r0~1–10 мкм) и потребуетсоответствующей дискретизации, а значит и совсем маленького шага по времени (Δt~10–8 cдля типичных условий эксперимента). Таким образом, преимущество явной схемы в свободевыбора шага по времени для решения задач с быстро протекающими тепловыми процессами,в нашем случае практически отсутствует. Помимо этого, ввиду наличия источника Q0 ,значимым становится порядок обхода узлов сетки D . Неявная схема, несмотря надостаточную сложность реализации, является безусловно устойчивой, то есть устойчивойпри любых соотношениях шагов по времени Δ t и координатам hx , hy .

Существуют такжекомбинированные варианты схем, например, схема с весами, представляющая собойвыпуклую комбинацию явной и неявной схем. Комбинированная схема с весом равным 0.5называется схемой Кранка-Николсон, которая, как и неявная схема является безусловноустойчивой, однако при превышении некоторой величины шага по времени может нарушатьусловие монотонности разностной схемы [57, Гл. 3].Ввиду вышесказанного, для численного решения задачи в постановке п. 2.2.2 выберемнеявную разностную схему, которую реализуем экономичным методом переменныхнаправлений (МПН) [71, Гл.9], [151, Гл.7], [163].Согласно МПН проводится расщепление на уровне временных слоёв, путём введенияпромежуточного временного слоя tτ +1/2 = tτ + Δ t / 2 , а также на уровне задач (по x и по y),когда на каждом промежуточном временном слое пространственный дифференциальныйоператор по одному (текущему) направлению аппроксимируется неявно, а по другому –явно.Проведёмрасщеплениесхемы(2.11)указаннымспособомисоставимсоответствующий алгоритм вычислений.

Запишем сначала схему (2.11) в упрощенном виде,переобозначив её коэффициенты:⎛ xcVτ +ij1 ⎞ τ +1cVτ +ij1 τyx τ +1y τ +1Uij + bijτ +1−α U i −1, j − α U i , j −1 + ⎜ β ij + βij +⎟⎟ Uij − γ ijU i +1, j − γ ij U i , j +1 =⎜ttΔΔ⎝⎠xijτ +1yijτ +1(2.12)где обозначения коэффициентов α ijx , γ ijx , β ijx = α ijx + γ ijx , α ijy , γ ijy , β ijy = α ijy + γ ijy , bijτ +1 = qlτij+1 / π r02соответствуют (2.11) и все эти коэффициенты положительны.53Вместо 6-точечного шаблона рис. 2.3а принимаем 7-точечный шаблон рис.

2.3б спромежуточным временным слоем τ + 1/ 2 . Проводим вычисления согласно следующемуалгоритму.Алгоритм2.1. (Расчет температурного поля в системе «источник-образец-подложка»)Шаг 1. Задать расчетную область и соответствующую ей пространственную сеткуD = {xi , y j } , xi = x0 + ihx , y j = y0 + jhy , i = 0,1,..., M ,j = 0, N . Задать временную сеткуE = {tτ } , tτ = t 0 + τ Δt , t 0 = 0 , τ = 0,1, 2,..., P .Шаг 2. Вычислить U ij0 = Vij0 = T0 согласно (2.6). Положить τ := 0 .Шаг 3. Вычислить решение на (τ + 1) -м временном слое: Uijτ +1 в два этапа.Этап 1. Перейти со слоя τ на промежуточный (τ + 1/ 2) -й слой с шагом Δ t / 2 . Длякаждого фиксированного j ( j = 0, N ) методом прогонки [71, Гл. 1] решить 3-диагональнуюСЛАУ ( i = 0, M ):τ +1/2−α U i −1, jxij⎛ 2cVτ +ij1/2⎞ τ +1/2⎛ 2cVτ +ij1/ 2⎞xx τ +1/2y τ+⎜+ β ij ⎟ Uij − γ ijU i +1, j = α ij U i , j −1 + ⎜− β ijy ⎟ Uijτ + γ ijyU iτ, j +1 + bijτ +1/2 ,⎜ Δt⎟⎜ Δt⎟⎝⎠⎝⎠где все коэффициенты соответствуют (2.12), (2.11).В результате выполнения этапа 1 получаем решение на промежуточном (τ + 1/ 2) -мвременном слое – сеточную функцию U ijτ +1/ 2 .Этап 2.

Перейти со слоя (τ + 1/ 2) на (τ + 1) -й слой с шагом Δ t / 2 . Для каждогофиксированного i ( i = 0, M ) методом прогонки решить 3-диагональную СЛАУ ( j = 0, N ):⎛ 2cVτ +ij1⎞ τ +1⎛ 2cVτ +ij1⎞yy τ +1x τ +1/ 2−α U i , j −1 + ⎜+ β ij ⎟ Uij + γ ij U i , j +1 = α ijU i −1, j + ⎜− β ijx ⎟ Uijτ +1/2 + γ ijxU iτ++1,1/2j + bijτ +1 ,⎜ Δt⎟⎜ Δt⎟⎝⎠⎝⎠yijτ +1где все коэффициенты соответствуют (2.12), (2.11).В результате выполнения этапа 2 получаем решение на (τ + 1) -м временном слое –сеточную функцию Uijτ +1 .Шаг 4. Если τ < P , то положить n := n + 1 и перейти к шагу 3, иначе завершитьвычисления и принять сеточную функцию U ijP как решение задачи.Заметим, что при решении трёхдиагональных СЛАУ (шаг 3) вычислительный процессметода прогонки корректен и устойчив в силу преобладания диагональных элементов [71,Гл.