Ю. А. БУЛЫГИН, А. В. КРЕТИНИН В. С РАЧУК, С В. ФАЛЕЕВ
РАСЧЕТ ТЕПЛОВОГО СОСТОЯНИЯ КАМЕРЫ ЖРД
УЧЕБНОЕ ПОСОБИЕ
ВОРОНЕЖ 1997

МИНИСТЕРСТВО ОБЩЕГО И ПРОФЕССИОНАЛЬНОГО ОБРАЗОВАНИЯ РОССИЙСКОЙ ФЕДЕРАЦИИ
Воронежский государственный технический университет
Ю. А. Булыгин А. В.Кретинин В.С.Рачук С. В. Фалеев
РАСЧЕТ ТЕПЛОВОГО СОСТОЯНИЯ КАМЕРЫ
ЖИДКОСТНОГО РАКЕТНОГО ДВИГАТЕЛЯ
Учебное пособие
ВОРОНЕЖ 19S"
УДК 621.454.2.018(075.8)
Расчет теплового состояния камеры жидкостного ракетного двигателя: Учеб, пособие/ Ю. А. Булыгин, А.В.Кретинин, В. С. Рачук, С. В. Фалеев; Под ред. В. П. Козелкова; Воронеж, гос. техн. ун-т. Воронеж, 1997. 90с.
В учебном пособии приводятся основные принципы построения одномерной математической модели горения в камере двигателя, основы газодинамического профилирования по результатам термодинамического расчета, математическая модель наружного охлаждения камеры с изложением практического алгоритма расчета на ЭВМ.
Учебное пособие предназначено студентам специальности “Ракетные двигатели*’ для использования в учебном процессе в рамках базовых учебных курсов “Теория двигателей”, “Термодинамика и теплопередача”, “Математическое моделирование”, а также для курсового и дипломного проектирования.
Табл. 11. Ил. 9. Библиогр.: 18 назв.
Печатается по решению редакционно-издательского совета Воронежского государственного технического университета
Научный редактор заслуженный деятель науки и техники РФ, д-р техн, наук, профессор В.П. Козелков
Рецензенты: кафедра гидравлики и энергетики ВГЛА;
д-р. техн, наук, профессор В. Н. Мелькумов
(Г)Булыгин Ю. А. , Кретинин А. В. , Рачук В.С., Фалеев С.В., 1997
© Оформление. Воронежский государственный технический университет, 1997
ВВЕДЕНИЕ
Постоянное стремление конструкторов-разработников к увеличение значений удельных параметров Стяга, удельный импульс, время работы) при минимальных массе и габаритах двигателя весьма обостряет проблему правильной оценки и выбора основных принципов проектирования жидкостного ракетного двигателя СЖРД). Для конкретной схемы двигателя, системы подачи и охлаждения, термогазодинамических параметров в камере сгорания и выходном сечении сопла, значения коэффициента избытка окислителя, допустимого уровня потерь в камере сгорания и сопле проводится серия проектных расчетов, в результате которых определяются энергетические и геометрические характеристики камеры. Последующее детальное проектирование двигателя в целом, а также отдельных его узлов и систем, может привести к необходимости корректировки данных проектного расчета. Особое внимание следует уделять расчету такого тёплонапряженного узла, как камера сгорания ЖРД, которая во многом определяет работоспособность и надежность конструкции летательного аппарата СЛА).
Задача овладения современными методиками термодинамического и теплового расчета при проектировании камеры ЖРД студентами специальности "Ракетные двигатели" решается в настоящем учебном пособии.
1. ТЕРМОДкШМИЧЕСКШ РАСЧЕТ КАМЕРЫ ХИДКОСТНОГО РАКЕТНОГО ДВИГАТЕЛЯ
Для определения удельных параметров и геометрических размеров ракетного двигателя необходимо знать температуру, давление, состав и скорости продуктов сгорания в камере сгорания в различных сечениях сопла. Эти данные необходимы также и для расчета охлаждения стенок камеры.
Высокие температуры горения в камере ЖРД вызывают сильную диссоциацию продуктов сгорания; которая,в свою очередь, влияет на их температуру и состав. Поэтому при определении температуры и состава продуктов сгорания в ЖРД учет диссоциации является обязательным. Кроме того, следует принимать во внимание и зависимость теплоемкости продуктов сгорания от температуры.
Исходными данными для термодинамического расчета горения и истечения продуктов сгорания являются:
химический состав горючего и окислителя;
соотношение между количеством горючего и окислителя в камере сгорания;
давление в камере сгорания;
давление на срезе сопла.
Расчет температуры и состава продуктов сгорания на выходе из камеры сгорания (перед реактивным соплом) производится при следующих допущениях:
смешение компонентов топлива является полным;
физическая неполнота сгорания отсутствует;
нет отдачи тепла в стенки камеры;
процесс сгорания протекает при постоянном давлении;
состояние продуктов сгорания на выходе из камеры сгорания является полностью химически и энергетически равновесным.
При расчете процесса расширения в реактивном сопле ЖРД принимаются следующие допущения:
процесс расширения химически и энергетически предельно равновесен;
догорание топлива в сопле отсутствует;
нет отдачи тепла в стенки сопла;
трения и газодинамических потерь в сопле нет.
При химической и энергетической равновесности, отсутствии теплообмена с внешней средой и газодинамических потерь процесс расширения является изоэнтропным ввиду его ббратимости.’ Обоснование допущений при расчетах горения и истечения в камере двигателя, получение исходной системы уравнений для расчета равновесного состава и методы решения этой системы подробно освещены в литературе [1] - (3J.
9 настоящее время для широкого круга исследовательских и проектных работ термодинамический . расчет горения и расширения продуктов сгорания проводится при помощи справочника C4J. Так как основная форма представления справочных данных в нем табличная для ограниченного набора топливных композиций, то его использование на начальном этапе обучения не позволяет студентам провести осознанный многовариантный термодинамический расчет камеры. При последующем курсовом и дипломном проектировании использование указанного справочного издания представляется методически обоснованным и эффективным за счет значительного сокращения объема вычислений при выборе оптимальных параметров камеры двигателя при выбранной топливной композиции.
В представленном учебном пособии приводятся сведения пс определению исходных данных для термодинамического расчета, по расчету состава продуктов сгорания и их температуры, по расчету процесса расширения продуктов сгорания в реактивном сопле. При этом в качестве иллюстрации излагаемого материала приводятся примеры для топлива, состоящего из следующих элементов: водорода, углерода, кислорода и азота, что ни в коей мере не ограничивает общности подхода к задачам термодинамического расчета камеры ЖР2.
1.1.	Состав компонентов топлива
Состав вещества в массовых долях отдельных элементов называется элементарным составом. Массовую долю конкретного химического элемента будем обозначать его символом с индексом "г” для входящего в состав горючего, с индексом “ок“ для входящего в состав окислителя. Например, Н^. - массовая доля водорода в
6
горючем, Нок - массовая доля водорода в окислителе и т. д.
Сумма массовых долей всех элементов, имеющихся в данном компоненте топлива, должна быть равной единице.
Если горючее или окислитель представляют собой индивидуальное химическое соединение, то массовые доли входящих в них элементов находятся по следующей формуле:
(1.1)
где Bj, - массовая доля к-го элемента;
ак - число атомов данного элемента в молекуле химического соединения;
- атомная масса этого элемента;
М - молекулярная масса соединения.
Если ограничиться элементами Н , С , 0 , N , то в общем случае химическая формула вещества имеет вид:
CRON	(1.2)
с п о п
Тогда элементарный состав можно представить как
с	h	16 о	14
С = ------- , Н = ----- ,	0 = ------- , N = --------- , (1.3)
И	И
здесь p = 12c + h + 16o +14 п; С, Н, О, N - массовые доли углерода, водорода, кислорода и азота соответственно.
Если топливо или его компонент представляет собой комбинацию нескольких веществ, то массовая доля отдельного (к-го) элемента
Ви .	ч
где BR - массовая доля к-го элемента в смеси;
gt - массовая доля отдельного Ci-гоЗ вещества в смеси;
Вп - массовая доля к-го элемента в i-м веществе.
Для смеси из веществ, состоящих из элементов Н , С , 0 , N ,
элементарный состав определяется следующим образом:
С’ЬЛ  "-2S, в, •	°,  » Ч9, ", <1®
Здесь gt - массовая доля отдельного Сi-го) вещества в смеси;
С , Н. , 0*, Nt - массовые доли элементов в отдельном С1-м) веществе.
Если топливо состоит из окислителя и горючего и известны соотношения компонентов х и элементарный состав обоих компонентов, то массовая доля отдельного (к.-го) элемента в топливе
вк
1 + X
С1.6)
Для топлива, состоящего из элементов Н, С, 0, N , элементарный
состав топлива определяется следующими соотношениями:
С + ц С . Н + х Н и о + X 0 ь N + X N . г * ок.	г ок	г ок	г ок
С = ----------; Н = ----------; 0 = ----------; N = ----------- . С 1.7)
1 + X	1 + X	1 + X	1 + X
Если компоненты представляют собой смеси индивидуальных веществ, то в ряде случаев целесообразно использовать условную химическую формулу данного компонента. Такую формулу можно построить разными способами. Например, удобно определять ее, исходя из числа атомов различных элементов, приходящихся на м=100 для рассматриваемого компонента. В этом случае условная химическая формула будет иметь вид
Сс Hh Оо Nn ,
Cl.8)
100 С	100 Н	100 0	100 N
где с = -------- ; h = --------- ; о = ---------- ; п = -------- .
12	1	16	14
Здесь С, Н, О, N - массовые доли соответствующих элементов в данном компоненте топлива.
Если вычисления ведутся не на 1 кмоль, а на 1 кг топлива, то вместо молекулярной формулы компонента записывают удельную химическую формулу:
Сс НК Оо Nn ,
Cl. 93
с где с = — 12
h _ о
Н = --- ; о = ------ ;
1	16
п
14
п =
Напрмер, для керосина с элементарным составом, содержащим 86,7% углерода и 13,3% водорода, условная химическая формула имеет вид С Н з , а удельная формула для 1 кг керосина будет с7'гаэ Н О,07220	0,133
1.2. Соотношение компонентов топлива
Расчет состава топлива сводится к определению соотношений между горючими и окислительными элементами.
Эти соотношения могут быть найдены из представления реакции окисления как обмен электронами на внешней электронной оболочке атомов. При образовании продуктов полного окисления (сгорания) горючие элементы СГЭ) отдают все свои внешние электроны, а окислительные элементы (ОЭ) дополняют число внешних электронов до восьми. Число электронов, отдаваемых или приобретаемых элементом при химической реакции ,* определяет валентность СД) Для горючих элементов Д = п , а для окислительных Док = 8 - п , где п - число электронов на внешней оболочке. В приложении 1 приводится число электронов п для ряда элементов.
Используя указанные определения, можно найти продукты полного окисления для различных комбинаций иэ окислительных и горючих
g
элементов, полагая Др = Док. Некоторые иэ них даны в приложении 1.
Для топлив, состоящих из элементов Н , С , 0 , N , продуктами полного окисления будут Н20 и СО* . Азот окисляется кислородом только при высоких температурах, при этом реакция его окисления эндотермическая. 0 продуктах полного окисления азот присутствует в молекулярном состоянии .
Пусть число свободных валентностей в ГЭ в 1 кг горючего равно:
Дг - I АОк Г £ г ------------,	C1.1OD
^ог
где У Д^ , у Д°к - сумма валентностей ГЭ и ОЭ в молекуле горючего соответственно.
Аналогично число свободных валентностей ОЭ в 1 кг окислителя равно:
У Док - У Д\
L ок 4 ок ------------— , (1.11)
^ок
где у Д^к , £ Д°£ - сумма валентностей ГЭ и ОЭ в молекуле окислителя соответственно.
Тогда теоретическое (стехиометическое) соотношение компонентов топлива:
1^-1
С1.12}
5	- у дг.
4 ok	ок
Величина для ряда топлив дана в приложении 1.
В тех случаях, когда число возможных элементов в топливе ограничено, целесообразно для определения величины хо пользоваться формулой, основанной на массовых долях элементов;
8 --------- С + 8 Н - 0 3 г г г у =------.	(1.13)
°	8
0 v-------с - 8 Н v
оК. з г	ок
Отношение количества окислителя % , в действительности приходящегося на один килограмм горючего, к количеству окислителя Хо , теоретически необходимому для полного окисления одного килограмма горючего, называется коэффициентом избытка окислителя и равно;
а = — (1.14)
При а < 1 топливо называется “богатым" по горючему с недостатком окислителя. Если же а > 1 , то топливо называется “бедным" по горючему (избыток окислителя). Если а = 1 , то топливо часто называют "стехиометрическим".
1.3. Состав продуктов сгорания
Состав продуктов сгорания определяется исходным составом топлива и зависит от температуры и давления, а также от полноты сгорания. Состав продуктов сгорания в камере ЖРД для четырех групп топлива приводится в приложении 1.
Термодинамический расчет ведется в предположении, что состав продуктов сгорания является равновесным при данных температуре и давлении. В состоянии химического равновесия количества отдельных газов в реагирующей смеси находятся в строго определенном соотношении, исходя из уравнения химического равновесия (уравнения константы равновесия). Для реакции;
аА + вВ = сС + dD ,	С1.15)
где А, В, С, D - регулирующее вещества;
а, в, с, d - стехиометрические коэффициенты
Уравнение химического равновесия имеет вид
pc pd C D К = --------- ,	(1.15)
Р ра рв А В
где Кр - константа равновесия;
Р , Рв, Рс, PD- парциальные давления соответствующих веществ.
Количество независимых уравнений, описывающих равновесие реагирующей смеси, равно:
Z = X - Y ,	(1.17)
где X - число отдельных газов, входящих в смесь;
Y - число химических элементов, составляющих эти газы.
Если в состав продуктов сгорания входят газы Н О, СО , СО, Н , 02, Na, ОН, NO, О, Н, N, то X = И, Y = 4, Z = 7 .
Расчет состава продуктов сгорания заключается в определении количеств отдельных газов в смеси продуктов сгорания. Число неизвестных равно числу возможных газов X .
Для решения этой задачи используются уравнения химического равновесия и материального баланса.
Уравнения материального баланса выражают закон сохранения материи и определяют равенство количества отдельных элементов в исходном топливе и в конечных продуктах сгорания. В общем случае уравнения материального баланса имеют вид
В. = А. V av< М4 ,	* 1.18)
к к 4 к1 i
где Вк - массовая доля к-го элемента топлива, определяемая по формуле С1.6); Ак - атомная масса к-го элемента;
afci - число атомов к-го элемента в i-м компоненте смеси
продуктов сгорания;
Мх - число молей i-го элемента.
Система 2 уравнений химического равновесия и уравнений материального баланса, согласно (1.17), включает X уравнений и достаточна для решения задачи.
Состав газообразных продуктов сгорания может быть выражен в молях Мх отдельных газов на 1 кг продуктов сгорания, в массовых долях g или в парциальных давлениях Р, .
Если состав смеси газов задан парциальными давлениями, то ее кажущаяся молекулярная масса определяется по формуле
Р! 2 Р1
(1.19)
соединение
Здесь - молекулярная масса i-го газа.
Для гетерогенной системы, в которой одно и то же присутствует как в газообразном, так и в конденсированном виде, для определения состава продуктов сгорания добавляется еще уравнение зависимости давления насыщенного пара Ps данного соединения от температуры:
Р = f ( Т ) .	(1.20)
S
Введение уравнения (1.20) обусловлено появлением дополнительного неизвестного - доля вещества в конденсированном виде.
1.4. Исходная система уравнений для расчета равновесного состава
Исходная система уравнений, необходимых для определения количеств отдельных составляющих смеси продуктов сгорания включает в себя уравнения химического равновесия и материального баланса. Уравнения химического равновесия описывают реакции диссоциации. В общем случае система уравнений диссоциации может быть составлена произвольно, но так, чтобы она полностью описывала состав смеси продуктов сгорания и химическое взаимодействие между отдельными составляющими смеси. В ряде случаев при этом исходят из уравнений диссоциаций молекулярных соединений на атомы.
Напрмер, для топлива состоящего уравнений диссоциации имеет
из элементов
О и Н система
вид
Н О
= Н +
1
----О
2 1
н о
= ОН +
1 -------Н
2 '
С1.21)
О
20 ,
н
2Н .
Вместо этой системы может система химических уравнений атомы:
быть диссоциации молекул
записана и
ей
эквивалентная и радикалов на
На0 = 2Н + О/
ОН = 0 + Н ,
0 = 20 , 2
Нг = 2Н .
С1.22)
Для системы (1.21) уравнения химического равновесия имеет вид
К
р
р р° >8 н О
2	2
РНгО
	К = D * 2	
		'он н
			2
		
		Р-л Л
		я^о
		
		
	К = 	4 О
		О
	к = •	“ О 2 4 И
