маслов. Практикум предназначен для проведения лабораторных занятий и са мостоятельной работы студентов факультета Ф, обучающихся по спе циальности Ядерные реакторы и энергетические установки
 Скачать 1.6 Mb.
|
|
2.4. Интегроинтерполяционный метод (метод баланса) построения разностных схем 2.4.1. Свойство консервативности разностной схемы В приведенных выше разностных схемах при аппроксимации пространственного дифференциального оператора производные в дифференциальном уравнении и граничных условиях просто заме- нялись конечными разностями. Такой подход не всегда приводит к успеху. Для более сложных задач, описываемых нелинейными уравнениями и уравнениями с переменными коэффициентами, за- мена производных конечными разностями может привести к схе- мам, которые будут иметь большую погрешность, либо вообще окажутся непригодными для счета. Исходные дифференциальные уравнения получаются на основе записи законов сохранения энергии, массы, количества движения для произвольного объема сплошной среды. Очевидно, что для по- лучения разностного решения, хорошо описывающего реальный процесс изменения искомых функций в количественном и качест- 43 венном отношениях, целесообразно потребовать выполнения этих законов и для разностного решения. Для непрерывного точного решения закон сохранения выполня- ется для произвольной области тела. Для разностного решения тре- бование выполнения закона сохранения имеет важную особен- ность, обусловленную дискретным разбиением тела. Поскольку разностное решение ищется в отдельных точках тела, то необходи- мо разбить тело на такое же число элементарных объемов, каждый из которых будет включать одну точку, а затем потребовать вы- полнения закона сохранения как для произвольного элементарного объема так и для любой области, составленной из этих элементар- ных объемов. Последнее требование будет выполнено, если обес- печить согласование тепловых потоков для любых соседних объе- мов, заключающееся в равенстве значений протекающих через об- щую границу тепловых потоков. Желательно точное выполнение сформулированных условий при конечном разбиении расчетной области, а не только при стремлении максимального размера элементарной области к нулю. Это позволяет получать правдоподобные решения даже на грубых сетках. Разностные схемы, при которых получаются численные реше- ния, удовлетворяющие законам сохранения, называются консерва- тивными. Консервативность схемы не обеспечить без принятия специаль- ных мер. Поэтому в большинстве случаев разностные уравнения получают не из аппроксимации операторов дифференциального уравнения, а из непосредственной аппроксимации самих соотно- шений теплового баланса, записанных для элементарных объемов. При этом для тепловых потоков на границах используются выра- жения, обеспечивающие выполнение условий согласования. Этот способ построения консервативных разностных схем называется интегроинтерполяционным методом или методом баланса. 2.4.2. Метод баланса Основные этапы применения метода баланса: 44 − область, в которой ищется решение, разбивается на элементар- ные объемы (элементарные ячейки), построенные вокруг каж- дого узла сетки; − для всех внутренних и граничных ячеек записываются уравне- ния теплового баланса, включающие значения тепловых пото- ков на границах ячеек; при записи уравнений баланса для яче- ек, прилегающих к границам, используют граничные условия; − аппроксимируются члены, входящие в уравнения теплового баланса, выражая их через значения сеточной функции; при этом выражения для тепловых потоков должны удовлетворять условию согласования. Поскольку число ячеек равно числу узлов пространственного разбиения, то в результате этих действий получается полная систе- ма алгебраических уравнений – разностную схема, при решении которой можно определить разностное решение. Рассмотрим построение консервативной разностной схемы в случае нестационарного уравнения для стержня с боковым тепло- обменом: ) , ( ) ( ) ( t x q T x x T x x t T c v v + α − ⎥⎦ ⎤ ⎢⎣ ⎡ ∂ ∂ λ ∂ ∂ = ∂ ∂ ρ . (2.28) Выберем неравномерную пространственную сетку { } N n j n x 1 = с ша- гом h n =x n+1 – х п . Элементарные ячейки для всех внутренних узлов х n построим, отступая от каждого узла на половину шага влево и вправо (рис. 2.2,а). Элементарная ячейка для узла х п представляет собой отрезок [x n-1/2 , x n+1/2 ], где 2 / 2 / 1 n n n h x x + = + , 2 / 1 2 / 1 − − − = n n n h x x 45 x n-1 x n x n+1 (h n-1 +h n )/2 x 1 x 2 h 1 /2 α v q 0 q 3/2 α 0 q v а) б) Рис. 2.2. Внутренняя (а) и граничная (б) элементарные ячейки разбиения Уравнение теплового баланса для ячейки [x n-1/2 , x n+1/2 ] на проме- жутке времени от t j-1 до t j имеет вид [ ) ( ) ( 2 / 1 2 / 1 1 2 / 1 2 / 1 2 / 1 2 / 1 1 dt dx T q q q dx T T c n n j j n n x x v v t t n n x x j j ⎥ ⎥ ⎦ ⎤ α − + + − = − ρ ∫ ∫ ∫ + − − + − − + − (2.29) Здесь выражение в левой части представляет собой количество те- пла, идущее на нагрев элементарной ячейки. Члены в правой части имеют следующий смысл 2 / 1 ) ( 2 / 1 ± λ − = ± n x n dx dT x q – тепловые пото- ки на границах ячейки, интеграл по отрезку – баланс тепловыделе- ния и бокового оттока от ячейки. Аппроксимируем левую часть (2.29) выражением 2 ) ( ) )( ( ) ( 1 1 1 2 / 1 2 / 1 − − − + ρ − ≈ − ρ ∫ + − n n n j n j n x x j j h h c T T dx T T c n n (2.30) В правой части равенства (2.29) для аппроксимации интегралов по пространственной переменной используем следующие выраже- ния. 46 При вычислении теплового потока с боковой поверхности будем считать, что температура не изменяется на отрезке [x n-1/2 , x n+1/2 ], т.е. 2 / 1 2 / 1 2 / 1 2 / 1 ∫ ∫ + − + − α ≈ α n n n n x x v n x x v dx T Tdx (2.31) Приближение для тепловых потоков 2 / 1 ± n q получим в предпо- ложении о малом изменении потока q(x) на соответствующих ин- тервалах. Очевидно, что при малых h поток мало изменяется даже в случае разрыва λ(х). Из закона Фурье имеем ) ( ) ( или ) ( x x q dx dT dx dT x q λ − = λ − = . (2.32) Проинтегрировав равенство (2.32) по отрезку [x n-1 , x n ], получим: ∫ ∫ ∫ − − − λ − λ − − ≈ − = − = n n n n n n x x x dx n x x x x q n n x x q dx T T dx dx dT 1 1 1 ) ( 2 / 1 ) ( ) ( 1 . (2.33) Таким образом, тепловые потоки через границы элементарной ячейки выражаются через разности температур в узлах: , , 1 1 2 / 1 2 / 1 1 1 2 / 1 2 / 1 − + + + − − − − − λ ≈ − λ ≈ n n n n n n n n n n h T T q h T T q (2.34) где 2 / 1 − λ n , 2 / 1 + λ n – эффективные теплопроводности отрезков [x n-1 , x n ] и [x n , x n+1 ], соответственно , 1 ) ( 2 / 1 1 ) ( 1 2 / 1 1 1 − λ + − λ − − ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ = λ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ = λ ∫ ∫ + − n n n n x x x dx n n x x x dx n n h h (2.35) Очевидно, что аппроксимация (2.34) удовлетворяет условию со- гласования потоков. 47 Для аппроксимации интегралов по времени в правой части вы- ражения (2.29) примем, что при расчете интегралов можно заме- нить изменяющуюся на отрезке [t j-1 , t j ] температуру на постоянное в каждой точке этого отрезка средневзвешенное значение , ) 1 ( ) ( 1 − σ − + σ ≈ j n j n n T T t T (2.36) где σ – параметр, который может принимать значение от 0 до 1. Тогда при вычислении интеграла по времени получаем следующее выражение: ]. ) 1 ( [ ) ( 1 1 − σ − + σ ∆ ≈ ∫ − j n j n t t n T T t dt t T j j (2.37) Подставляя выражения для тепловых потоков в уравнение ба- ланса (2.29), получим разностную схему [ ] , ) 1 ( ) 1 ( ) ( 1 1 1 1 1 2 / 1 1 1 1 2 / 1 1 1 2 / 1 1 2 / 1 1 − − − − − − − − + + − − − + + − σ − + σ α − + ⎥ ⎥ ⎦ ⎤ − λ + ⎢ ⎢ ⎣ ⎡ − λ σ − + ⎥ ⎥ ⎦ ⎤ − λ + ⎢ ⎢ ⎣ ⎡ − λ σ = ∆ − ρ j n j n n vn n j vn n j n j n n n j n j n n n j n j n n n j n j n n j n j n n u u l l q h u u h u u h u u h u u t u u cl (2.38) где 2 ) ( 1 − + = n n n h h l , , ) ( 1 2 / 1 2 / 1 ∫ + − α = α n n x x v n vn dx x l 48 ∫ ∫ − + − ∆ = j j n n t t x x v n j vn dt dx t x q t l q 1 2 / 1 2 / 1 ) , ( 1 Уравнения (2.38) записываются для всех внутренних точек тела. При σ = 0 получаем явную схему, при σ = 1 – неявную, при σ = ½ – схему Кранка – Николсона. В случае непрерывности функций λ(х), α v (х), q v (x, t)интегралы от них обычно заменяют простейшими квадратурными формулами. Для q v и α v чаще всего полагают ) ( , ) 2 / ) ( , ( 1 n v vn j j n v j vn x t t x q q α ≈ α + ≈ − , (2.39) т.е. заменяют соответствующие интегралы формулой прямоуголь- ников. Эффективные теплопроводности 2 / 1 ± λ n вычисляют одним из следующих способов: а) ) ( 2 / 1 2 / 1 ± ± λ = λ n n x ; (2.40) б) 2 / )] ( ) ( [ 1 2 / 1 n n n x x λ + λ = λ ± ± ; (2.41) в) )] ( ) ( /[ ) ( ) ( 2 1 1 2 / 1 n n n n n x x x x λ + λ λ λ = λ ± ± ± . (2.42) Если функции λ, α v , q v имеют разрывы между узлами, то для по- вышения точности разностной схемы, как правило, следует вычис- лять интегралы точно. Особенно это существенно в случае много- мерных задач, когда приходится вести расчет при достаточно гру- бых сетках. Аппроксимация граничных условий. Возьмем для определенно- сти элементарную ячейку [0, h 1 /2], прилегающую к границе х = 0 (рис. 2.2,б). При записи закона сохранения энергии для элементар- ной ячейки используем чисто неявную схему (при σ = 1), а также выражения (2.34) для тепловых потоков. Тепловой поток, выходящий из ячейки через границу х = h 1 /2, равен 49 1 1 2 / 3 1 2 1 2 / 3 2 / 3 2 1 ) ( где , − ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ λ = λ − λ ≈ ∫ x x j j x dx h h T T q , (2.43) тепловой лоток, рассеиваемый в среду на границе, ра- вен j I T q 1 0 α = ;выделяемая внутренними источниками мощность – 2 / 1 1 h q q j v II = ;рассеиваемый с боковой поверхности тепловой по- ток – 2 / 1 1 1 h T q j v III α = ; расходуемая на нагрев элементарного объ- ема мощность t T T h c q j j IV ∆ − ρ = − ) ( 2 1 1 1 1 Из закона сохранения энергии следует 0 0 2 / 3 = − − + + − − IV III II I q q q q q q или 0 2 1 1 1 1 1 1 1 0 1 0 1 1 2 2 / 3 = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ∆ − ρ − α − + + α − − λ − t u u c u q h q u h u u j j j v j v j j j . (2.44) Первые три слагаемых (2.44) совпадают с простейшей аппрок- симацией граничного условия (2.11), полученной простой заменой производной конечной разностью. Дополнительные слагаемые учитывают действие внутренних источников, теплообмен с боко- вой поверхности и затраты теплоты на нагрев элементарной ячей- ки. Эти слагаемые пропорциональны h 1 поэтому при h 1 → 0 обе аппроксимации граничного условия становятся идентичными. 50 Можно показать, что погрешность аппроксимации граничного ус- ловия уравнением (2.44) – о( 2 1 h ), а уравнением (2.11) – о( 1 h ). Аналогичным образом строится разностная аппроксимация гра- ничного условия при х = l. Система уравнений (2.38) для внутренних точек п = 2, ..., N – 1 и уравнений типа (2.44) для граничных точек представляет собой консервативную неявную схему численного решения задачи (2.28), (2.2). 2.5. Метод прогонки При использовании неявных разностных схем, например (2.38), на каждом временном слое необходимо решать систему алгебраи- ческих уравнений с числом неизвестных N, которое может быть достаточно велико. Одним из лучших прямых методов решения систем линейных алгебраических уравнений общего вида является метод последова- тельного исключения Гаусса с выбором главного элемента. Расчет по формулам этого метода требует примерно N 3 арифметических операций, поэтому при достаточно больших N потребуются значи- тельные затраты машинного времени. Особенность системы (2.38), (2.44) состоит в том, что в каждое уравнение для внутренних точек входят по три неизвестных, номе- ра которых отличаются на единицу, а в первое и последнее уравне- ния для точек п = 1 и п = N – по два «соседних» неизвестных. Если учесть такой специфический вид построенной нами системы разно- стных уравнений, то эффективность алгоритма ее решения можно существенно повысить. Запишем систему уравнений в следующем каноническом виде: для граничной точки п = 1 0 1 1 1 2 1 = + + d u b u a j j ; (2.45) для внутренних точек п = 2, ..., N – 1 51 0 1 1 = + + + − + n j n n j n n j n n d u c u b u a ; (2.46) для граничной точки п = N 0 1 = + + − N j N N j N N d u c u b . (2.47) Выражения для коэффициентов а п , b n , c n , d n легко получить из соответствующих уравнений разностной схемы. Трехдиагональный вид матрицы системы уравнений (2.45) – (2.47) позволяет упростить вычисления по методу Гаусса так, что- бы не проводить операции с нулевыми элементами. Тем самым объем вычислений удается значительно уменьшить. Модификация метода Гаусса для системы уравнений с трехдиагональной матри- цей называется методом прогонки. Алгоритм расчета по методу прогонки: 1) определение коэффициентов f 1 и g 1 1 1 1 1 1 1 / , / b d g b a f − = − = ; (2.48) 2) определение коэффициентов f n и g n при п = 2, ..., N – 1 1 1 1 , − − − + + − = + − = n n n n n n n n n n n n f c b g c d g f c b a f ; (2.49) 3) определение u N ) /( ) ( 1 1 − − + + − = N N N N N N N f c b g c d u ; (2.50) 4) определение u n для п = N – 1,...,1 n n n n g u f u + = +1 . (2.51) Для решения системы (2.45) – (2.47) по методу прогонки требу- ется примерно 9N арифметических действий, т.е. значительно |