CC BY
62
свою деятельность) за 2002 г. Ставка исполнения составляет 4,96 %.
Напомним, что контракт, определяющий верхнее граничное значение процентной ставки, называется кэпом. Кэп состоит из серии кэплетов, представляющих собой контракты со сроками погашения платежа в момент времени Т+тна выплату разницы между рыночным значением ставки ЯТ в момент времени Т и значением ставки Х в течение периода т, основная сумма при этом составляет Ь. Главным условием является то, что данная разность должна быть положительна, то есть стоимость (ра-уо$) кэплета в момент Т+т составляет [2]:
шахЬ[(Л—Х),0].
На практике калибровка модели Халла-Уайта осуществляется путем выбора значений параметров, характеризующих поведение процентной ставки и волатильности таким образом, чтобы они соответствовали рыночным ценам опционов. Эмпирические значения параметра а принадлежат интервалу [0, 0,1], значение же параметра а принадлежат интервалу [0,01, 0,03] [5]. При калибровке начальные значения параметров а и а принимались равными 0,1 и 0,01 соответственно.
Параметры определяются при помощи метода наименьших квадратов таким образом, чтобы сум-
ма квадратов отклонений цен, полученных по модели Халла-Уайта, и рыночных цен опционов, была минимальной. В работе при выводе выражения применяли методы финансовой математики и теории опционов. Для проверки адекватности полученных результатов использовалась рыночная информация, основные расчеты и калибровка выражения были произведены с использованием пакета Mathematica.
В результате калибровки модели были получены следующие значения параметров: (а—>0,108114, о"—>0,0112018}. Зная параметры модели, можно рассчитать цену свопциона, дающего держателю право осуществлять платежи при ставке процента 6,2 % согласно условиям 3-летнего свопа, действие которого начинается через 5 лет. Волатильность процентной ставки свопа составляет 20 %. Платежи осуществляются раз в полугодие, основная сумма долга составляет 100 USD.
Применив (2), получим значение стоимости свопционного контракта - 2,0035. Рыночная стоимость данного свопционного контракта составляет 2,0038 (данные предоставлены компанией «Эконо-физика-Томск»). Таким образом, полученное выражение позволяет оценить стоимость свопционного контракта, при этом отклонение модельной цены от соответствующей рыночной не превышает 1 %.
СПИСОК ЛИТЕРАТУРЫ
1. Black F., Scholes M. The pricing of options and Corporate Liabilities // Journal of political economy. - 1973. - № 6. - P. 637-659.
2. Wilmott P. Derivatives: the theory and practice of financial engineering. - N.Y.: John Willey & Sons, 2000. - 742 p.
3. Hull J., White A. One-factor Interest-Rate Models and the Valuation of Interest-Rate Derivative Securities // The Journal of Financial and Quantitative analysis. - 1993. - № 2. - P. 235-240.
4. Hull J. Options, futures and other derivatives. - Toronto: University of Toronto, 2002. - 680 p.
5. Marshall J.F., Kapner K.P. Understanding swap finance. - OH: South-Western, 1990. - 250 p.
6. Keith R. An introduction to derivatives. - N.Y.: John Willey & Sons, 2000. - 198 p.
УДК 533.6.011
ИТЕРАЦИОННЫЙ МЕТОД РЕШЕНИЯ ОДНОМЕРНЫХ ДИВЕРГЕНТНЫХ УРАВНЕНИЙ ГАЗОВОЙ ДИНАМИКИ
В.М. Галкин
Томский политехнический университет E-mail: [email protected]
Предлагается итерационный метод решения дифференциальных уравнений, имеющих дивергентный вид и описывающих одномерное стационарное течение с переходом через скорость звука. Метод основан на использовании априорной информации о монотонном возрастании числа Маха вдоль сопла. Сравнение с другими методами проводилось на точном решении, а также при расчете двухфазного течения.
Введение
В [1] был предложен и далее в [2] уточнен итерационный метод решения одномерных стационарных уравнений газовой динамики, который за счет
использования априорной информации о монотонном изменении числа Маха внутри рассматриваемой области по скорости сходимости на порядок превосходит метод установления на основе схемы Маккормака [3] и имеет простой алгоритм. В
этом методе в качестве переменных используется полная энтальпия и энтропийная функция, а наличие в уравнениях правых частей позволяет учесть неравновесные процессы.
Однако, как отмечено в [4], методы, использующие уравнения, записанные в виде законов сохранения массы, импульса и энергии, имеют преимущества перед недивергентной формой записи. Кроме того, запись в дивергентной форме имеет более простой вид по сравнению с записью, используемой в [1]. Поэтому в данной статье рассматривается дальнейшее развитие метода [1]. Предлагаемый ниже подход, обладая простотой метода [1] и несколько уступая ему в скорости, использует дивергентную форму записи уравнений.
Математическая постановка задачи
Рассмотрим течение в сопле с переходом через скорость звука. Ставится задача нахождения параметров течения р, и, Р, удовлетворяющих одномерным стационарным уравнениям для идеального совершенного газа при наличии неравновесных процессов. Уравнения записаны в дивергентном виде [5]:
r = A(pU2 + P), R2 = ApUH,
dApU dx
= 0,
dA(pU2 + P) dx
dApUH = AC = A°2
= A'P + ACj,
(1)
L =-
YR)2
2R2Y2 -1)
(3)
(4)
При этом второе и третье уравнения в (1) примут вид:
= АР+ЛСи ^ = АС2. (5)
ах 1 ах 2
Обратный переход от Я1, К2, О к переменным р, и, Р производится по следующим формулам:
U =±
yR
- 2
(Y-1) R2
G (y+1)-G (y+ 1)) G(Y + 1)
p=G •P=A-p
(6)
где знак +(-) соответствует сверхзвуковому (дозвуковому) течению.
Из условий (2) следует, что существует единственная точка %=%,, в которой
(M2 -1)2 = min.
Очевидно, что при этом выполняются соотношения:
G = L(x*), x* = arg(min{L(x), x e (xa; xA)}), (7)
где: И=Ру/(р(у-1))+и2/2 - полная энтальпия; р, и, Р, у - плотность, скорость, давление и показатель адиабаты газа; А - площадь поперечного сечения сопла; х - продольная координата, принадлежащая рассматриваемой области [ха;х4]; С1 и С2 -известные, в общем случае нелинейно зависящие от параметров газа правые части уравнений движения и энергии, связанные с неравновесными процессами. Штрих обозначает производную по «х». Все величины безразмерны. Переход к безразмерным переменным производился путем отнесения: х - к половине размера минимального сечения, А - к площади в минимальном сечении шт(А), и -к критической скорости К, р - к критической плотности р., Р - к произведению р,К2, И - к квадрату критической скорости и.2.
Полагается, что задана площадь сопла А(х); на входе в сопло заданы граничные условия в виде И=И(ха), S=S(xa), где S=P/рr - энтропийная функция; С1(ха)=С2(ха)=0. В качестве априорной используется информация о том, что существует только одна внутренняя точка х., в которой число Маха М=1; внутри рассматриваемой области число Маха монотонно возрастает от дозвуковой до сверхзвуковой величины и не равно нулю:
М(х) ф 0, х е[ха; хь],
М(ха) < 1, М(х4) > 1, М(х.) = 1, х. е (ха; хД (2)
Вместо первого уравнения в (1) воспользуемся его интегралом Ари=О, где О - неизвестная константа. Введем обозначения:
dL
dx
= 0.
(8)
Таким образом, если в рассматриваемой области известно распределение Я1(х) и К2(х), то соотношения (4), (7) и (6) позволяют определить критическое сечение х., расход О и величины р, и, Р.
Поскольку распределение Я1(х) и К2(х) находится из решения задачи Коши для обыкновенных дифференциальных уравнений (5), то рассмотрим нахождение граничных условий для этих уравнений. Предположим, что расход О известен. Так как имеются И(ха), S(xa), С1(ха)=0, С2(ха)=0, то в сечении ха можно найти и, р и Р из соотношений:
^ ( л г ^ (хя)+- н (х и*-1= о
p =
AU
P = PYS (xa).
(9)
Подставляя найденные значения в (3), получим граничные условия Я1(ха) и К2(ха) для уравнений (5).
Численный алгоритм
С учетом выше изложенного для решения системы уравнений (1) предлагается следующий алгоритм:
1. Вводится расчетная сетка х1 и сеточные функции р, и, Р, Кц, К, ¿=0,1,...,к, где к - число точек сетки. Обозначим через верхний индекс у номер итерации.
2. Задаетсяу=1, О и начальное приближение р, и,, Р , удовлетворяющее (2).
3. С использованием (3) находится шах{Л1;, /=0,...,£} и присваивается всем Л^.
4. Переход к следующей итерации: у'=/+1. По значениям р, Ц, Р вычисляются правые части для уравнений (5).
5. Из (9) находятся граничные условия Л1(ха) и Л2(ха) для уравнений (5).
6. Из задачи Коши для обыкновенных дифференциальных уравнений (5) находятся Л1; и Л2; используется схема Эйлера второго порядка точности.
7. Для ускорения сходимости с использованием Л1; выполняется нижняя релаксация: Л()=Л("1)+ю(Л1,—Л(-1)), где ] - номер текущей итерации, - значение с предыдущей итерации, Щ) - будет использоваться на текущей итерации, Л1! - вычислено на текущей итерации, а - параметр релаксации, 0<а<1.
8. Из (7) определяется х,. Так как при использовании сеточных функций значение х, будет принадлежать дискретному множеству {х0,...,х£}, то более точное х, находится следующим образом. Пусть Х1=ш1п{/=1,...,^}. Параболическая интерполяция функции Ь по трем точкам и (8) дает: х = (х. - х1+1 )Lf_1 + (х.+1 - х1_1 )ЬГ + (х1_ 1 - х )11+1
* 2((х- х1+1)Ц_1 + (х.+1- х,.1)ь1 + (х..1- х )11+1) ' Далее, используя х,, уточняем О:
с = 4.1
(х* - х.)(х* - х1+1)
(х/-1 — х1 )(х/-1 — х1 +1)
+ г (х* х.-1)(х* х. +1) + ^ (х* х.- 1)(х* х )
+Ьт ' \+Ь1 +1 ,
(х. х/-1)(х/ х. + 1) (х. + 1
— х1.1)(х+1— х)'
9. Из уравнений (6) вычисляются р, Ц, Р; если х<х,, выбирается дозвуковое решение, в противном случае - сверхзвуковое.
10. При необходимости следующей итерации производится переход на пункт 4.
Сравнение с точным решением
Здесь и в следующем разделе рассматривалось течение в радиусно коническом сопле Лаваля, которое с учетом обезразмеривания описывалось зависимостью А(х)=у 2(х), где:
3,125, х < х, х = -0,5 -3,258^(0,785398) 2, 125 + - (х- х)2, х е (х;х2], х2 = х + 8т(0,785398) у(х) = 1,625 - х + 2х3, х е (х2; х3], х3 = -0,625»т(0,785398) , 1,625- ^0,6252 -х2, хе (хз;х4], х4 = 0,625»Ь(0,26 1799) 1,625 - 0,625 со8(0,26 1799) + (х - х4 )tg(0,26 1799), х > х4
ха=-4, хь=2, у=1,4, число точек сетки £=40, начальное поле параметров (пункт 2) находилось с использованием соотношения из [6]:
. . .. / , \(у+1)/(2(у-1))
т1п( А). = М' г+1
2 + (7-1) М2
С = нр (-
4 = н 0 [ N - А -
5'
ргБ'
(1 -7)5)'' 1 -7 5'
С =ин р IN - А.-
2 Н 1 НА (1 -7)5
здесь Н=АМ(х-х)/(х-х,)У+1-Ьг]; S=So[bl(x-xa)+1]; Н=Н1[[(8,/8)т/АГ0; Аа=А(ха); Ь2=0,5; Ь1=0,2; Н0=(7+1)/(2(7-1)); ^0=1//. Критическое сечение задано: х,=1.
Использование указанных правых частей позволяет найти точное решение в виде:
\н0
7+1 | _ т1п(Н)
м
и=М
2 + (7-1)М'
1/2
N
2(7-1) Н 2 + (7- 1)М:
Р =
т1п( Н)
Аи '
Р = 5ру.
На рис. 1 изображено положение точки х, в процессе решения: 1 - предлагаемым методом, при а=0,1; 2 - итерационным методом [1]; 3 - методом установления с использованием явной схемы Маккормака [3]. Для всех методов решение сходится к точному решению х,=1. Предлагаемый метод, немного уступая итерационному методу [1], на порядок превосходит метод установления по скорости сходимости.
Число итераций
]_I_I
200
300
400 500
Для получения точного решения в качестве правых частей использовались соотношения из [7]:
Рис. 1. Точное решение. Сходимость разных методов
На рис. 2 показано влияние параметра а на сходимость итерационного процесса в предложенном методе. Необходимо отметить, что при малых правых частях уравнений (1) можно использовать большое значение параметра а, однако при этом увеличивается вероятность появления осцилляций. Для ускорения сходимости можно первоначальное значение а изменять через несколько итераций. На рис. 2 номеру 1 соответствует а=0,06; 2 - а=0,08;
[0,1; у < 10
3 - а=0,1; 4 - а = < . В последнем слу-
[0,2; у > 10
чае видно заметное ускорение сходимости.
На рис. 3 приведен профиль сопла и полученное распределение числа Маха, которое для всех рассматриваемых методов сошлось к точному решению, при этом х,=1.
Рис. 2. Точное решение. Влияние со на сходимость предложенного метода
вязкость газа 5-10-5 Пас при Т0; весовая доля второй фазы 0,4; число Прандтля 0,7; теплоемкость вещества второй фазы 1420 Дж/(кгК); молекулярная масса смеси 30 кг/кмоль; показатель адиабаты газа 1,1; вторая фаза монодисперсная, состоящая из А12О3; диаметр частиц второй фазы 10-5 м.
Кроме этого в описанном выше алгоритме изменены пункты 3 и 7:
- пункт 3 - по формулам (3) вычисляются Д^и Щ};
- пункт 7 - выполнялись релаксации: Д^Д^+о^-ЯГ) и Я^ЯГ+о^-ЯГ), а о изменялось следующим образом:
[0,9; у < 5 [0,45; у > 5' Полученные результаты продемонстрировали сходимость к решению х*«0,137, аналогичную представленной на рис. 1. Процесс сходимости показан в таблице, где для разных номеров итерации ] приводятся разности х* на двух соседних итерациях.
Таблица. Сходимость разных методов. Значения аЬз(х*!г1>-х*!(>)
ф =
Номер итерации Предлагаемый метод Метод [1] Метод установления [3]
5 3,2-10-2 3,310-2 1,4-10-2
10 1,310-5 5,7-10-5 8,1-10-3
15 7,2-10-8 8,510-7 5,910-3
20 3,310-9 1,110-9 2,910-3
500 2,610-9 0 5,7-10-5
920 1,2-10-10 0 9,4-10-6
Рис. 3. Профиль сопла Лаваля и распределение числа Маха вдоль него
Сравнение для двухфазного течения
С использованием упомянутых методов были проведены расчеты двухфазного неравновесного течения в сопле Лаваля. Начальное поле параметров (пункт 2) находилось аналогично предыдущему разделу. Конденсация, испарение, коагуляция и дробление частиц второй фазы не рассматривались. Коэффициенты межфазного взаимодействия и особенности расчета параметров второй фазы приведены соответственно в [1] и [8]. Использовались следующие параметры: давление торможения 50.105Па; температура торможения Т=3000 К; динамическая
Видно, что если в методе установлении с использованием явной схемы Маккормака [3] до совпадения 6 цифр после запятой требовалось около 900 итераций, то в предлагаемом методе и в методе [1] было достаточно 15 итераций.
Заключение
Изложенный итерационный метод решения уравнений газовой динамики унаследовал простоту и скорость сходимости ранее предложенного метода, но в отличие от последнего использует уравнения в дивергентной форме. Это позволяет рекомендовать приводимый метод для расчета одномерных стационарных неравновесных газодинамических течений, имеющих единственный переход через скорость звука с монотонным возрастанием числа Маха вдоль рассматриваемой области.
СПИСОК ЛИТЕРАТУРЫ
1. Галкин В.М. Итерационный метод решения одномерных уравнений газовой динамики // Известия Томского политехнического университета. - 2002. - Т 305. - № 8. - С. 130-136.
2. Галкин В.М. О методе решения одномерных стационарных уравнений газовой динамики. // Математическое моделирование. - 2003. - Т 15. - № 11. - С. 30-36.
3. MacCormack R.W The effect of viscosity in hyperbolicity impact cratering // AIAA Paper. - 1969. - № 354. - 17 p.
4. Роуч П. Вычислительная гидродинамика. - М.: Мир, 1980. - 616 с.
5. Годунов С.К., Забродин А.В., Иванов М.Я., Крайко А.Н., Прокопов Г.П. Численное решение многомерных задач газовой динамики. - М.: Наука, 1976. - 400 с.
6. Лойцянский Л.Г. Механика жидкости и газа. - М.: Наука, 1987. - 840 с.
7. Галкин В.М. Пример точного решении и тестовые расчеты для одномерных стационарных уравнений газовой динамики // Математическое моделирование. - 2005. - Т 17. - № 1. - С. 3-9.
8. Глазунов А.А., Рычков А.Д. Исследование двухфазных неравновесных течений в осесимметричных соплах // Известия АН СССР Механика жидкости и газа. - 1977. - № 6. - С. 86-91.