а для системы (1.22)		“ н 2 ра р Н 0
	К ’ = р 2	Р
		л Р Р о гн
	К =	р
	К =	Л он _Д_ Ро
	к =	2 ‘ н
		н 2
(1.23)
(1.24)
Схема разложения молекул на атомы (1.21) и С 1.22) оказывается не очень удобное в тех случаях, когда в продуктах сгорания отсутствует какой-либо из элементов в атомарном состоянии, хотя он и входит в состав других соединений. Например, в продуктах сгорания топлива, состоящих из элементов Н, С, О, N , практически отсутствует свободный углерод, входящий в состав поутих газов, как
СОа и CO .
Константа равновесия для идеальных газов зависят от температуры и приводятся в приложении 5 методических указаний £53.
Если в таблицах отсутствуют константа равновесия.
соответствующие принятой схеме реакций, то они могут быть получены комбинацией имеющихся в таблицах. Например:
рн ро'в рн *н ро • ро °'8 к» ’ MU	М М v f w	Р
к = -----2---2_ = ---— ------- Г ---5- ]	= -----_*--- ;
₽l Рн о ₽н ₽н о *0	К К0’8
2	Н	2	Р4 Р3
£1.25)
Рон Рн'“ Рон *н РО РН °’’ К ’ им м	им Н v 4* Н	р
К = ---------— =--------------[ ----2_ ]	= -----□---- .
3 Рн о Рн Ро Рн 9 Рн	КК®'8
2	М О 2	Р^ р^
Рассмотрим составление исходной системы уравнений для расчета равновесного состава диссоциированных продуктов сгорания на примере топлива, состоящего из элементов 0 и Н .
Ожидаемая температура сгорания такого топлива при а < 1 согласно (53 (приложение 5) составляет Тк > 3000 К . При таких температурах в состав продуктов сгорания входят газы На0, 0а, Н , ОН, О, Н (см. приложение 4). В данном случае исходная система уравнений включает в себя Z=X-Y =6-2=4 уравнения химического равновесия и два уравнения материального баланса. Уравнения химического равновесия представлены (1.23) и (1.24) . Уравнения материального баланса на основании (1.18) могут быть записаны в следующем виде:
во = 16 £ мнао + 2 «о/ «ОН + «О 3 ’
(1.26)
Вн = 1 £ 2 Мн,. + 2 М„а+ Мон + И,, ) или в парциальных давлениях:
В = 16 С P„ n + 2 Pn «• PftU + P„ ) -22- ; О	НО	О ОН О P
‘	a	™	Cl. 27)
BH = 1 C 2 P + 2 PH ♦ Po„ ♦ P„ )	.
a	•	PCM
Для исключения сомножителя Мсм / Р^ целесообразно перейти к относительным уравнениям материального баланса:
16 ( Ри + 2 Рл + Pnu + РЛ )
В	Н2°	0Н 0
—2-= ----------------------------------- . С 1.28)
Вн 2 Рн о + 2 Рн 4 Рон + Рн а	а
Так как число уравнений при этом уменьшается, то исходная сиситема должна быть дополнена еще одним уравнением - уравнением является уравнение баланса парциальных давлений (закон Дальтона):
РСМ = 2 Р1
(1.29)
Для данного топлива (1.29) имеет вид:
рсн = рно *ро ♦ РН * РОЯ * РО ‘ Р„  С1-ЗЮ а а а
Таким образом, для определения шести неизвестных парциальных давлений имеются шесть уравнений, с помощью которых может быть найден состав продуктов сгорания при известных элементарном составе исходных веществ, давлении и температуре продуктов сгорания.
Система уравнений имеет только одно единственное решение, состоящее из системы корней только положительного и действительного значения. Следовательно, задача сводится к отысканию именно этой единственной системы корней полученной системы уравнений, но может решаться различными методами.
например, с использованием программы [61.
1.5.	Полная энтальпия топлива и продуктов сгорания
Общий запас энергии в топливе или продуктах сгорания оценивается величиной полной энтальпии J , равной сумме физической энтальпии I и химической энергии 0х , т.е.
J = I ♦ Qx .
(1.31)
Величина I, входящая в С 1.31), представляет собой то количестве теплоты, которое затрачивается на нагревание данного вещества ст некоторой начальной температуры То до той температуры Г . при которой оно используется в двигателе . За химическую энергию 1х принимают тепловой эффект реакции образования данного вещества из исходных.
Значения полной энтальпии компонентов топлива в кЛх/кг при стандартных условиях приведены в приложении 7 источника i5j.
Полная энтальпия топлива определяется полными энтальпиями горючего и окислителя и количественным соотношением последних. Полная энтальпия одного килограмма топлива раздельной подачи при соотношении компонентов % равна
Полная энтальпия продуктов сгорания одного килограмма топлива может быть определена по формуле
где Р. - парциальное давление i-rc газа s поодуктах огспания.
u - молекулярная масса i-го газа в продуктах сгорания.
Значения полной энтальпии i-го газа берутся из приложения 8 источника £5j в ккал/кмоль и переводятся в кДж/кмоль по соотношению
[кДж/кмоль] = 4,187 [ккал/кмоль] .
(1.34)
Как показано в [7], для ЖРД с различными системами питания при принятых допущениях полная энтальпия продуктов сгорания при температуре ?к равна энтальпии исходного топлива в баке при начальной температуре.
1.6.	Определение параметров продуктов сгорания на входе в сопло
Исходным уравнением при определении температуры сгорания в ракетных двигателях является уравнение теплового баланса:
В уравнение теплового баланса входят все искомые величины термодинамического расчета камеры ЖРД, т. е. температура Тк и все парциальные давления Р отдельных газов, составляющих продукты сгорания.
Определение температуры сгорания следует вести одновременно с расчетом состава продуктов сгорания в камере ЖРД. Последовательность расчета следующая.
Для заданных компонентов топлива определяется элементарный состав по формулам (1.1) или (1.3), стехиометрическое соотношение по формулам (1.12) или (1.13), действительное соотношение из формулы (1.14) и по формуле (1.32) полная энтальпия топлива. По приложению 5 находится ожидаемая температура продуктов сгорания Тк1 и задаются еще два значения этой температуры Тк2 и Ткз, лежащими в районе Тк1 и отличающимися друг от друга на 100-150 К.
Рис. 1.1.
По известному составу продув смеси:
Для каждой температуры определяется состав продуктов сгорания в соответствии с п.1.4 -.
По формуле (1.333 для температур Tfc, . Tka . вычисляются значения энтальпии продуктов сгорания J , J , Jj3 и строится график зависимости энтальпии продуктов сгорания от температуры Срис. 1.1}. Пользуясь соотношением (1.35} и графиком на рис. 1.1 определи-. ется теоретическая температура продуктов сгорания на входе в сопло Т. . к
сопло Тк .я находят энтрслию
1
Sv = ------------5 Р< С S, - 1,986 Ln Р, ) ,
к	р Z. 1	1	1
и к
с 1.36:
где - абсолютное значение энтропии i-го газа при стандартном давлении и данной температуре, значения которой в ккал/кмоль»К приведены в приложении 9 источника [5].
По трем расчетным точкам строят график зависимости энтропии продуктов сгорания от температуры Срис. 1.2} и определяют энтропию
Затем строят графики зависимости парциальных давлений i-тых газов Рх , входящих в смесь, от температуры и определяется состав продуктов сгорания для Т, Графическое определение Р, для температуры проводится аналогично графику на рис. 1.2.
го
Лс известным Р. определяют кажущуюся молекулярную массу омеги продуктов сгорания по формуле (1.19) и газовую постоянную:
JL = ----- (1.37)
1.7.	Расчет параметров продуктов сгорания на выходе из сопла
Свойства, которые заданы в п.1.1 процессам горения и расширения, позволяют знать два термодинамических параметра состояния в точке завершения процесса. Так, для моделей процессов расширения в реактивном сопле ЖРД один из параметров состояния -энтропия - известен из расчета точки завершения предыдущего процесса - горения, другой - задается непосредственно - давление на срезе сопла. Два термодинамических параметра в условиях равновесия полностью определяют термодинамическое состояние любой СИСТвМЫ.
Расчет изоэнтропического расширения и истечения продуктов сгорания в реактивном сопле ЖРД производится в следующем порядке.
Для заданных компонентов топлива по приложению 6 определяется средний показатель изоэнтропического процесса к расширения продуктов сгорания и вычисляется величина ожидаемой температуры в выходном сечении сопла:
Т = Т. ( —— 1 k	(1.38)
1 Рк
В области ожидаемой температуры продуктов сгорания Та1 выбирается два значения температур Таа и Таэ , отличающихся от Та1 на 100-150 К. Для каждого значения температуры определяется состав продуктов сгорания. По формуле (1.36) для каждого значения температур Т , Таа , Таа вычисляются значения энтропии продуктов сгорания Sai , Sa> , SaJ и строится график зависимости энтропии
смеси от температуры на срезе сопла Срис. 1.31.
Используя условие изоэнтропично-сти течения продуктов сгорания в реактивном сопле Sk = Sa и график на рис. 1.3, определяется теоретическая температура продуктов сгорания на выходе из сопла Та .
Рис. 1.3
Для найденного значения температуры Тк определяется величины парциальных давлений газов Р£ составлявших продукты сгорания, аналогично п. 1.6.
По известным парциальным давлениям определяется энтальпия смеси J по формуле С1.33), кажущаяся молекулярная масса ма смеси продуктов сгорания по формуле С 1.19) и газовая постоянная Ra этой смеси по формуле С 1.37) .
22 II. РАСЧЕТ ОСНОВНЫХ ПАРАМЕТРОВ КАМЕРЫ ДВИГАТЕЛЯ
Исходными данными для расчета основных параметров камеры жидкостного ракетного двигателя, расчетная схема которой приведена на рис. 2.1, являются результаты термодинамического расчета параметров продуктов сгорания в камере сгорания С сечение к-к) и на выходе из сопла (сечение а-аЗ. Термодинамический расчет производится с помощью главы I.
Рис. 2.1. Расчетная схема камеры сгорания ЖРД
В сечениях к-к и а-а по данным термодинамического расчета известны температура Т, полная энтальпия J, газовая постоянная R, молекулярная масса ц продуктов сгорания.
2.1. Определение теоретических параметров камеры
Теоретические параметры камеры сгорания определяются в предположении ее изобаричности, изоэнтропичности процесса расширения продуктов сгорания в сопле и равновесности процессов в камере двигателя. Теоретические параметры в дальнейшем обозначаются со штрихом сверху, например Va и т.д.
Скорость истечения продуктов сгорания на выходе из сопла
V' = 44,72 J I - J n	к а
(2.13
Плотность продуктов сгорания в выходном сечении сопла
(2.23
Средний показатель равновесного расширения продуктов сгорания
Параметры продуктов сгорания в критическом сечении сопла
п + 1
(2.43
Ркр п + 1
(2.53
< = 1 2 ---------D----R T.
kp	n + । к к
(2.7)
Значение газодинамической функции полного импульса потока
где \ = V' / а. - коэффициент скорости потока на срезе сопла. А А )Ср
Значение газодинамической функции расхода п	п
[П + 1	п-1	, г П - 1 ,а 1 п-1
--------! Ч I 1------------------Ч I С2-93 2 J а 1 п + 1 а J
Удельный импульс тяги на расчетном режиме С Р* = Рн )
(2.10)
где В
гу
R
Удельный импульс тяги в пустоте ( Рн = 0 )
В
2	п-1
п + 1 3
Z С X а
)
(2.11)
Секундный расход топлива через камеру двигателя
Р
ш_ = —~ ( 1 - т_. ) , Т j	г г
(2.12)
где тгг = тгг / тт “ относительный расход генераторного газа (для схем с дожиганием генераторного газа mrr = 0 , для схем без дожигания величина mrr определяется по методике, изложенной в [3J и [41). Тяга Р является заданной. Для двигателей первых ступеней ракет обычно задается тяга у Земли (Рн_о) и подсчет ведется по 1уН-о . определяемой формулой (2.10) при Ра = Рн_о Для двигателей верхних ступеней задается тяга в пустоте Рп и расход rrtj. определяется по величине I , вычисленной по формуле (2.11).
Расходный комплекс
13 =--------—
V
(2.13)
n + 1
2
где
V =
п + 1
Площадь критического сечения сопла и его диаметр
р- 13 fcp р к
(2.14)
Чр
4 _____La.
я
(2.15)
Площадь выходного
сечения сопла и его диаметр
F" = Л ,	C2.163
a P v a
4 F d = Л ------2-	(2.17)
a	я
Геометрическая степень расширения сопла
2.2. Определение основных параметров камеры
Реальные процессы в камере двигателя происходят, с заметным отличием от принятых схем полного адиабатного и изобарного горения и изоэнтропического однородного и равновесного течения в сопле. Эти отличия вызывают уменьшение удельного импульса тяги.
Действительный удельный импульс тяги
I = I' л р . (2.19) у у гс
В формуле (2.19) рк - коэффициент, учитывающий потери удельного импульса из-за несовершенства рабочего процесса непосредственно в камере сгорания. Величина рк зависит от качества организации процессов смесеобразования в камере и колеблется в пределах 0,96 -0,99; для камер большой тяги рк = 0,98 - 0,99 .
Входящий в формулу (2.19) коэффициент рс учитывает потери удельного импульса из-за несовершенства процесса истечения газа через сопло. Его величина может быть определена по зависимости
₽с = 1 - ?р - ?тр -
С2?20)
где f - потери удельного импульса из-за рассеивания; f - потери в сопле из-за трения; ?н - потери в сопле из-за химической неравновесности; ?з - потери в сопле из-за многофазное™. Значения-величин ?р , £ , ?н могут быть рассчитаны по следующим соот ношениям:
г	R П
£ =	13,2 - (89п - 98) , -2-	(1 - Z2) 1О'г ; (2.21;
Р I	* D J
кр
= ?’ Г 1 - ( 1,19 + sTp Тр [ I
0,014 п - 1
е-(: -z>
(2.221
[(-
0,28 + Т'0,эз W_____________
1,06п - 0,86
1 + 1,88-1 п-1 - 1075
R -R	R J
max kp 0>8 _ (0г018 п - 0,0175) С—2- + 24
R.
<Р
Здесь hQ - приведенная длина релаксационной зоны;
Pq = 1,0332x10э Па ; R - радиус сечения сопла, в котором
/3 = 0 ; (3 - угол непараллельное™ истечения на срезе сопла. Значения показателей а , b и длины релаксационной зоны h
приведены для ряда топлив в табл. 1.
Таблица 2.1
Компоненты топлива	 । 	 ! а и ।	ак	t	,	« |	ho,	м	| 1	»	1 t 2 !	b
0 > Н	3.8	1.3 1С-3	С. 4	*
2.
О ♦ Т - 1 2
АТ + НДМГ
Продолжение таблицы 2.1
0,8	6,7 10"4	0,35	0,8
0,9	5,3 1О'в	0,25	0,5
Потери в сопле из-за многофазности продуктов сгорания зависят от многих факторов: массовых долей жидкоЗй и твердой фазы, дисперсности, формы частиц и распределения их по сечению камеры и др. Для достаточно надежного определения коэффициента требуются дополнительные исследования.
Действительный расход топлива, необходимый для
обеспечения заданной тяги
р	р
— ( 1 - m 1 = —---------------- ( 1 - го ) .	(2.26)
I 1 гр J I р	гг
у	у *с
Расход окислителя через камеру сгорания
К m ср
(2.27)
сечению
где К - среднее массовое соотношение компонентов по ® ср
камеры сгорания.
Расход горючего через камеру сгорания
Шг = Шу - mok	(2.28)
Расход топлива через ядерные форсунки
оц, =	( 1 - шпр) ,	(2.29)
где швр = швр / Шу - относительный массовый расход топлива через пристеночные форсунки, равный 0,05-0,20 При этом в камерах
29 малой тяги доля расхода через пристеночный слой будет большей, чем в камерах с большой тягой.
Для ориентировочных расчетов расход величины может быть найден из эмпирического соотношения
mn = С - 1,79 + 1,84 10‘2 J Р ) пр
Расход топлива через пристеночные форсунки
^пр
Расход компонентов через ядерные форсунки
; V = тТ. ток, ’ Я	Я	я
(2. 32:
я
где К = а К m я я то
массовое соотношение компонентов в ядре камеры
сгорания;	Кто - стехиометрическое соотношение компонентов;
- коэффициент избытка окислителя, равный для углеводородных
топлив 0,75-0,95.
Расход компонентов через
пристеночные форсунки
К
_	m д Р -I	m - m	m
°knp’ 1 ♦ К	V V oknp’
mnp
где К = а К массовое отношение компонентов в пристеночном ** тпр пр то	г
слое; апр ” коэффициент избытка окислителя в пристенсчнсм слое, равный 0,30-0,45.
Действительная площадь критического сечения сопла

:2.34)
30 здесь - коэффициент восстановления давления торможения на входе в сопло, определяемый из уравнения
Р. г	П - 1	— -|
а, = —=	( 1 + X? ) ( 1-----------X?	, (2.35)
• Р L k	п + 1	J
1
где имеем п - средний показатель изоэнтропы расширения в интервале от до Р^р , определяемый по справочнику [4] для критического сечения; Л.к - коэффициент скорости на входе в сопло, определяемый по графику рис. 2.2. При этом относительная площадь камеры сгорания Fk=FkzFkp находится в пределах 2-6. Рекомендации по выбору Fk приводятся в C2J , [9] .
Коэффициент снижения удельного импульса вследствие неизобаричности камеры сгорания
где п - средний показатель изоэнтропы расширения в интервале от Рко До Ра , определяемый по справочнику (4] при условии течения продуктов сгорания в ядре потока.
Действительная геометрическая степень расширения сопла
Площадь выходного сечения сопла
F = F F. . а а кр
(2.37)
(2.38)
Площадь сечения цилиндрической части камеры сгорания
F
к
= Fu к кр
Диаметры характерных сечений камеры двигателя
(2.40)
Рис. 2.2.График зависимости от Ffc
2.3. ПОСТРОЕНИЕ ГАЗОДИНАМИЧЕСКОГО ПРОФИЛЯ КАМЕРЫ ДВИГАТЕЛЯ
Газодинамический профиль образуется из состыкованных контура сопла и профиля камеры сгорания.
Для построения профиля цилиндрической камеры сгорания с плоской головкой необходимо определить объем и диаметр цилиндрической части камеры, форму и длину докритической части сопла.
Объем камеры сгорания
V, = х ,	(2.41)
kc kp пр
где Lnp - приведенная длина камеры сгорания, значения которой для различных топлив приведены в табл. 2.
Таблица 2.2
Топливо RP-1+0 С Н N +HNO Н +0 С Н 0Н+0 NH +F CH NO 2	2 2 2	3	2	2	22	2	32	3
Продолжение табл. 2.2
, м 1,5-2,5	1,5-2,0	0,25-0,5 2,5-3,0 1,0-1,5	7,5
(2.43)
(2.44)
Для ориентировочных расчетов величина приведенной длины может быть определена с помощью соотношения
L = 1,5x10’ / J 10 Рь ,	(2.42)
пр	к
здесь Рк - давление в камере сгорания.
Поскольку за объем камеры сгорания Vkc принимается объем от днища до критического сечения сопла
V = V + V SC ВСЦ С докр * то длина цилиндрической части камеры V - V вс с довр L" ° ' где Vc - объем докритической части сопла. Объем входной части сопла зависит от ее формы, которая может быть разной, но она должна обеспечивать безотравное течение продуктов сгорания и равномерное поле скоростей в критическом сечении.
Коническая сужающаяся конечная часть:
v<• vс F* ’	* F«» > • (г-4®
здесь Ьдожр - длина докритической части сопла, определяемая по формуле
где р ~ угол раствора докритической части сопла, равный 60+90° С рис. 2.3).
Рис. 2.3.Схема построения докритической части сопла
Для снижения потерь проводится плавное сопряжение входной части с цилиндрической частью камеры сгорания и юбкой сопла. Радиус скругления профиля критической области Rt = С 0,7 - 1,0 ) D . а радиус сопряжения на входе Ra = С 0,1 + 1,0 ) Rk . Радиусы дуг окружностей Rt и Ra выбирают исходя из условий сопряжения так, чтобы сопло имело плавные очертания.
Профилирование докритической части сопла обычно проводится на основании опытных данных, например, по двум сопряженным радиусам Rt = Dip и Ra (рис. 2.4). При этом с повышением Рк радиус кг следует брать большим, так как при меньшей кривизне контура входной части сопла пристеночный слой более устойчив. Величину ₽г рекомендуется определять из соотношения [21
R
R = -г- = 0,23 1О'Э Ри	(2.47)
Если Р. < 4 МПа. то R = 1, если Р. > 30 МПа. то R = 5. к	к
При принятой форме входной части сопла
V =	У^Н/СЗЬ Я + Гсу2* у + 4Dh/C6L )1|, (2.46)
сдокр кр ьх [ к J	ьх J L 7	7	эх J
где L = 0.5D, 4(2 - кГЁ” )г - [(R - 1) ГЁ~ - з].	С2.43)
s х	к г	-«к	I	"• < J
Рис. 2.4.Схема построения камеры и сопла ЖРД
Контур расширяющейся сверхзвуковой части (юбки) сопла может быть построен приближенным методом, основанным на результатах решения вариационной задачи о нахождении оптимального профиля сопла. При этом закритическая часть сопла приближенно аппроксимируется параболой.
Длина сверхзвуковой части сопла
Ч =	-
С2.51)
здесь Ц - безразмерная длина сверхзвуквой части сопла, определяемая из выражения
= Ео | 0,32 + 0,66 ехр[ а. С 1 - z .	(2.52J
х . 0004*0.1 h
где L = 3.16 R	; а = 2.73 + 0,9 ( п - 1
о	а	1
К = 0,7 - 0,33 ( Fa )°'ae;

Порядок построения контура сверхзвуковой части сопла:
1. Вычерчиваем контур сопла в области критического сечения по дуге окружности А-Ап (рис.2.4) с радиусом R( = 0,45 Rfcp .
2. По известным безразмерным значениям длины сопла Га и
Рис. 2.5.График зависимости R от La
3. Зная L , Ra , /Зт , /За найдем точку А , проводя под уголом (3 касательную к дуге А-А , и точку С по известным L и R а а
4. Для построения параболы провыодим из точки С под углом прямую Cf до пересечения с касательной Anf Разбив отрезки Anf и CF на 5-7 частей и соединив соответствующее точки 1.2. 3 и т.д. прямыми, строим огибающую параболу АпС . Линия АА^С и будет искомым контуром сверхзвуковой части сопла.
РАСЧЕТ НАРУЖНОГО ОХЛАЖДЕНИЯ
3.1.	Расчет конвективных тепловых потоков
Основные расчетные соотношения
Конвективный тепловой поток вдоль стенки камеры ЖРД вычисляется на основе решения уравнений турбулентного пограничного слоя. Непосредственно около стенки образуется пограничный слой, в котором скорость резко изменяется в поперечном направлении от скорости в ядре потока до полного торможения на стенке. Математические модели турбулентных пограничных слоев очень сложны, поэтому в инженерных расчетах используют приближенные полуэмпирические методики. Тепловой поток представляется точным соотношением
g= a_pviy-CJ„-J-_3,	(3.13
Т гх ог СТ где - безразмерный коэффициент теплоотдачи;
Рх - характерная плотность потока, переменная вдоль контура;
- скорость в ядре потока;
Jor- энтальпия замороженного потока при температуре торможения;
JCT- энтальпия потока при температуре стенки.
При наличии внутреннего охлаждения предполагается, что пограничный слой полностью находится внутри пристеночного слоя.
Для расчета по формуле СЗ.1) нужно знать коэффициент теплоотдачи о^, с вычислением которого связаны основные сложности расчета конвективного теплообмена. В. М. Иевлевым был предложен метод, основанный на применении некоторых априорных зависимостей для простых случаев течения несжимаемой жидкости в цилиндрическом канале вида:
^(г.^.РгЗ,	(3.23
где г, гт~ безразмерные параметры, характеризующие состояние динамического и теплового пограничных слоев в сечении;
С₽м
Рг=----- число Прандтля;
X
С - теплоемкость продуктов сгорания СПС);
м - вязкость ПС;
X - теплопроводность ПС.
Зависимость (3.2) при ее распространении на случай сжимаемого газа можно оставить практически без изменений, если рх вычислять по специальному соотношению
рх=----------------------—:--------------,	с 3.3)
г1+Тст	j. 3+ТСТ 9/32 °-*в
L г 4-1	1-4	16 -I
где рог- плотность ПС в ядре потока;
‘	относительная скорость потока;
/ &
W = / 2 ------ (ЯТЛ_) - максимальная скорость потока;
пах
Тст=Тст/Т0Г- относительная температура стенки;
?ст- температура стенки со стороны газа;
Тог- температура ПС в ядре потока;
Я- газовая постоянная ПС.
Тогда закон теплообмена для сжимаемого потока может приближенно аппроксимироваться зависимостью
0^=0.01352 z'0-07ezT'0’07ePr-*0-98.	(3>4)
Анализ показывает, что отношение параметров пограничного слоя zT/z можно представить в виде
т г	х ст	-|
— =	1.769 ---------з--------------------- .	(3.5)
2	l-T^O.1/32
Величина zT/z изменяется вдоль камеры сравнительно слабо, и при приближенных расчетах ее принимают равной этому отношению в критическом сечении, т.е. при
/3=/ Е7Г и выбранном Тст.
Параметр zT, характеризующий пограничный слой в данном сечении, определяется решением интегрального соотношения энергии. Это решение при TCT=const и zT/z=const дает следующую формулу:
С х
гт= -Г о
(3.6)
Здесь Reo~ характерное число Рейнольдса:
^ОГ^тах Аср
(3.7)
где рк- давление в камере;
Мог~ вязкость ПС при температуре в ядре потока;
С- величина, постоянная при k=const и TCT=const:
С- 4.42
(1+Т„_)1 •в2(3+Т„)°-1в
(3.8)
1-/3“ f> = ----------------------------------------,	(3.9)
’ „...Г. «• 2(1+T„)J L 4(3+T) vl	U1
где D= O/djtp" относительный диаметр камеры;
D- внутренний диаметр камеры;
с^р- диаметр критического сечения; х- относительная координата вдоль контура сопла.
Вычисления по формулам (3.3-3.9) дают наиболее точные
результаты, нс требуют огромных вычислительных затрат, занимающее даже при применении ЭВМ существенное время. Поэтому иа практике часто пользуются упрощенными методиками.
Первая приближенная формула
Б формулу СЗ.4) подставим выражения (3.5) и СЗ.б) для комплексов z и упростим. Используя одновременно формулу для рх, окончательно получим выражение для конвективного теплового потока: о.ев
-рк S ср В--------------------------------[кВт/м2), СЗ.10)
Г Jp dx] D* вРг°-ов
L 0 J	Р
Коэффициент В учитывает переход от вязкости рог к вязкости М1000 ПРИ температуре 1000 К, входящей в комплекс S, зависит от Т„_ и слабо- от показателя адиабаты k. В интервале Т =<0.1.......
0.8) с достаточной точностью приближается полиномом
В= С8. 3605+0.46072-Т„+0. 34524 Т* ) -ю’*	(3.11)
С1	G1
и изменяется в пределах С8.4+9.0)-IO'3 1/К0,10в.
1-/3*
р = -------------------------------------.	С 3.12)
/З2	9/32
I1 2C1+Tct)J L1 4C3+TOT)J
0.13 О.429 ^ог'^ст^ ’^1 000 ОГ S=----------------------------------------СЗ. 13)
О . 429	О. 300	0,19
В СТ0Г+Тст) СЗТ0Г+Тст) функция теплофизических параметров газа, которая зависит от рода топлива, соотношения компонентов в пристеночном слое и температуры стенки Тст. Аналитические выражения для комплекса S очень сложны и плохо приближают его, поэтому даже при расчете на ЭВМ приходится пользоваться интерполяцией по набору затабулированных полуэмпирических данных либо непосредственно рассчитывать энтальпии.
Вторая приближенная формула
С достаточной точностью можно положить в СЗ.10) ♦	2
--------------% 1-/3 .
гх	Л"
I ^3^1 L о J
Тогда формула (3.10). принимает вид
2	0.8В
1-13 р„	S
з= в- —• ------------— • ---------- [кВт/м2],
СЗ.14)
СЗ.15)
где 1-/3 =тСк) - газодинамическая функция.
В обеих приближенных формулах распределение тепловых потоков 1.82
вдоль сопла, в основном, определяется множителем 1/D
3.2. Расчет лучистого теплового потока
Основными источниками излучения являются трех- и большей атомности газы. Излучение одно- и двухатомных газов мало по сравнению с ними и обычно при расчетах не учитывается. На практике вполне достаточная точность получается при расчете излучения, соответствующего среднему по составу и температуре газу. Величина лучистого теплового потока определяется законом Стефана-Больцмана:
Т2 4 Т
q = с -С [(—)*- С—)4],	С3.16)
л п 0 L100J 1(ХГ J
где Г*- температура излучающего тела, К;
Г - температура поглощающего тела, К;
сп~ приведенная степень черноты системы;
коэффициент излучения абсолютно черного тела.
При расчетах обычно учитывают только наиболее излучающие гаэы-водяной пар и углекислый газ, для которых известны эмпирические
излучение
(3.17)
черноты стенки.
(3. 18)
зависимости:
eCpl.D;
£н o= eCp,pl,T\ 2
где p- плотность газа;
I- геометрическая длина луча;
Г- температура газа (Г в формуле (3.16)).
Кроме того, пренебрегают температурой стенки по сравнению с температурой продуктов сгорания и учитывают вторичное стенки в газ:
г ^г ♦ Луг 4-.
*л= £ст. эф. Со [£гС—) -	]'
где Яст Эф = £ст/(1-С1-£стХ1-Лг)]- эффективная степень стенки;
Лг^г- поглощавшая способность газа при температуре
Иногда приближенно принимают _А =(£._+!)/2.
Степень черноты продуктов сгорания
£г=5со +£н о“ссо ен О' 2	2	2	2
Эта формула учитывает уменьшение излучательной способности газов за счет перекрытия некоторых участков спектра излучения.
Геометрия оболочки также влияет на лучистый теплообмен. В оболочках, форма которых отлична от полусферы, длина пути лучей в разных направлениях различна. Излучение газа, заключенного в оболочку произвольной формы, при расчетах лучистого теплообмена можно заменить эквивалентным излучением газовой полусферы. Радиус такой полусферы, равный средней длине пути луча ^Эф, определяется из следующего соотношения:
LA=0.9 ---’
эф г
где И- объем газа, ограниченный оболочкой;
Г- площадь поверхности оболочки.
Практически доля лучистого теплового потока невелика и для его вычисления применяют следующий прием. Вычисляют значение лучистого теплового потока в камере сгорания на достаточном удалении от смесительной головки <9ЛоЛ а затем экстраполируют результат на остальные сечения камеры исходя из правила:
непосредственно у смесительной головки в зоне смешения
*n=0 25-<W
в докритической части при D>1.2
g =g_ [1-12.5(1.2-D)];	С3.19)
в докритической части при DS1.2 9л=<?Ло-
в критическом сечении дл=0.5дЛо;
в закритической части при D<5
Чл=дЛо/С20*У.	СЗ. 20)
D	2	3	4	>5
	0.12	0.05	0.03	0
