Диссертация (1145283), страница 36
Текст из файла (страница 36)
Дальнейшее уменьшение шагаздесь также не привело к существенному выигрышу в расчетах, а лишь значительно увеличило время счета по этой схеме.Разработанный комплекс программ по схемам 4.3.5.1. – 4.3.5.10. (схемы3.5.2–3.5.5, 3.5.7, 3.5.8 приведены в приложении С) и численные расчеты поним задачи Коши (4.30), (4.34) (или эквивалентной системы (4.48)) показали,что для исследуемой задачи ни одна из этих схем не приводит к успеху.
Аименно, начиная с некоторого достаточно короткого промежутка времени, всеэти схемы дают срыв. Поэтому возникла необходимость в создании другихвычислительных алгоритмов решения уравнения (4.30) (и эквивалентной системы (4.48)). Достаточно эффективной в рассматриваемой задаче оказаласьпредложенная ниже модифицированная явная схема [125].4.3.6. Новая модифицированная явная схема численного решенияжесткого нелинейного неавтономного обыкновенногодифференциального уравнения второго порядкаРазрешим уравнение (4.30) относительно R̈ и запишем его в видеR̈(t) + Ṙ(t)2 f1 (R(t)) + Ṙ(t)f2 (R(t)) = f3 (R(t), t) ,(4.53)h(t) h(t)2 − 4h(t) + 6Am2h(t)−3h(t)+3,, f2 (R(t)) =f1 (R(t)) =2R(t)R(t)2238Ap Φ(t)As (2 − h(t))+,h(t)R(t)2h(t)R(t)1h = h(R(t)) = 1 −.(1 + RK3 )1/3f3 (R(t), t) = −Разностный аналог уравнения (4.53) запишем следующим образом:Rn+1 − 2Rn + Rn−1 + (Rn+1 − Rn−1 )2 f1 (Rn+1 )/4+(4.54)+(Rn+1 − Rn−1 )τ f2 (Rn )/2 = τ 2 f3 (Rn , tn+1 ).Это квадратное уравнение относительно Rn+1 , физический смысл имеет лишьположительный корень, равныйRn+1b=− +2b = −2Rn−1 +2c = Rn−1+b2−c41/2,4f2 (Rn )+ 2τ,f1 (Rn )f1 (Rn )f3 (Rn , tn+1 )f2 (Rn )4(Rn−1 − 2Rn )− 2τ Rn−1− 4τ 2.f1 (Rn )f1 (Rn )f1 (Rn )Найденное явное решение уравнения (4.54), позволяет записать простойалгоритм перехода к следующему временному слою.
Величина y = Ṙ определяется по формулеRn+1 − Rn−1.2τПредложенная схема имеет второй порядок аппроксимации по τ .yn+1 =В модифицированной явной схеме 3.6. отклонения ∆R численного решения задачи Коши (4.30), (4.34) от известного точного решения во все моментывремени значительно меньше, чем во всех рассмотренных схемах. Величина∆R для разных моментов времени приведена в таблице 4.5. Шаг по времени,как и в представленных выше схемах, составил величину τ = 10−5 .Таблица 4.5t0.000150.00030 0.50000 1.00000 1.50000 56.8850 68.0650∆R 5.62e-13 2.00e-12 2.61e-62.00e-65.20e-6 -2.43e-3-0.462394.3.7.
Выводы из проведенных расчетов решения прямой задачи вначале процесса расширенияРасчеты по созданному комплексу программ решения жестких нелинейных неавтономных обыкновенных дифференциальных уравнений применительно к уравнению (4.30) позволяют сделать следующие выводы.1. Для всех вычислительных алгоритмов, как и следовало ожидать, устойчивость счета достигалась только при весьма малых значениях шага интегрирования τ . Величина шага не превышала величины обратно пропорциональной максимальному (по t) числу обусловленности матрицыЯкоби (4.49) и составляла величину τ ≤ 10−5 .2. Для всех вычислительных алгоритмов требовалась большая точность вычислений.3. Наименьшие отклонения численного решения уравнения (4.30) от известного из обратной задачи точного решения достигались в предложенномалгоритме 3.6., для него и точность расчета была выше и расчет оказалсяустойчивым на значительно большем интервале времени процесса расширения. Конечно, успех предложенного алгоритма базируется на использовании конкретного вида уравнения (4.30).
Метод не является универсальным, однако в рамках рассматриваемого класса задач о расширениижидкого сферического слоя, он обладает неоспоримым преимуществом.4. Ни одна из опробованных схем не позволила для набора параметров задачи, представляющих практический интерес, получить численное решениеуравнения (4.30) на всем интервале времени [0, tk ].Проведенные исследования позволили построить вычислительный алгоритм численного решения прямой задачи на первом этапе в широком диапазоне исследованных параметров процесса. Решение задачи на остальноминтервале времени представлено в четвертом параграфе.2404.4.
Асимптотическое решение приближенногосингулярно-возмущенного уравнения, моделирующего динамикурасширения жидкого слоя при t > t∗4.4.1. Решение прямой задачи на конечной стадии процессарасширенияКак отмечалось, при t → tk функция h(t)h(t) = h(R(t)) = 1 −1(1 +K 1/3R3 )=H(t)R̂(t)стремится к нулю. Это приводит к появлению малого функционального параметра при старшей производной в уравнении (4.30). Известно [118], чтосвойства таких уравнений близки к свойствам сингулярно-возмущенных систем, однако в отличие от последних, здесь малым является функциональныйпараметр, который сам зависит от искомого решения. Некоторых успеховудалось достичь в интегрировании уравнения (4.30) в конце процесса расширения жидкой оболочки в работах [129].
Там была использована концепцияметода функционального параметра, сформулироваyная в цикле работ А. Н.Панченкова [131].Более эффективным оказался подход, не содержащий трудностей обоснования, присущих методу малого функционального параметра. Он позволиланалитически рассчитать окончание процесса расширения жидкого слоя с достаточной точностью.Заметим, что в исследуемой задаче есть малый постоянный параметр δ:δ=K,Rk3(4.55)равный отношению объема жидкости 34 πK к объему сферы конечного радиуса43πRk3 .
Практический интерес представляют лишь те задачи, в которыхвыполняется условиеK ≪ Rk3 → δ ≪ 1.Выделим в функции h(R) параметр δ, записав ее в виде: 3 !−1/31Rk.=1−1+δh(R) = 1 −R(1 + RK3 )1/3(4.56)241В начале процесса величина безразмерного радиуса R(t) много меньше егоконечного значения Rk и величина γγ=δRkR3не мала даже для весьма малых значениях δ. Однако, в конце процесса приR → Rk выполняется условие:γ → δ ≪ 1.Разложим функцию h(R) (4.56) в степенной ряд по γ и для γ ≪ 1 огра-ничимся членами первого порядка малости по γ. В результате приходим кследующему преставлению:δγh(R) ≃ =33RkR3.(4.57)Перейдем при t → tk от решения точного уравнения (4.30) к решению прибли-женного уравнения, которое получается из (4.30) заменой функции h(R) ееприближенным значением (4.57).
В первом приближении по γ это уравнениеимеет вид:Ṙ6As R δAs 3Ap Φ(t)R2δ R̈ + δ3Am 2 = − 3 + 2 +.RRkRRk3(4.58)Оно аппроксимирует исходное уравнение (4.30) с точностью до ε1 (t), равным:δ Rk3.ε1 (t) =3 R(t)3При заданной малой величине ε∗1 значение радиуса Ra , начиная с которого допустим переход от (4.30) к (4.58), определяется равенством: 3 1/3 1/3δRkKRa ==.3ε∗13ε∗1(4.59)Допустимы только управляющие функции Ф(t), которые обеспечиваютмонотонное возрастание функции R(t) на всем интервале [0, tk ].Момент времени ta , для которого радиус достигает значения Ra (4.59),определяется в ходе решение прямой задачи на начальном этапе из условия:R(ta ) = Ra .(4.60)242Для характерных значений безразмерного комплекса K, представляющих практический интерес при создании космических зеркал, при ε∗1 порядка10−3 интервал [ta , tk ] был значительно больше начального интервала [0, ta ].Асимптотическое решение приближенного уравнения (4.58)Приведем асимптотическое решение задачи Коши для уравнения (4.58)на некотором интервале [t∗ , tk ] при начальных данныхR(t∗ ) = R∗ ,Ṙ(t∗ ) = Ṙ∗ .(4.61)Алгоритм выбора момента времени t∗ приведен далее (п.
4.4.2). Величины R∗ ,Ṙ∗ в момент времени t∗ рассчитываются в ходе решения прямой задачи дляуравнения (4.30) с начальными данными (4.34) методами, рассмотренными впараграфе 4.3 (схема 3.6).Утверждение4.◮Асимптотическоерешениеприближенногосингулярно-возмущенного уравнения (4.58) на интервале [t∗ , tk ] приначальных данныхR(t∗ ) = R∗ ,Ṙ(t∗ ) = Ṙ∗имеет видR(t) = b0 (t) + tC1 + δ(C2 + b1 (t)),Ṙ(t) = ḃ0 (t) + C1 + δ ḃ1 (t),b0 (t) =2As,Ap Φ(t)b1 (t) = Rk3 (b20 b̈0 + 3Am ḃ0 − As )/(6As b20 ),C1 = Ṙ∗ − ḃ0 (t∗ ) − δ ḃ1 (t∗ ),C2 = (R∗ − b0 (t∗ ) − t∗ C1 )/δ − b1 (t∗ ). ◭Доказательство. ⊲ Уравнение (4.58) содержит малый постоянный параметр при старшей производной и является сингулярно-возмущенным. Исследования асимптотических решений задачи Коши для таких уравнений восходят к работам А.Н.Тихонова [132] и А.Б.
Васильевой [133]. Следуя этимработам, запишем асимптотическое решение уравнения (4.58) в виде суммы243рядов по «быстрому»времени τ = t/δ и по «медленному» времени t (δопределено формулой (4.55) )R(t) = δr(τ ) + b(t),r(τ ) =∞Xkδ rk (τ ),b(t) =k=0∞Xδ k bk (t).k=0Выразим первую и вторую производные Ṙ, R̈ через производные функцийr(τ ) и b(t)Ṙ =dr db+= r′ + ḃ,dτdt1 d2 r d2 b 1 ′′+= r + b̈.δ dτ 2 dt2δВ первом приближении по δ для функции R(t) и ее производных спраR̈ =ведливы следующие выражения:R = b0 (t) + δ(r0 (τ ) + b1 (t)),Ṙ = ḃ0 (t) + r0′ (τ ) + δ ḃ1 (t),1R̈ = r0′′ (τ ) + b̈0 (t) + δ b̈1 (t).δ(4.62)Здесь и далее точкой обозначено дифференцирование по t, штрихом –по τ . Представим функцию R(t) (4.62) следующим образом:R = b0 (1 + η),(4.63)η = δ(r0 + b1 )/b0 ,b0 = b0 (t),b1 = b1 (t),r0 = r0 (τ ).Для малых величин η функции R2 и 1/R2 , входящие в уравнение (4.58),можно разложить в степенные ряды.
В первом приближении по η эти разложения имеют вид:11≃(1 − 2η),R2b20R2 ≃ b20 (1 + 2η).Запишем уравнение (4.58), используя разложения (4.62), (4.64),r0 (τ )′′ + δ b̈0 (t) + 3Am δr0 (τ )′6Asḃ0 (t)+3Aδ=−b0 (t)+mb0 (t)2b0 (t)2Rk3(4.64)244+6AsAs6As3Ap Φ(t)2δ−b(t)+δr(τ)−δb1 (t)+00Rk3b0 (t)2Rk3Rk3+(4.65)6Ap Φ(t)6Ap Φ(t)b(t)r(τ)δ+b0 (t)b1 (t)δ.00Rk3Rk3Приравняем коэффициенты в левой и правой частях уравнения (4.65)при одинаковых степенях δ, причем отдельно коэффициенты, зависящие от τ ,и коэффициенты, зависящие от t.