Учет снижения лучистого теплового потока из-за частичного поглощения в холодном пристеночном слое производят по формуле

(3.21)
где коэффициент у зависит от доли расхода в завесу:
«зав-СО.15+0.2) ча
• 1+0.13) m 'Ъав-0-1'"1
у * 0.7+0.6;
у % 0.8+0.7;
V * 0.9+0.95
(3.22)
Теплофизические характеристики и количество водяного пара и углекислого газа находятся предварительно из расчета рабочего процесса в камере.
3.3. Расчет размывания завесы
Для вычисления конвективного теплового потока необходимо знать значение коэффициента избытка окислителя в каждом расчетном сечении. Рассмотрим случай, когда смесительной головкой образуется по
яс завесы с km, отличным от ядра потока. На начальном участке жидкость испаряется, отбирая тепловой поток на свое испарение, и температура стенки со стороны камеры гарантированно не превышает температуры кипения жидкости. Лучистый тепловой поток также полностью поглощается жидкостью. Дальше тепловой поток определяется составом Си реакционной способностью) газов в пристеночном слое, т.е. соотношением компонентов в нем.
Расчет испарения пленки завесы, особенно турбулентного перемешивания завесы, базируется на эмпирических данных. Применяемые методики дают только весьма приближенные значения в завесе. Считается, что завеса перемешивается только с пристеночным слоем некоторой условной толщины, и справедливо
k	-к	(3.23)
тСТ тЗаВ т тЗЭВ s
где К - соотношение компонентов непосредственно у стенки; шст к начальное соотношение компонентов в завесе;
тзав к - соотношение компонентов после полного смешения завесы с m пристеночным слоем условной толщины;
коэффициент полноты турбулентного перемешивания:
С=1-ехрС-М-х*),	(3.24:
ст где М* (0.05^0.2)—; ^зав х- расстояние от начала завесы, отнесенное к условней толщине пристеночного слоя; /пст~ расход в условном пристеночном слое; расход в завесе; k /	G1
к = -------------------;
m ^^Чав^с? к - соотношение компонентов в ядре потока. Характер изменения по длине стенки имеет экспоненциальный характер. Точное вычисление по приведенной методике связано с ап-
"Ъхл
V =_________
ОХЛ „ г
?охл F
(4.4)
где F - площадь проходного сечения охлаждающего тракта; вычисляется по тривиальным формулам в зависимости от геометрии тракта;
п^хл- массовый расход охладителя.
Эти числа характеризуют соответственно гидродинамическое и теплофизическое подобие течения охладителя. Зависимость Nu=Nu(Re,Pr) находится экспериментальным путем и, аппроксимируется функцией вида
n m Nu=A-Re -Рг ,
(4.5)
где А, п, т- постоянные коэффициенты.
При вычислении всех критериев надо учитывать зависимость теплофизических характеристик охладителя от его температуры.
Для удобства использования выразим из формул (4.2L4.5), приняв, например, 4=0.023, п=0.8, т=0.4 (формула Нуссельта-Краус-
сольда):	^охл^охл^ а=0.023	К, 0.2 <*г	(4.6)
где	с р о.4 о.е К=(-) X . м	(4.7)
Для ряда охладителей комплекс К затабулирован в зависимости от температуры. За характерный линейный размер охлаждающего тракта принимают его гидравлический диаметр, вычисляемый по формуле
4-F
d_=---,	(4.8)
П
где П- периметр охлаждающего тракта. Учитывается весь смоченный периметр, даже если принимается, что часть тракта не участвует в теплообмене.
На практике расчеты ведутся последовательно на участках охлаждающего тракта длиной Дх. На каждом таком участке принимают, чтс происходит охлаждение прямого цилиндра боковой площадью AS, через стенку которого идет тепловой поток величины, средней между величинами теплового потока в начале и конце участка.
4.2. Теплоотдача оребренной стенки
В настоящее время большинство трактов охлаждения ЖРД характеризуется наличием густой сети продольных ребер. Наиболее распространены фрезерованные каналы, гофрированные проставки и трубчатые тракты Срис.4.1).
Рис. 4.1
Примем единую систему обозначений геометрических данных охлаждающего тракта и для линейных размеров укажем размерность величин, используемых в формулах:
D- диаметр стенки со стороны газа (м3;
б- толщина огневой стенки или стенки трубки (мм3;
h- высота тракта охлаждения (м3;
hp- высота ребра [mJ;
(5р- толщина ребра [мм3;
t - шаг ребер по нормали к оси камеры См3;
13 - угол закрутки ребер [рад];
t =t-COS/3- шаг ребер по нормали к стенке ребра или шаг по нормали мехду гофрами [м3;
а- ширина тракта [ м]
Ь- ширина основания гофра, которым он припаивается к стенке Ем];
п - количество ребер по периметру сечения;
бц- толщина внешней стенки [мм];
Dx=D+2c5- диаметр сечения с учетом толщины стенки См];
X - коэффициент теплопроводности материала ребра;
Хст- коэффициент теплопроводности материала огневой стенки.
Присоединенную к ребру поверхность наружной стенки, к которой ^епло может передаваться через ребра, представляют в виде некоторого эквивалентного ребра, являющегося продолжением основного. Его высота равна половине ширины канала между ребрами по нормали к .основному ребру, а толщина- удвоенной толщине наружной стенки (рис. 4.1):
h=a/2;	ба=2би.	(4.9)
Эффективный коэффициент теплоотдачи от стенки к охладителю зависит от конструкции тракта. Увеличение теплоотдачи учитывают введением коэффициента эффективности оребрения рр (здесь и далее индекс “р“ соответствует основному, а "э‘‘- эквивалентному ребру):
5 общем случае 1
г) -------
р cos (3
г 2аЖ где и = -----
h th(^zh)
( г—----------
t (мюр
6 р
(4.11)
t
р
(4.
Ср- коэффициент, учитывающий повышение эффективности основного ребра за счет присоединения наружной стенки. Если теплоотдачей наружной поверхности пренебречь, то (р=1, в остальных случаях (р>1. Эта величина вычисляется по разному для разных геометрических трактов.
Для прямого ребра
pp-tMpW3
Pg-tWMh)p
с = ---------------------------.	(4.13)
р 1+ u/fv	-thCpW-
Р v	р	«3
Для ребра-гофра
Mp-thCph)g
2pg-th(ph)p
С = ------------------------------.	С 4.14)
р 1+ pp/C2pg)- thCMhJp-thCpWg
Здесь, аналогично (4.12)
где kg- теплопроводность материала неогневой стенки.
В случае трубчатого тракта полагают (р=1 и учитывают присоединенную наружную поверхность, как простое удлинение ребра на половину ширины канала: hp=h+a/2, <5э=2бр.
Следует заметить, что начиная с (ph)>3 значение th(ph)^l, и, следовательно, влияние наружной оболочки на эффект оребрения практически исчезает. Это имеет место либо при высоких и нетеплопроводных ребрах, либо при высоких значениях с^. Следовательно, эффективность оребрения за счет теплботдачи в охладитель через наружную стенку повышается, в основном, при коротких, толстых, теп
лопроводных ребрах, омываемых нетеплопроводным охладителем.
Геометрические характеристики трактов различной конструкции выделяются по следующим формулам.
2ля щелевого тракта С без реберЗ:
2<5+/г
F= п D h ( 1 + -у ),	dj,^.
2ля тракта с продольными или винтовыми ребрами:
<5	1
h	пр, ^гл------------—.
N	i+ —
Для тракта с продольными или винтовыми гофрами: tfp	6
1-6/h
2h --------------------------------
. о h~6	/ t +6 -гь г
i+ г-----	/ Г n р 1
t-б / i+ Г----------------1
N р	L 2(h-<5) J
Для трубчатого канала:
1
Г= a-h-n ,	d = 2h -------.
р	l+h/2
(4.15)
(4.16)
(4.17)
(4.18)
4.3. Стандартный порядок расчета проточного охлаждения
В настоящее время применяется методика поверочного расчета проточного охлаждения ЖРД, содержащая следующие шаги.
Шаг 1
Камера разбивается на ряд участков, для каждого из которых отдельно проводится расчет геометрии и параметров рабочего процесса. Составляется сводка геометрических данных о каждом сечении. Задастся начальной температурой стенки камеры С обычно, 1000 К, если ма
териал стенки это позволяет!.
Шаг 2
Пригодится расчет лучистого теплового потока в стенку КС.
Шаг 3
Проводится расчет конвективного теплового потока, например, по методике приведенной в главе 3. Полный тепловой поток в стенку равен <?= 9К*9Л-
Шаг 4
Вычисляют температуру нагрева охладителя в тракте. Подогрев охладителя на t-м участке находится по формуле
АГ = ---------------AS.,	С 4.18)
1 пц -с , г 1 1 pi
где ntl- массовый расход охладителя,
ASt- площадь боковой поверхности прямого цилиндра высотой, равной длине i-гс участка, и диаметром, средним между D и
На каждом участке, начиная от сечения входа до сечения выхода охладителя, суммируют
T»=T»-,+ATl-	U2O:
Теплеэыкссть охладителя ср вычисляется при средней температуре на участке
-х
V--------г-----’ гж.-.*4Г.'г-	!4г-;
Температура жидкого охладителя в сечении выхода не должна превышать допустимую температуру, определяемую из условия отсутствия кипения в жидкости при данных температуре и давлении.
Шаг 5
Определяют коэффициент теплоотдачи, коэффициент эффективности оребрения и скорость жидкости в тракте по формулам (4.4^4.14).
Шаг б
Определяют действительные тепловые потоки и температуру стенки со стороны газа.
Если имеется значение теплового потока, вычисленное при температуре 10С0 К, то из С4.1) можно записать, что
^к ГогЛт
^77 " Гог-1000
С4.22)
для любой Тст. Подставив это соотношение в формулы методики и разрешив относительно qR, получают выражение для уточненной температуры стенки:
ж	^л	ог
б 1	Т-1000
^Ki ооо ^К1 ооо	ОГ
--- + -----
St «хЛ
Т=----------------------------------------------------. (4.23)
тх	1	1000 Гог 1000
"б i 77 т* + тл_-юоо "тТ ЧС1 ооо	Ж ОГ	ОГ
^ОТ Vp
На каждом участке проверяют, что температура стенки со стороны газа не превышает допустимую с точки зрения температурной стойкости материала.
Вычислив действительные значения температуры стенки на каждом участке их подставляют в соотношение С4. 22) для уточнения конвективных тепловых потоков.
Лучистый тепловой поток обычно не пересчитывают, т.к. он мало зависит от температуры стенки.
Если требуется большая точность, уточнение температуры стенки повторяют несколько раз, но обычно ограничиваются первым приближением.
Шаг 7
Определяют температуру стенки со стороны жидкости по формуле, получаемой из (4.1)•
6
^ст ж=^ст г- 7	С4.243
лст
где Тст - температура стенки со стороны жидкости;
Тст г- температура стенки со стороны газа.
Теплопроводность материала стенки X вычисляют при температуре
Т +Т ст ж ст г Т=--------------.	С 4.25)
2
Расчет ведут методом последовательных приближений, пока
I 1^-
При проведении итераций следует проверять сходимость процесса, т.е. условие | T1+S-T1+J |<|	|.
На каждом участке охлаждения проверяют, что температура стенки со стороны жидкости не превышает температуру кипения охладителя при данном давлении больше, чем на 100Л50 градусов. Если охладитель склонен к разложению, допустимый перегрев устанавливают еже меньше.
Шаг 8
Рассчитывают гидравлическое сопротивление охлаждающего тракта Сем. ниже).
Эта методика ориентирована на минимальное количество ручных вычислений. Так, вместо формулы (4.24) рекомендуется пользоваться •некоторыми графиками, приведенными, например, в [10, 13). Понятно, что такой подход не способствует получению результатов высокой точности. Кроме того, в формуле (4.19) приходится часто делить друг на друга числа одного порядка, что существенно уменьшает точность, даже при использовании ЭВМ. Исходя из этого, предлагается другая методика расчета охлаждения, направленная на использование ЭВМ на всех стадиях расчета и обеспечивающая меньшую погрешность.
4.4. Алгоритм расчета наружного охлаждения на ЭВМ (18)
Далее с краткими пояснениями, без вывода и ссылок на предыдущие разделы приведены формулы, используемые для составления программы
расчета охлаждения ЖРД. Общая идея алгоритма остается прежней- задаться распределением температуры на стенке камеры и последовательно вычислять тепловой поток в стенку, подогрев жидкости и уточнять температуру стенки, пока ее изменение не станет меньше заданной точности.
Для рассмотрения алгоритма зададимся следующими посылками, данными в описании шага 1 в разд. 4.3.
Камера разбита сечениями с достаточно мелким шагом, выбираемым исходя из требований к точности расчета; для каждого сечения, располагающегося на расстоянии xt от смесительной головки, задан диаметр Di.
Газ течет от 1 сечения к последнему.
Охладитель течет от i-го сечения к j-му, причем для определенности положим j<i.
Алгоритм расчета может содержать следующие шаги.
Шаг 1
Вычисление газодинамических параметров потока при заданных рк, зависимости Л4Сх4). Получаем в каждом сечении W. , \ , ядра зав	х 1	11
TKi , р*, р* и прочие необходимые параметры. Удобнее всего вести расчет с помощью газодинамических функций. Так, Х.£ находится из решения любым численным методом трансцендентного уравнения ?Ш = CD/d )z. Заметим, что в качестве интервала поиска решения на участке 1хорсКр[ следует задавать COu.lt, в критическом сечении Х=1, в закритической части X. и ] lu.X. С,
Температура потока вычисляется, исходя из соотношения компонентов:
аядра ^ядра* ^ав Лзав ----—----------;----------,	(4.27)
'^ядра* ^ав
где \ядра зав)" Расход компонентов в ядре потока и на завесу,
“(ядра зав)" коэффициент избытка окислителя в ядре и завесе.
С достаточной точностью можно задаться а
лДра
Шаг г
Задание начальной температуры стенки. Часто рекомендуется задавать температуру стенки всей камеры одним значением, например, 1000 К. Для ускорения сходимости итерационного процесса уточнения тепловых потоков лучше задать ее в интервале [500^1000] К с использованием, например, интеграла вероятности:
1	6х1	6х1
ТСТо=( 1600^2100)------ехр(------(640*840) ехр(- --------
/	2	2
«/2л
I Х1 хкр где 6х=---------------------
1 тах(х -х„_,х -х,._) а Кр о Кр
(4.28)
(4.29)
Последняя формула, в среднем, уменьшает количество итераций на 2*4, что заметно при проведении расчетов на "медленных" ЭВМ, особенно не оснащенных сопроцессором.
Шаг 3
Вычисление удельного теплового потока.
Конвективный удельный тепловой поток в стенку камеры рассчитывается по второй приближенной формуле Иевлева;
S
-------- [кВт/м2]. О . 38 Рг
При инженерных расчетах часть множителей (S,B,Pr) принимается зависящей от коэффициента избытка окислителя в пристеночном слое и от температуры стенки или от относительной температуры стенки Тст=
а множители D, cL_, p-Xl-fl2) = т(А.) определяются геомет-ч* 1 U1	л
рией камеры сгорания и расходом продуктов сгорания. Практически удобно разделить последнюю формулу на два независимых сомножителя:
2	0.88
1~р Рк
S
В = В-----------.	(4.31)
q	о. ве
Рг
Комплекс зависит только от геометрии камеры и используемого топлива и может быть рассчитан только один раз. В^ является функцией используемого топлива и температуры стенки и должен пересчитываться с каждым уточнением температуры.
В формулу для В^ при вычислении критерия S необходимо подставлять Озазесы' вычисляемое по методике разд. 3.3.
Лучистый тепловой поток также рассчитывается во всех сечениях по известным методикам, либо задается в одном сечении камеры сгорания и аппроксимируется по правилам, приведенным в раз. 3.2.
Шаг 4
Вычисление подогрева охлаждающей жидкости. Подогрев жидкости между сечениями i и у рассчитывается по формуле:
1 дтх= --------- f q dS,	(4. 32)
^охл’сохл
где %хл” массовый расход охладителя в сечении; сохл” теплоемкость охладителя (ср в случае газа).
Интеграл
Я= J q dS	(4.33)
можно приблизить следующим образом:
q. +Q D +D
Q= -——- л |х.-х |.	(4.34)
2	2	1 J
Здесь поверхность приближается прямым цилиндром диаметра, среднего между диаметрами сечений, и тепловой поток рассчитывается по средней величине:
D1+DJ r----------s-------~
Q=------------n /(x-xJ +(D +D.) /4 .	(4.35)
2	2	1 J 1 J
Здесь точно вычисляется боковая поверхность усеченного конуса, образованного линейным приближением профиля камеры между i и j сечениями, а тепловой поток берется средним.
Интеграл (4.33) может также быть проинтегрирован при аппроксимации профилей камеры и удельных тепловых потоков полиномами более высоких степеней. Однако при достаточно частом разбиении камеры сечениями, особенно в области критического сечения, где велика кривизна стенки камеры и большая производная удельного теплового потока по осевой координате, предлагается с удовлетворительной точностью использовать формулы (4.34) и (4.35).
Теплоемкость охладителя, как правило, принимается зависящей от его температуры. Поэтому необходимо производить расчеты по формуле (4.32) итеративно, каждый раз заново вычисляя теплоемкость при новой температуре охладителя.
Для каждого охладителя существует критическая температура, нагрев выше которой приводит либо к вскипанию, либо к химическому разложению охладителя. На этом шаге расчета необходимо контролировать перегрев охладителя выше критической температуры.
Шаг 5
Вычисление коэффициента теплоотдачи от стенки к охладителю производится по формуле:
^охл^охл^ а= 0.023-----------------К-т) .
0.2	'р
Методика расчета величин, входящих в эту формулу, описана в разд. 4.1.
Для ряда охладителей комплекс К затабулирован в зависимости от температуры и может быть вычислен методами интерполяции. Это существенно уменьшает время счета и, в ряде случаев, может привести к увеличению точности, особенно если значения ср, Хи р также вычисляются интерполяцией по графикам или приближаются аппроксимационными зависимостями.
Шаг 6 *
Вычисление температуры стенки со стороны охладителя.
Известно, что 90ХЛ=а0ХЛ*-Тж) и 90ХЛ=<?. откуда легко находится гст ж.
Шаг 7
Вычисление температуры стенки со стороны газа. Тепловой поток через стенку камеры определяется по формуле
Чт
^7~(Гст.г-Гст.г	“-36’
° ст
Эта зависимость справедлива для плоских стенок, но обычно используется и для тонких цилиндрических, для которых справедливо бгмг« R, где R- радиус стенки. Более точной является формула
С1
Х„ СТ„ -Т„ _) ст СТ. Г СТ. X
где 2г- внутренний диаметр стенки, равный ;
2₽- внешний диаметр стенки, равный (.D^6QT).
Применение последней формулы дает в диапазоне • разумных толщин стенки и диаметров камеры очень небольшой выигрыш в точности. Ею можно рекомендовать пользоваться при расчете микроЖРД со сравнительно толстыми стенками и маленькими диаметрами камеры и сопла.
Из этих формул находится температура стенки со стороны газа. Поскольку зависит от температуры материала, процесс определения температуры итерационный и ведется последовательным вычислением J/2) и _, пока не будет достигнута заданная точность. При вычислениях следует контролировать сходимость процесса.
Шаг 8
Вычисляется удельный тепловой поток, анологично описанному на шаге 3 Если абсолютная погрешность вычисления теплового потока на предыдущем шаге
l?k.t-Qkl б9= ----------,
где тепловой поток на k-1-й итерации;
qk - тепловой поток на k-й итерации, в каком-либо сечении больше допустимой величины, шаги алгоритма 4^8 повторяются.
Шаг 9
Повторить шаги с 4 по 8 для всех сечений участка охлаждения.
Шаг 10
Для расчета прочих параметров системы охлаждения, например, гидравлического сопротивления, можно вычислить критерий Re и другие необходимые величины во всех сечениях.
Шаг 11
Вывод результатов на экран или в файл.
Этот алгоритм является ЭВМ-ориентированным, что отражается на его структуре и содержании. Прежние методы, ориентированные на ручные расчеты с применением графиков, таблиц и номограмм, организовывались так, чтобы производить минимум подсчетов вручную в ущерб точности. Использование ЭВМ позволяет провести множество итераций за сравнительно короткое время (исчисляемое секундами), что существенно повышает точность расчета.
4.5.	Инженерная методика расчета гидравлических потерь в тракте охлаждения
Гидравлические потери складываются из потерь на преодоление сил трения по длине канала и местных сопротивлений, возникающих при обтекании выступов и изменении направления движения потока, а также при резких изменениях режима течения в тракте.
Сопротивление трения обусловлено вязкостью охладителя и составляет 8CL90 % от общих потерь. Потери давления на участке безотрывного течения в канале охлаждения рассчитываются по формуле
г
I ^ОХЛ^ОХЛ
г
(4.38)
где I- полная длина канала; для закрученных каналов необходимо учитывать угол закрутки /3;
коэффициент сопротивления трения. Зависит от режима течения охладителя и определяется по формулам: 64/Re- для ламинарного потока;
Д 68	0.29
С = 0. И Г — ♦ — I - для турбулентного потока. L Re J
Существуют и другие зависимости ((A'd^fle), где Д- высота шероховатостей на стенке, зависящая от чистоты обработки поверхности.
5. ЧИСЛЕННЫЕ МЕТОДЫ ДЛЯ ПРОВЕДЕНИЯ РАСЧЕТА ВЕЛИЧИНЫ КОНВЕКТИВНОГО ТЕПЛОВОГО ПОТОКА
Существуют различные подходы к облегчению многочисленных расчетов, потребных для вычислений по приведенным выше методикам. Раньше широко использовались многочисленные графики и номограммы, позволяющие практически мгновенно найти результат, правда, с весьма низкой точностью, зависящей от многих факторов: типографского качества издания, аккуратности и остроты зрения вычислителя и т.п. При использовании ЭВМ разумно отказаться от такого подхода и возложить задачи объемных вычислений на нее, применив соответствующие численные методы расчета по тем или иным алгоритмам.
5.1.	Численное интегрирование методом Симпсона
Для вычисления делителя Jp3dx в первой приближенной формуле Иевлева необходимо взять определенный интеграл от функции, включающей относительный диаметр D(x)=D(x)/d , заданный таблично, итерационно уточняемую температуру ТС7=Тст/Гог, которую невозможно точно выразить аналитически. Таким образом, найти первообразную рз не представляется возможным и нужно использовать методы численного интегрирования. Как правило, интегрируемую функцию для этого приближают некоторой кривой на множестве нескольких точек, для которых интеграл легко находится. Т.Симпсон (1710-1761) предложил использовать кривую второго порядка и вывел для нее общую формулу интеграла для кривой с равноотстоящими на расстоянии h узлами:
1 = 1 1=1
(5.1)
Здесь /	- значение функции в точке Сх +х4 )/2.
1-1 /2	11-1
В случае переменного шага, который более применим для данного приложения, формула (5.1) принимает вид
•Точность формулы Симпсона весьма высока и, при условии непрерывной производной f до четвертого порядка соблюдается неравенство
т V b“a	4
Д7< max(/IV)------max(h.).
a,b 2880 a,b 1
Аналогичные (5.1 ) формулы получены для приближения подынтегральной функции прямоугольниками или трапециями. Их точность существенно ниже:
п	/ » b~d г
-пахСЮ, M<maxCf )--------------maxChJ. (5.2)
i=i х’1/2 a,b 1	a,b 24 a,b 1
n Л -1 +A	t . b”G 2
/= £--------max(h ),	Д/< max(/<3>)-------maxChJ. C5.3)
i =1	2 a,b 1	&,b	12 a,b 1
Следует отметить, что формулы (5.2) и (5.3) имеют второй порядок точности и являются точными для кривых первого порядка, а формула Симпсона- четвертый порядок точности и является точной для кубических многочленов. Для расчета рз на стенке камеры сгорания вполне применима формула прямоугольников, тогда как на стенке сопла формула Симпсона дает значительно более высокую точность.
Примеры программ для вычисления определенного интеграла методами прямоугольников и Симпсона приведены в приложении 3.
5.2.	Метод наименьших квадратов
Для вычисления комплекса 5(Тст) необходимо описать аналитически вид известного графика, полученного экспериментальным путем. Часто в таких случаях искомую функцию приближают полиномом некоторой степени, а для поиска его коэффициентов используют метод наименьших квадратов (МИК). Суть метода состоит в том, чтобы, приближая функцию /(х) полиномом Рю(х), где
m
Р (,х) = а +а х+а х+... хт= Г а. х1 ,	(5. 4)
т	0	12	т *-• i
о
найти такой вектор коэффициентов а=<ао,ах ,а ,... ,ат>, при котором среднеквадратичное отклонение
/1 n
<5CF ,/) = / - E CP (x)-/(x))2
m	n-1 i=o m
было бы минимальным. Поиск таких коэффициентов осуществляется решением системы линейных арифметических уравнений (СЛАУ)
d
--- 6СР ,/)=0 для i=0,. . . ,ш.	г da. m
i
Матрица системы таких уравнений называется матрицей Грама и для полиномиального приближения имеет вид:
	п+1 п	п Exi 1 =0 п	п Ех? 1 =0 п Е<	п Ё xj*+1	п Ь,/.	
г =	1 =0	1 =0	1 =0	1 =0	i =о	. (5
	п Е 1 =0	п Е<‘ 1 =0	Ёх^+г ... 1 =0	Ёх?т 1 =0	1 =0	
Для формирования этой расширенной матрицы достаточно вычислить только первую строку и два последних столбца, а остальные элементы не являются "оригинальными” и заполняются с помощью циклического присвоения.
Хотя МНК кажется весьма привлекательным своей простотой реализации, он таит в себе некоторые серьезные недостатки. Например, уже при т>5 система С5.7) очень часто оказывается плохо обусловленной, что приводит к большим вычислительным погрешностям. Поэтому без специального выбора базиса полинома не рекомендуется увлекаться повышением его степени- это может привести к составлению "почти линейно зависимой” системы, и параметры модели будут сильно искажены ошибками вычислений. Метод также чувствителен к погрешностям задания значений функции в узлах, которые неизбежны при использовании экспериментально полученных данных.
Программа МНК, приведенная в приложении, содержит подпрограмму решения СЛАУ методом Гаусса, которая может быть использована отдельно.
Особо стоит вопрос рационального вычисления значения полинома. Непосредственный алгоритм типа Рт(х) = а0+ах х+агх+... -t-o^x"1 является очень неэкономичным, так как требует nt возведений в степень, и на практике часто заменятся на разложение Горнера:
Р (х) = а +х (а +х (а +.. .х (а +а х)..).	(5. 8)
m	Ota	m-i т
5.3.	Сплайн-интерполяция
Часто встает задача вычислить значение функции, заданной таблично на некотором множестве, в точке, не совпадающей ни с одной течкой области определения. Такая проблема решается методами интерполяции. то есть приближенной замены сложной функции (например, Рг(а )) на более простую, которую нетрудно вычислить при любом значении аргумента. Известно много методов интерполяции- Лагранжа, Ньютона, Эйткена, каноническим полиномом и т.п. Для инженерных методик давно применяется метод сплайнов, ведущий свое начало от гибких лекал, которыми чертежники приближают сложные графические контуры. Гибкая линейка устанавливается таким образом, что проходит через все узловые точки, и по ней проводится гладкая кривая. Уравнение деформированной балки, как известно из курса сопротивления материалов, удовлетворяет равенству
pIV(x)=O.	(5.9)
Функцию р(х) и будем использовать для аппроксимации зависимости /Сх). заданной в узлах интерполяции хо,... ,хп значениями /о,..,/п. £оли в качестве р(х) выбрать полином, то его степень должна быть не выше третьей, и он будет записан в виде
(х)=а Сх-х )+с. (х-х. )2+d. (х-х4 )3,	(5.10)
j 1	i-i 1	i-i	1	i-i
где a,fc,c,d- коэффициенты сплайна, определяемые из условий: 13 равенство в узлах интерполяции (условие Лагранжа):
23 непрерывность первой и второй производных в узлах,
3) условия “свободных концов линейки1*:
<р Сх Сх )=0, Г1 о rn п
или другие граничные условия, определяющие поведение функции за пределами отрезка задания узлов [153.
При составлении СЛАУ по указанным выше условиям получают трех-диоганальную систему уравнений относительно ci:
Г с =0; 1
=	1=2.....п;
Сп+1=0’
которую решают методом прогонки [16], и формулы для остальных коэффициентов имеют вид;
ь1 =( A "A _t	-Сс1 +1 +2c1	/3,	(5.12)
^Сс^ЖЗ^).
Затем по известным коэффициентам можно вычислить значение сплайна в любой точке х, ближайшей к некоторому узлу х. по формуле (5.10). Если точка х лежит за пределами задания функции, то имеет место экстраполяция функции и (в случае использования условия свободного конца линейки в соответствии с условием 3) проводится линейная экстраполяция.
Этот метод также не свободен от недостатков, как и все прочие методы интерполяции, включающие условия Лагранжа. При некоторой погрешности задания значений функции в узловых точках “гибкая линейка” начинает приближать совершенно невероятные кривые колебательного характера с большими забросами от заданных значений. Поэтому при приближении неточно заданных функций с гарантированным характером получаемой аппроксимации рекомендуется использовать МНК, а не интерполяционные методы.
5.4.	Многомерная интерполяция
Комплекс S зависит от двух переменных для каждого рода топлива:
от аот и Тог. Таким образом, встает задача двумерной интерполяции экспериментально полученной зависимости.
Задача многомерной интерполяции гораздо сложнее, чем для функции одной переменной. Для выполнения именно интерполяции, а не экстраполяции, необходимо, чтобы точка, в которой вычисляется приближенное значение, находилась внутри области задания исходной функции. Само определение того, является ли точка внутренней для некоторой, в общем случае, невыпуклой области, нетривиально и требует определенных усилий. Однако, для этого приложения известно, что задана на прямоугольной (выпуклой) области аргументов, и задача упрощается.
Для многомерных случаев разработаны интерполяционные формулы Лагранжа и Ньютона нескольких аргументов, но они весьма объемны, требуют больших вычислительных затрат и очень чувствительны к погрешностям задания значений функции в узловых точках. Поэтому на практике чаще всего применяют последовательную одномерную интерполяцию по каждой переменной. В ручных расчетах использовалась линейная интерполяция, а применение ЭВМ дает возможность по одной или нескольким переменным использовать более точные методы.
Следует отметить, что если для переменной использовалась сплайн-интерполяция для вычислений значений SCa) при различных температурах, то не следует использовать тот же метод для вычисления SCT), так как накопившаяся погрешность, скорее всего, приведет к списанному выше эффекту.
5.5.	Решение нелинейных уравнений
При вычислении тСХ) необходимо предварительно найти коэффициент скорости X. Он находится из газодинамической функции
F	di 2
q(X) = — = [ — 1 = Q.	(5.13)
Значение Q в любой точке может быть найдено исходя из геометрии КС. Встает вопрос о решении нелинейного (трансцендентного) уравнения с(Х)=С. Для поиска эго решения нужно предварительно указать область определения этой функции. На участке камеры сгорания и в докритической части сопла л изменяется от 0 до 1, а в закритичес-кой части сопла X больше 1 и меньше X max
где
(5.143
Будем искать решение уравнения
/(Х) = д(Х) - Q = О
(5. 15)
методом секущих. Суть метода в том, что в двух близких точках находятся значения и /(\+1) и по ним строится прямая. Следующее приближение Хк+2 является точкой пересечения полученной прямой с осью абсцисс. Таким образом получаем общую формулу метода секу-
щих:
< =< “ ----------------/(Хь).
k+1 * /(xk)-/(xk_a) k
(5.16)
Процесс прекращается при |Xk+i -Xk|<с, где с- заранее заданная точность вычислений.
Метод секущих удобен тем, что начальные приближения Хо и X* могут находиться с разных сторон от корня, что позволяет в докритической части камеры задать Хо=0+с, \=1-е, а в закритической части задать Хо=1+г,	В критическом же сечении Х=1. Этот метод
существенно превосходит метод дихотомии (деления отрезка поиска корня пополам) в скорости сходимости, хотя требует несколько больших объемов вычислений' на каждом шаге для поиска следующего приближения.
Литература
1.	Алемасов В. Е., Дрегалин А. Ф., Тишин А.П. Теория ракетных двигателей/ Под ред. В.П. Глушко.- М.: Машиностроение, 1989. 464 с.
2.	Основы теории и расчета жидкостных ракетных двигателей. В 2 кн. Кн. 1. Учеб, для авиац. спец, вузов/ А.П. Васильев, В.М. Кудрявцев, В.А. Кузнецов и др.; Под ред. В.М. Кудрявцева. - М.: Высшая школа, 1993. 383 с.
3.	Шевелюк М. И. Теоретические основы проектирования жидкостных ракетных двигателей. - М.: Оборонгиз, 1960. 684 с.
4.	Термодинамические и теплофизические свойства продуктов сгорания: Справочник в 10 т. / Под ред. В.П. Глушко. -М.: ВИНИТИ АН ССР, 1971 - 1979.
5.	Глушаков А. Н., Рачук В. С., Рубинский В. Р., Фалеев С. В. Термодинамический расчет камеры жидкостного ракетного двигателя: Методические указания по выполнение курсовой работы по дисциплине "Теория двигателей" для студентов специальности 130400 дневного обучения. - Воронеж : ВГТУ, 1995. 33с.
6.	Глушаков А.Н., Кретинин А.В. Программа вычисления корней исходной системы уравнений для расчета равновесного состава ПС. -Воронеж : ВПИ, 1993. 16 с.
7.	Мелькумов Т. М., Мелик -Пашаев Н. И., Чистяков П. Г. Ракетные двигатели. - М.: Машиностроение, 1968. 511 с.
8.	Дубинкин Ю. М. Кондрусев В. С., Фрейдин А. С. Тепловой расчет и проектирование камер ЖРД. - Куйбышев : КуАи, 1985. 75 с.
9.	Добровольский М.В. Жидкостные ракетные двигатели.- М.: Машиностроение, 1968 . 396 с.
10.	Булыгин Ю. А. , Кретинин А.В., Скачилов В. Н. Расчет наружного охлаждения ЖРД; Методические указания к выполнению курсовой работы по теплопередаче для специальности "Ракетные двигатели". - Воронеж: ВПИ, 1992. 35 с.
И. Орлин С.А., Пелевин Ф.В. , Поляев В.М. Расчет теплового состояния камеры двигателя.- М.: МВТУ, 1987. 16 с.
12.	Березанская Е.А., Курпатенков В.Д., Надеждина Ю. Д. Расчет конвективных тепловых потоков в сопле Лаваля.- М.: МАИ, 1976. 76 с.
13.	Березанская Е.А., Курпатенков В.Д., Надеждина Ю. Д. Расчет наружного проточного охлаждения.- М.: МАИ, 1977. 52 с.
14.	Поляев В. М., Майоров В. А., Васильев Л. Л. Гидродинамика и теплообмен в пористых элементах конструкций летательных аппаратов.-М.: Машиностроение, 1988. 168 с.
15.	Амосова. А., Дубинский Ю. А., Копченова Н. В. Вычислительные методы для инженеров. - М.: Высшая школа, 1994. 544 с.
16.	Мудров А.Е. Численные методы для ПЭВМ на языках Бейсик, Фортран и Паскаль.- Томск: МП “Раско", 1991. 272 с.
17.	Пелевин Ф. В., Леонтьев С.Н. Интенсификация теплообмена в кс-льцевом канале.- М.: МГТУ, 1994. 16 с.
18.	Затонский А.В., Орлин С.А., Пелевин Ф.В. Расчет теплового состояния камеры ЖРД с использованием ЭВМ. - М.: МГТУ, 1996. 93 с.
70 ОГЛАВЛЕНИЕ
Введение............................................ 3
1.	Термодинамический расчет камеры КРД..................4
1.1.	Состав компонентов топлива...........................5
1.2.	Соотношение компонентов топлива......................8
1.3.	Состав продуктов сгорания...........................10
1.4.	Исходная система уравнений для расчета равновесного состава.....................................12
1.5.	Полная энтальпия топлива и продуктов сгорания.......17
1.6.	Определение параметров продуктов сгорания на входе в сопло............................................18
1.7.	Расчет параметров продуктов сгорания на выходе из сопла.................................................20
2.	Расчет основных параметров камеры двигателя.........22
2.1.	Определение теоретических параметров камеры.........22
2.2.	Определение основных параметров камеры..............26
2.3.	Построение газодинамического профиля камеры.........31
3.	Расчет наружного охлаждения.........................36
3.1.	Расчет конвективных тепловых потоков................36
3.2.	Расчет лучистого теплового потока...................40
3.3.	Расчет размывания завесы............................42
3.4.	Расчет радиационного охлаждения.....................44
4.	Поверочный расчет наружного охлаждения..............45
4.1.	Основные расчетные соотношения......................45
4.2.	Теплоотдача оребренной стенки.......................47
4.3.	Стандартный порядок расчета проточного охлаждения.... 50
4.4.	Алгоритм расчета наружного охлаждения на ЭВМ........53
4.5.	Инженерная методика расчета гидравлических потерь в тракте охлаждения......................................5S
5.	Численные методы для проведения расчета величины
конвективного теплового потока......................61
5.1.	Численное интегрирование методом Симпсона...........61
5.2.	Метод наименьших квадратов..........................62
5.3.	Сплайн-интерполяция.................................64
5.4.	Многомерная	интерполяция.............................85
5.5.	Решение -нелинейных	уравнений........................36
Литература............................................68
Оглавление............................................7С
Приложение............... ............................72
72
ПРИЛОЖЕНИЕ 1
Число электронов во внешней оболочке и валентность химических элементов
Элемент	Н	Li	Mg	Be	в	А1	С	Si	0	F	Cl
п	1	1	г	г	3	3	4	4	6	7	7
Д	1	1	г	2	3	3	4	4	2	1	1
Продукты полного сгорания
Г. э.	н i	Li	Mg	Be	В	Al	c
О. э.							
0	Н 0 2	Lis0	MgO	BeO	в 0 a a	Al 0 a a	co a
F	HF	LiF	MgFa	BeFa	BF 3	A1F 3	CF 4
С1	НС1	LiCl	MgCla	BeCla	BC13	A1C1 a	CC1 4
Средний показатель изоэнтропы для ракетных топлив при Р «15,0 Ила; 6 » 300 и при оптимальном соотношении компонентов
Окислитель	Горючее	а	К
J2 S	Этиловый спирт 95 %	0,9	1,140
	Керосин	0,8	1,146
	НДМГ	0.9	1,136
	NH	1,0	1,188
	з а		
	н	0,7	1.194
	г а		
	Н (6=3000)	0,7	1,214
	га		
AK + 27 % AT	Керосин	1.0	1.148
AK * 20 % AT	НДМГ	0.95	1,176
AT	НДМГ	0.95	1.156
AT	Аэрозин-50	0.95	1.174
98 % H 0 a a	Керосин	1.0	1.151
Теоретическое соотношение компонентов
Окислитель	Горючее		j х	1
Кислород	Водород Керосин	। НДМГ Гидразин Пентаборан i	п 94	j 3 39	* 1	1. 77 13 7 4
Азотный тетраксид	Керосин НДМГ Гидразин Тонка	4' 26 j я
Азотная кислота	Керосин Тонка НДМГ Пентаборан Гидразин	А 3^	1 1б7 3.34	1 4.80	1 1.97	j
96% Азотная кислота	Керосин	5.57	|
Тетранитрометан	Керосин	6. S2 j
Перекись вЛюрода	Керосин Гидразин Пентаборан	7.20	I 6^45	
90% Перекись водорода	Керосин	-г
Фтор	Водород Гидразин Пентаборан Аммиак Керосин	► * (ХХа) 1(\)00 ( >UH\) (a)(1) ( И СП .г* ООО
Дифторид кислорода	Керосин Гидразин НДМГ	с ^с: . — 3. 60
Трехфтористый хлор	Гидразин Пентаборан	СХ) (0 гпг (\io6
Состав продуктов сгорания для групп топлива ЖРД
; Химические	значения температуры 11U в КС Тк, к				
	1300-dUUU	]		1	300U-3QOO I		1	>3UU0
& Лт5вЛХ5П х ш :• группы топлива	значения коэффициента избытка окислителя а				
	>1	<1		<1	<1
	На О	НаО	HaO	HaO	HaO
'ч	СОа	СОа	СОа	COa	COa
3	Оа	СО • На	CO Оа На* of H*	CO Oa. Ha* OH* O’ H*	CO °a. Ha* of H*
1	'-4	На U	. НаО	HaO	HaO	HaO
	СОа	СОа	COa	COa	COa
J	Оа	СО	CO	CO	CO
N	Na	Na. На* С . сн**	N». Ha* NO* oa* OH* ° z н! N *	Na, NO* Oa* OH* 0 Z H Z N *	Na. bZ • NO* OH^ s N *
'	н	01а НС1 Оа	01а HC1 Оа	UOa;HaO Na;OH* Oa;NO*	UOa; HaO CO;OH* OH*.	UOa;HaO CO;OH* OH*.
	N	НаО	НаО	ClajCl	Ha;Oa*	Ha;Oa*
1	*41	Na СОа	Na СОа СО	HS1	NO*. Cla* F	« «
	Ка	Fa	UOa; CO	UOa; 00	UUz;UU
	HF	HF	HaO; Ha	HaO;Ha	HaO;Ha
	CF*	CF*	Na; S'	Na; HF	Na;HF.
5	СОа Оа НаО . Na	СОа Ha НаО Na	* ж ou.o* ci • *	T	O*;QH* №
Примечание: в первом приближении парциальные давления газов, отмеченных могут быть приняты равными нулю.
Характеристики двухкомпонентных топлив С б =70?
Окислитель	Горючее		Тк, к	Та, К	рем	рг, кг^м3 j
Оа	На	4,U0	ЗУ ( (	13ЬЬ	10,08	...	ж...	1
	Керосин	3,70	ЗЬ8Ь	8457	ЗЬ. 7У	1ОЬ ”	]
	Nan*	0,У8	340b	1У47	30, Ь4	4 ’ ’>*4^**	’
	НДМГ	1,70	ЗЬ08	3380	83,14	
	Аэрозин ЬО	1,37	3500	—		1 1	_	1
	Хайдин	 Гуо	зьоо	—	— : -	
	ЬеНа	1 ,33	ззьз	38Ь7	1У,У1	85	
	А1Нз	 0,7У	’4301"	3333	ЗУ,4Ь	1 50- _
	НвН©	8,18	4150	ЗУЬУ	ЗЬ.ЗЬ	ЗУ
	ННз	1,41	зоьь	"1710	1У.50	8У0
	ОаНвОН	" 1,78	3411	330b	34, ЬО	УУО "
Na 0*	На	ь,зь	^Ь40	1105	11,30	зьз
	Na И*	1,33	3347	"1703	31, ьз	1817
	НДМГ	3, Ь'7	341b	1УЬЬ	34,73	1170
	Аэрозин ЬО	з,зь	ЗЗЬЗ	30У8	—	-
	ммг	3,Ь0	ззьО	—	—	—
	Хайдин	3,30	3370	—	—	—
	Керосин	4,00	3438	301Ь	ЗЬ,8Ь	13УЬ
	ЬеНа	1,08	3831	3403	1У,17	УЗО
	• А1Нз	1,1а	417У	ЗОУО	30,40	14Ь7
	Ив На	3,00	ЗУ 13	зььЗ	37,31	1084
ШЮз	На	5,14	3474	1V43	13,о4	348	!
	NlHt	 Гьо-	3031	1530	31 ,ЗУ	13Ь4
	НДМГ	з, оО	3147	174Ь	34,38	133 г
	Гонка-ЗЬО	4,Ь0	Зооо	—	ЗЬ,50	
	Керосин	"4,88"	3147‘‘	1838	ЗЬ,73	
	ЬеНа	3,33	ЗЗЬ8	8815	81,88	
	АГНз	1,04	ЗУУЗ	8У14	38, bl	14У5
	ЬвНв	"8,00	ЗЬ88	"88Г1 ‘	35,15"'	‘w
на0а	На	7,аа	341У	' 1050”	11,71	43Р
	Керосин	Ь,ЬУ	300b	174Ь	33, Ь7	1342
	Na На	3,03	ЗУ37	1ЬЗЗ	1У,77'"	1 1 ;
	НДМГ	3,00	3147	174b	34,38	" 1^/5 _•<	j
	ЬеНа	1,44	ззо5	ЗЬЬ8	17;оу 	У8С
	А1Нз	0,У8	3834	3707	34,33	14bb
	Из Но	3,33	ЗУЬУ	80У4'	1У,5У	103'	i
76
ПРИЯОЛЕНкЕ 2
ТЕПЛОФИЗИЧЕСКИЕ СВОЙСТВА ГАЗОВ .
Газ	Т.1К.1	gio* ЛН/М8]	ХЮ’ДВт/См-ЮГ	Ср,tКдж/Скг-К)3
	40U	19. ЗУ	23.00	О.УзУ
	50С	27.06	40.53	1.076
	800	33.41	54.22	1.170
	1000	39.15	67.02	1.237
	4.200	44.29	78.84	1.284
	1400	49.07	89.73	1.317
i	1500	53.61	99.86	1.343
	1800	57.94	109.41	1.364
	2С0С	62.09	118.49	1.380
	2200	66.09	127.20	1.393
	2400	69.96	135.62	1.405
	2600	73.71	143.74	1.415
	2800	77.35	151.70	1.424
	3000	80.90	159.44	1.431
	400	22.17	33.33	1.049
	600	29.61	46.77	1.088
	800	35.84	58.99	1.140
	1000	41.53	70.72	1.186
	1200 '	46.82	81.86	1.221
	1400	51.78	92.45	1.248
	1600	56.48	102.50	1.268
	1800	60.97	112.04	1.284
	2000	65.28	121.17	1.296
	2200	69.43	129.84	1.305
	2400	73. 45	138.17	1.314
	2500	77.35	146.21	1.320
	2800	81.14	154.00	1.326
	3000	84.83	161.53	1.331
		17.62	47.98	1.904
	i O'JiJ	26.45	76.20	2.018
	><' r<	34.88	106.39	2.153
1	i •	42.50	137.08	2.294
	\ 120C	49.25	167.52	2.422
	1420	55.22	197.29	2.538
	i -500	60.59	Й26 22	2.638
	: isoc	65.55	254. 40	2.723
		70.29	284.38	2.793
j	i 2200	74.95	308.88	2.851
	1 2400		335.71	2.899
	250C	84 58	362.59	2.940
	' 2800	_ -	390. 52	2.977
	! зон		414.07	3.004
Газ	Т, [KJ	р-10в ДН/м2]	ХЮ’ДВт/См-К)]	Ср,[Кдх/Скг-Юj
Ha	400 600 800 1000 1200 1400 1600 1800 2000 2200 2400 2600 2800 3000	10.7У 14.10 17.02 19.68 22.14 24.46 26.66 28.76 30.78 32.73 34.61 36.44 38.22 39.96	&(. 14 297.11 363.14 428.37 493.73 559.38 625.04 690.18 754.46 817.60 819.40 939.81 999.02 1057.05	14.600 14.671 14.818 15.113 15.506 15.942 16.377 16.788 17.156 17.489 17.790 18.065 18.314 18.546
Na	400 600 800 1000 1200 1400 1600 1800 2000 2200 2400 2600 2800 3000	21.62 28.31 34.20 39.56 44.54 ' 49.21 53.65 57.89 61.95 65.88 69.67 73.36 76.94 80.44	32.41 44.38 55. 73 66.57 77.04 87.01 96.59 105.72 114.39 122.68 130.59 138.21 145.62 152.78	1.045 1.076 1.123 1.169 1.206 1.234 1.256 1.273 1.286 1.297 1.306 1 Q1 -5 1'. за: 1 325
Теплофизические свойства жидкостей
Жидкость	Темпер, т, скз	Плотность [кг/м®]	Динамич. вязкость и-0.001 [Па-с]	Коэф, теплопров.	Удельная теплоемкость
Керосин	273 323 373 423 473	825 788 750 712 676	3.00 1.08 0.47 0.20 0.12	0.156 0.144 0.133 0.121 0.109	1.84 2.03 2.22 2.41 2.60
Спирт	273 323 373 423 473	1"У40‘ 750 660 570 480	2.12 .0.76 0.34 0.16 0.10	0.186 0.193 0.200 0.207 0.214	2. 64 3.13 3.98 4.69 5.36
Азотная кислота	273 323 373 423	1000 1420 1340 1260	1.6У 0.67 0.40 0.27	б. 334 0.324 0. 312 0.295	1.84 1.88 1.95 2.01
Вода	273 323 373 423 473	1000 990 960 915 870	1.80 0.56 0.27 0.18 0.136	0. 6б8 0.644 0.682 0.678 0.663	4.19 4.19 4.19 4.19 4.19
Жидкий Оа	83 " 103 123 143	1175 1079 964 783	0. 203 0.178 0.108 0.096	0.1W 0.1341 0.1082 0.0826	1.67 1.75 1.90 2.67
НДМГ	273 283 293 303 313	783	1.135 1.096 1.065 1.034 1.002	0.2012	2.73
Теплофизические свойства ^металлов
Материал	Температура	Удельная теплоемкость ч>.с^Н	Коэффициент теплопроводности
Алюминий	273 373 473 573 673 773	0.87У 0.942 0.990 1.038 1.059 1.101	201.2 204.7 229.2 230.3 318.7 374.5
Медь	293 373 473 573 673 773 873 973 1073 1173	0.381 0.399 0.410 0.422 0.433 0.445 0.456 0.465 0.473 0.482	ЗУ5 392 382 373 363 354 344 336 329 321
Титан	2УЗ 373 473 573 673 773 873 973 1073 1173.	0.528 0.544 0.565 0.586 0.601 0.616 0.616 0.616 0.620 0.624	15.1 15.7 16.3 16.9 17.45 18.0 18.6 19.2 19.9 20.6
Никель	273 373 473 673 873 1073	0. 457 0. 470 0.486 0.516 0.540 0.558	67.5 62.8 58.2 52.3 57.3 65.1
Сталь 1Х18Н9Т	373 473 573 673 773 873 973 1073	0. ЬОе!	15. 0 17.6 19.2 20.8 22.3 23.8 25.5 27.6
Латунь Л62	698	0.387	133.7
Материал	Температура	Удельная теплоемкость Ср,Г КДж ] L КГ К J	Коэффициент теплопроводности Hrrf
Бронза	848	0.396	104.67
Молибден	згз 373 473 573 673 773 873 973 1073 1173 1273 1373 1473 1573 1673 1773 1873 1973 2073 2173 2273	О. 247 0.272 0.283 0.294 0.305 0.327 0.3397 0.3391 0.3386 0.338 0.3374 0.3368 0.3362 0.3356 0.335	134. О 138.16 133. 79 129.41 125.04 120.67 116.29 111.92 107.55 103.18 98.81 98.14 97. 47 96.8 96.13 95.46 94.79 94.12 93.45 92.78 92.11
Зольфрам	273 373 473 573 673 773 873 973 1073 1173 1273 1373 1473 1573 1673 1773 1873 1973 2073 2173 2273	0.134 0.136 0.138 0.1401 0.142 0.144 0.1456 0.1479 0.1499 0.1518 0.1557 0.1569 0.1582 0.1595 0.1612 0.163 0.1648	129. 46 127. 85 126.23 124.62 123.0 121.39 119.77 118.15 116.54 118.15 113.56 112.20 110.75 109. 31 108. 27 107.23 106.19 103.48 110.85 118.23 125.60

Зависимость «н 0 от рн 0 • I и Т
ПРИЛОЖЕНИЕ 3
Примеры реализации используемых численных методов Данные примеры взяты из [161 и адаптированы для языка MS
FORTRAN 5.0, используемого на ПЭВМ. С
СХХХХХХХХХХХХХХХХХХХХХХХХХХХХХХХХХХХХХХХХХКХХХХХХХХХХХХХХХУ
С	Метод Симпсона
Схххххххххххххххххххххххххххххххххххх-ххххххххххххххххххххх»
с EXTERNAL F COMMON Z 1 WRITEC*,*) 7 N,ZO,Zg,H?’
C
C N - КОЛИЧЕСТВО ИЗМЕНЕНИИ ИНТЕРВАЛА ИНТЕГРИРОВАНИЯ
С zo,zg - ОТРЕЗОК ИЗМЕНЕНИЯ ПАРАМЕТРА ФУНКЦИИ FCX) С Н - ШАГ ИЗМЕЕННИЯ ПАРАМЕТРА
С
READC*,*) N,ZO,Zg,H K=czg-zo)/H+i.5 z=zo DO 2 1=1,К
CALL SIMPC0.0,3.1415g265/2,N,F,S) ! МЕТОД ИНТЕГРИРОВАНИЯ WRITECx,*) Z,S
2 Z=Z+H GOTO 1 END
SUBROUTINE SIMPC A,B,N,F,S) ’ МЕТОД СИМПСОНА H=CB-A)/C2*N) S=FCA)/2
X=A DO 11 1=1,N X=X+H S=S+2*FCX) X=X+H F1=F(X)
И
S=S+F1
S=C2*S-F13*H/3.
RETURN
ENT'
FUNCTION FCX)
COMMON Z
F=SINCXj
F=SQRT(1-Z*F*Fj
RETURN
ENL
! ПОДЫНТЕГРАЛЬНАЯ ФУНКЦИЯ
Схэеннннннннннннененеееееееееееее*************************** С	Метод средних прямоугольников
Сххмхххж-мжхххххххххххххххххххмххххкхххххххххххххжхххххххххх С
EXTERNAL F
COMMON P.Z	'
1 WRITEC*.*) ’ N,P,Z0,Z9,H?’
2	N -	КОЛИЧЕСТВО ТОЧЕК
2	P -	ПАРАМЕТР ДЛЯ ПОДЫНТЕГРАЛЬНОЙ	ФУНКЦИИ PCX)
С	ZO-	НАЧАЛО ОТРЕЗКА ИНТЕГРИРОВАНИЯ
:	Z9-	КОНЕЦ ОТРЕЗКА ИНТЕГРИРОВАНИЯ
2	Н -	ШАГ ИНТЕГРИРОВАНИЯ
PEAD'*.*} N.P.Z0.Z9.H 3=3 14259255 0=1 Т
K='Z5-Z0j <4*1.5
ОС 2 1=1. К
CALL Р.ЕСТСО. .B,N,F,S) ! МЕТОД ИНТЕГРИРОВАНИЯ WRITE!x.*j Z.C*S
GOTO 1
END
SUBROUTINE RECTCA,B,N,F,S) ! МЕТОД ПРЯМОУГОЛЬНИКОВ H=CB-A)/N X=A+H/2 s=o. DO 11 1=1,N S=S+FCX)
11 X=X+H S=S*H RETURN END
FUNCTION FCX)	! ПОДЫНТЕГРАЛЬНАЯ ФУНКЦИЯ
COMMON P,Z F=COSCZxSINCX)-P*X) RETURN END
C
Cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxmxxxxxx С	МНК co степенным базисом
Схххххххххххххххххххххххххххххххмхххххххххххххххххххххххххх c REAL XC50),FC50),AC1O,113,CC1O)
1 WRITE(x,x) ’N,M,X0,X9,H?’
C
C	N -	КОЛИЧЕСТВО ТОЧЕК
С	M -	СТЕПЕНЬ ПОЛИНОМА
С	ХО-	НАЧАЛО ОТРЕЗКА ВЫДАЧИ	РЕЗУЛЬТАТА
С	Х9-	КОНЕЦ ОТРЕЗКА ВЫДАЧИ	РЕУЛЬТАТА
С	Н -	ШАГ ВЫДАЧИ РЕЗУЛЬТАТ
С
READC*,*) N,M,X0,X9,H CALL TABCN.X.F) ! ТАБЛИЦА ДАННЫХ XCN),FCN) CALL GRAMCN,M,X,F,A) ! МАТРИЦА ГРАМА ACM,М+1) CALL GAUSSCM,A,C) ! МЕТОД ГАУССА ДЛЯ СЛАУ do г 1=1,м г	WRITEC*,3) I.CCI)
3	FORMATC1X,'С',12,'=’,1РЕ13.6)
К=СХ9-Х0)/Н+1.5 Х1=ХО DO 4 1=1,К CALL FICM,C,X1,P) ! АППРОКСИМИРУЮЩАЯ ФУНКЦИЯ WRITEC*,*) XI,Р 4	Х1=Х1+Н
GOTO 1 END
SUBROUTINE TABCN,X,F) • ТАБЛИЦА ДАННЫХ XCN),FCN) REAL XC50),FC50) DO 11 1=1,N WRITEC*,12) I,I 11	READC*,*) XCD.FCI)
12	FORMATCS,’ X’,12,',F’,12,'?’)
RETURN END
SUBROUTINE GRAMCN,M,X,F,A) !МАТРИЦА ГРАМА ACM,M+1) REAL XC50),FC50),AC10,11) DO 22 J=1,M S=0. R=S Q=S DO 21 1=1,N P=XCI)**(J-1) S=S+P R=R+P«FCIJ 21	Q=Q+P«XCI)**CM-1)
A(1,B=S
ACJ,M+1)=R
22	ACJ,M)=Q
DO 23 1=2,M
DO 23 J=1,M-1
23	ACI, Л=АС1-1, J+13
RETURN
END
SUBROUTINE GAUSSCN,A,X) ! МАТРИЦА ГАУССА ДЛЯ СЛАУ
REAL АСЮ, 111,XC1O)
N1=N+1
DO 32 K=1,N
K1=K+1
S=ACK,K3
DO 31 J=K1,N1
31 ACK, Л=АСК, Л/S
DO 32 I=K1,N
R=ACI,K)
DO 32 J=K1,N1
32 ACI, Л=АС1, Л-АСК, J)*R
XCN3=ACN,N1)
DO 34 I=N-1,1,-1
S=ACI,N1)
DO 33 J=I+1,N
33	S=S-ACI, ЛмХСЛ
34	XCI)=S
RETURN
END
SUBROUTINE FICM,C,X1,P) ! АППРОКСИМИРУЮЩАЯ ФУНКЦИЯ
REAL CC1O)
P=CCM1
DO 41 I=M-1,1,-1
41 P=CCI)+P*X1
RETURN
END
XXXMXXMMXXWXMXWWKXXMXXXXXXXXXXXXXXXXXXXXXMXXXXMXMMXXXXXXXX Интерполяция сплайнами хххххххххххххххххххххххххххххххххххмхххххххххххххххххххххх
REAL XC1OO),FC1OO),CC1OO) WRITECx.x) ’N,X0,X9,H?’
N - КОЛИЧЕСТВО ТОЧЕК
X0.X9 - ОТРЕЗОК ВЫДАЧИ РЕЗУЛЬТАТА Н - ШАГ ВЫДАЧИ РЕЗУЛЬТАТА
READCx.x) N,X0,X9,H
CALL TABCN.X.F) ! ФОРМИРОВАНИЕ ТАБЛИЦЫ CALL CSCN,X,F,C) ! КОЭФФИЦИЕНТЫ СПЛАЙНОВ К=СХ9-Х0)/Н+'1.5 Х1=ХО do г i=i,к CALL SPCN,X,F,C,X1,Р,Р1,Р2) ! ВЫЧИСЛЕНИЕ СПЛАЙНОВ WRITECx.x) Х1,Р,Р1,Р2 Х1=Х1*Н GOTO 1 END
SUBROUTINE TABCN.X.F) ! ФОРМИРОВАНИЕ ТАБЛИЦЫ PEAL ХС100),FC 100) DC 11 I=1.N+1 WRITE! x, 12) I, I 1	P.EADCx.x? XCD.FCI)
2	FORMAT!®.’X’,12,’,F’,12. "7’)
RETURN END
SUBROUTINE CSCN,X,F,C) ! КОЭФФИЦИЕНТЫ СПЛАЙНОВ
REAL XC100),FC 100),CC100),КС 100)
KC2)=0.
CC2)=0.
DO 21 1=3,N+l
J=I-1
M=J-1
A=XCI)-XCJ)
B=XCJ)-XCM)
R=2*CA+B)-B*CCJ)
CCI)=A/R
21	КС I)=C3. *CCFCI)-FCJ))/A-CFCJ)-FCM))/B)-B*KCJ))/R
CCN+1)=KCN+1) DO 22 I=N,3,-1
22 CCI)=KCI)-CCI)*CCI+1)
RETURN
END'
SUBROUTINE SPCN,X,F,C,X1,P,P1,P2) ! ВЫЧИСЛЕНИЕ СПЛАЙНОВ REAL XC100),FC100),CC100) 1=2
31 IFCX1.LE.XCD) GOTO 32 1=1+1
IFCI.NE. N+l) GOTO 31
32	J=I-1
A=FCJ)
B=XCJ)
Q=XCI)-B
R=X1-B
P=CCI)
D=CCI+1)
B=CFCI)-A)/Q-CD+2*P)*Q/3.0
D=CD-P)/Q*R
P1=B+R*C2*P+D)
P2=2mCP+D)
P=A+R*CB+R*CP+D/3.0))
RETURN
END
Метод секущих
EXTERNAL F
WRITEC*,*) ’ XO.X.E?’
XO, X - НАЧАЛЬНЫЕ ПРИБЛИЖЕНИЯ КОРНЯ E - ТРЕБУЕМАЯ ТОЧНОСТЬ
READC*,*) ХО.Х.Е
CALL SECANT CXO,X,E,F)
WRITEC*,*) ' X=’,X
GOTO 1
FORMAT C’X=’,G13.6)
END
SUBROUTINE SECANT CXO,X,E,F)	! МЕТОД СЕКУЩИХ
P=X-XO
:=fcxc:
X=X*R
IFCABSCR). GT. E) GOTO 11
RETURN
END
FUNCTION FCX)	! РАССМАТРИВАЕМАЯ ФУНКЦИЯ
F=C0SC3.1415926*X*X)
RETURN
END'
Булыгин Юрий Александрович
Кретинин Александр Валентинович Рачук Владимир Сергеевич Фалеев Сергей Владиславович
РАСЧЕТ ТЕПЛОВОГО СОСТОЯНИЯ КАМЕРЫ ЖИДКОСТНОГО РАКЕТНОГО ДВИГАТЕЛЯ
Учебное пособие
Редактор’ А. Д. Гуляева
ЛР N020419 от 12.02.92. Подписано в печать 5.02.97.
Формат 60x84/16. Бумага для множительных аппаратов.
Усл. печ. л. 5,6. Уч.-изд. л. 5,4. Тираж 100 экз. Заказ N
Воронежский государственный технический университет 394026 Воронеж, Московский проспект, 14
Участок оперативной полиграфии Воронежского государственного технического университета