Identification and simulation mathematical model of thermo and physical loading of a small-caliber artillery barrel

全文:

Введение

Рост требований к интенсификации режимов стрельбы (увеличению количества выстрелов в очереди и количества очередей выстрелов, уменьшению интервалов времени между очередями выстрелов), непрерывность разработок по созданию новых скорострельных артиллерийских орудий, обусловливают исключительную важность корректировок методов и способов описания процесса нагрева ствола, ограничивающего боевые свойства ААО [1]. Поэтому научно-техническая задача формализации температурного поля ствола, циклично подвергающегося высоким термомеханическим нагрузкам, представляется приоритетной задачей проектирования и исследования ААО.

При решении введенной в рассмотрение проблемы необходимо определить, во-первых, в какой мере будут действенны традиционные методы решения обобщенной задачи нагрева инженерных конструкций и, во-вторых, какие новые подходы можно предложить для достижения цели, поставленной в работе. Так, в статьях [2–4] предложены схемы экспериментальных исследований и методы обработки выходных данных, обеспечивающие повышение точности определения температуры тела и расширение диапазона измерений; в публикации [5] представлен уникальный вариант тепловой модели, разработанной на основе аппарата теории вероятностей; в работе [6] численными решениями многомерной задачи теплопроводности определены температурные поля оребренных стенок различной конфигурации; трудом [7] предложены инструменты моделирования температурного поля в узлах газовых турбин, максимально учитывающие совокупность параметров в многофакторных граничных условиях пограничного слоя; в статье [8] получены корреляционные регрессионные зависимости оптимальных экстремумов нагружения стволов стрелково-пушечного артиллерийского вооружения. В качестве примеров работ подобной тематики в области авиационной артиллерийской науки можно привести публикации [9–13].

В анализируемых подходах находит применение как простой метод перебора различных вариантов и их комбинаций, так и группа сложных методов класса «предиктор – корректор», предусматривающих экстраполирующие и исправляющие формулы. Естественно, указать какую-либо определенную, строго регламентированную методику, дающую истинные результаты по нагреву тел разнообразных конструкций и условий эксплуатации нельзя. Тем не менее поиск решений, адекватно сопоставляющих сложность объекта со сложностью математических или экспериментальных средств (закон адекватности объекта и регулятора) актуален в ближайшей и долгосрочной перспективе. По принятой терминологии [14; 15] непосредственное построение отдельных блоков модели подразумевает процесс имитации. В случае использования известных заимствований и установлении соответствующей организации параметров, в работе различима проблема идентификации модели.

Синтез дифференциальных операторов идентификации и имитации температурного поля ствола

При решении задач теплопроводности используется основное дифференциальное уравнение теплопроводности Фурье, в общем случае описывающее распределение температуры в твердом теле [16]. Для ствола, имеющего сложную геометрическую форму составного цилиндра, уравнение теплопроводности Фурье целесообразно представить в цилиндрической системе координат, а так как ствол является телом вращения и симметричен относительно продольной оси, рассмотрение задачи определения температурного поля ствола ограничено двумерной формой в одной из радиальных плоскостей [16; 17]:

$\frac{\partial T}{\partial \text{τ}}=a\left(\frac{{\partial }^{2}T}{\partial r^2}+\frac{1}{r}\frac{\partial T}{\partial r}+\frac{1}{r}\frac{{\partial }^{2}T}{\partial {\text{θ}}^2}\right),$                                                                  (1)

где T – температура ствола; τ – время; a – коэффициент температуропроводности стали ствола; r, θ – радиус вектор и полярный угол, соответственно, цилиндрической системы координат.

Коэффициент температуропроводности материала в уравнении (1), характеризующий скорость изменения температуры тела, особо важен для описания именно нестационарных тепловых процессов и рассчитывается по формуле:

$a=\frac{\text{λ}}{с\text{ρ}},$                                                                                                            (2)

где λ, с, ρ – коэффициенты теплопроводности, удельной теплоемкости и плотность, соответственно, стали ствола.

Реальное явление нагрева и охлаждения ствола абстрагируется некоторыми допущениями, упрощающими, но не снижающими качества построения модели и основанными на обобщении опыта, аккумулированного в прикладных исследованиях, отмеченных ранее:

Принятие крайнего допущения при исследовании вопроса о безопасном нахождении очередного патрона в разогретом стрельбой стволе, диктует объективную необходимость прямого учета в формуле (2) зависимости коэффициентов теплопроводности λ и удельной теплоемкости c ствольной стали от температуры ствола T путем процедурных запросов требуемых значений из массива данных [18]. Вместе с тем нелинейность теплофизических коэффициентов материала ствола, несомненно, приводит к существенному усложнению задачи и увеличению времени машинного счета. Поскольку цилиндр является осесимметричным, то из соображений целесообразности рассмотрено его продольное сечение. Исходя из первого допущения, в целях дифференциации расчетных процедур, при известной связи цилиндрических координат с прямоугольными декартовыми координатами [17], уравнение (1) заменено двумя одномерными дифференциальными уравнениями, одно из которых записано в цилиндрических, а другое – в прямоугольных декартовых координатах. В результате получена система уравнений вида:

$\left.\begin{array}{l}\frac{\partial T}{\partial \text{τ}}=a\left(\frac{{\partial }^{2}T}{\partial r^2}+\frac{1}{r}\frac{\partial T}{\partial r}\right);\\ \frac{\partial T}{\partial \text{τ}}=a\frac{{\partial }^{2}T}{\partial x^2},\end{array}\right\}$                                                                                    (3)

где x – ось абсцисс прямоугольной декартовой системы координат.

Дуализм условий в формуле (3) объясняется потребностью в согласовании точностных характеристик применяемых физико-математических методов с располагаемыми характеристиками машинного времени, памяти и разрядной сетки вычислителя. Такой прием предполагает последовательное подключение подпрограмм-процедур, содержащих вызовы более упрощенных, по сравнению с полной, математических моделей, которые описывают локальные изменения температурного поля ствола в поперечном r и продольном x направлениях.

В приложении к моделированию теплофизического нагружения ствола первое уравнение системы (3) определяет нестационарную температуру ствола T исключительно в радиальном направлении r. Второе уравнение системы (3) определяет эту же величину в продольном направлении по оси x. В целом система уравнений (3) характеризует пространственно-временнóе изменение температуры цилиндра в любой точке, объединяя все типы явления теплопроводности независимо от физических свойств тела и условий взаимодействия с окружающей средой. Для выделения изучаемого явления к системе уравнений (3) присоединены условия однозначности, специфичные для объекта и предмета исследования. В условия однозначности включены требуемые интерпретации геометрических, краевых (совокупность начального и граничных), физических условий однозначности. На рис. 1 представлена схема осевой симметрии ствола наиболее распространенной авиационной пушки ГШ-301.

 

Рис. 1. Схема осевой симметрии ствола авиационной пушки ГШ-301

Fig. 1. Scheme of axial symmetry of the GSh-301 aircraft gun barrel

 

При самой большой для 30-миллиметровых пушек скорострельности на один ствол, данное ААО с типовым осколочно-фугасно-зажигательным снарядом ОФЗ-30 обеспечивает режимы стрельбы с максимальной тепловой нагрузкой ствола. Ось x на рис. 1 совмещена с продольной осью ствола, и рассмотрение области, расположенной только по одну сторону от продольной оси ствола, вполне достаточно.

Начальные условия однозначности определяются заданным распределением температур в стволе для фиксированного времени, предшествующему рассматриваемому и принимаемому за начальный момент времени τ = 0. Уравнение температурного поля ствола для этого момента времени τ записано в виде:

 T(r, x, 0) = T₀(r, x).                                                                                                 (4)

Так как начальная температура ствола T₀ в любой точке (r, x) в момент времени τ = 0 равна температуре окружающей среды, уравнение (4) имеет вид:

 T(r, x, 0) = T_в = const.                                                                                             (5)

Граничные условия определяют значения параметров теплопередачи на границах ствола. Ствол авиационной пушки ГШ-301 является нетеплоизолированным, поэтому граничные условия заданы в виде температур окружающей среды и законов теплообмена между средой и поверхностями ствола, в зависимости от конструктивных параметров и условий применения ААО (условия III-го рода) [16].

Интенсивность теплообмена между пороховыми газами (далее – газы), образующимися в процессе горения пороха, стволом и воздухом зависит от физико-механических процессов, протекающих у границ ствола. Эти процессы достаточно точно описываются известными уравнениями основной задачи внутренней баллистики для автоматических артиллерийских орудий [1; 19; 20] и уравнением теплового баланса [16].

Для каждой точки внутренней границы ствола граничное условие, в котором входящие величины зависят от времени τ и координаты x, записано через уравнение теплового баланса:

${\text{α}}_{\text{г}}\left[T_{\text{г}}\left(\text{τ}\right)-T\left(x\text{, τ}\right)\right]=-\text{λ}\frac{\partial T}{\partial r},$                                                                                         (6)

где α_г – коэффициент теплоотдачи от газов к каналу ствола; T_г – температура газов в канале ствола.

Уравнением (6) определяется конвективный тепловой поток от газов, направленный вглубь стенки ствола к его внешней поверхности в данной точке.

Воздействие движущегося снаряда на канал приводит к дополнительному нагреву ствола. При использовании в автоматической стрельбе медных ведущих поясков средняя температура трения рассчитывается по эмпирической формуле, учитывающей текущую по координатам и времени скорость газов в канале ствола v [20].

Граничное условие для каждой точки внешней границы ствола, на которой действует конвективный тепловой поток от воздуха, образовано аналогично уравнению (6)

${\text{α}}_{\text{в}}\left[T_{\text{в}}-T\left(x\text{, τ}\right)\right]=-\text{λ}\frac{\partial T}{\partial r},$                                                                                              (7)

где α_в – коэффициент теплоотдачи от внешней поверхности ствола к воздуху.

На дульном срезе ствола действует конвективный тепловой поток, формируемый температурой газов у дульного среза T_д. В этом случае уравнения (6) и (7) записаны для окончания процесса выстрела

${\text{α}}_{\text{г}}\left[T_{\text{д}}\left(\text{τ}\right)-T\left(r\text{, τ}\right)\right]=-\text{λ}\frac{\partial T}{\partial х}.$                                                                                             (8)

Для интервалов времени между очередями выстрелов уравнения (7) и (8) имеют вид:

${\text{α}}_{\text{в}}\left[T_{\text{в}}-T\left(r\text{, τ}\right)\right]=-\text{λ}\frac{\partial T}{\partial х}.$                                                                                                (9)

Так как канал ствола является замкнутым, т. е. непосредственное распространение теплового излучения во внешнюю среду отсутствует, а величина потока теплового излучения от внешней стенки ствола в окружающую среду незначительна и по предварительным расчетам составляет менее 1 % от температуры ствола T на его внешней поверхности, то расчеты теплового потока лучистой энергии в работе не проводятся.

При очередях выстрелов канал и внешняя поверхность ствола имеют достаточно высокие температуры, поэтому необходим учет особенностей конструкции опорного образца ААО. Моделирование процесса функционирования штатной системы охлаждения авиационной пушки ГШ-301 реализовано введением локального коэффициента теплоотдачи.

Учет влияния гильзы, находящейся во время выстрела в казенной части, на распределение температурного поля ствола необходимени Поскольку гильза, располагаемая на расстоянииx = 0,175 м длины ствола, под действием давления газов плотно прижимается к стенке патронника, можно принять допущение о том, что такой контакт является идеальным, т. е. перенос тепла осуществляется лишь посредством теплопроводности. Так как толщина стенки гильзы сравнительно мала, предположен ее мгновенный прогрев до температуры газов в канале ствола TНа основании вышеуказанного сформулировано граничное условие на поверхности патронника в характерных точках ствола, где происходит непосредственный контакт гильзы со стенкой (условие I-го рода) [16]

$T(r=\mathrm{0,}x=\mathrm{0...}\mathrm{0,175})=T_{\text{г}}\left(\text{τ}\right),$                                                                                          (10)

Физические процессы конвективного переноса тепла при автоматической стрельбе и в перерывах между очередями выстрелов носят достаточно сложный характер. Совершенный математический учет теплообмена на стенках ствола посредством коэффициентов теплоотдачи от газов к каналу ствола α₁ и от внешней поверхности ствола к воздуху α₂ представляется самостоятельной задачей, базирующейся, как правило, на применении аппарата теории подобия процессов теплообмена в артиллерийских орудиях [19; 20].

Таким образом, при заданных дифференциальных уравнениях процесса (3), геометрических условиях однозначности, показанных на рис. 1, физических условиях однозначности, двумерная, нестационарная, нелинейная задача по определению температурного поля ствола представлена системой уравнений, полностью определяющей краевую задачу (4)–(10):

$\left.\begin{array}{l}\frac{\partial T}{\partial \text{τ}}=a\left(\frac{{\partial }^{2}T}{\partial r^2}+\frac{1}{r}\frac{\partial T}{\partial r}\right);\\ \frac{\partial T}{\partial \text{τ}}=a\frac{{\partial }^{2}T}{\partial x^2};\\ {\text{α}}_{\text{г}}\left[T_{\text{г}}\left(\text{τ}\right)-T\left(x\text{, τ}\right)\right]=-\text{λ}\frac{\partial T}{\partial r};\\ {\text{α}}_{\text{в}}\left[T_{\text{в}}-T\left(x\text{, τ}\right)\right]=-\text{λ}\frac{\partial T}{\partial r};\\ {\text{α}}_{\text{г}}\left[T_{\text{д}}\left(\text{τ}\right)-T\left(r\text{, τ}\right)\right]=-\text{λ}\frac{\partial T}{\partial х};\\ {\text{α}}_{\text{в}}\left[T_{\text{в}}-T\left(r\text{, τ}\right)\right]=-\text{λ}\frac{\partial T}{\partial х};\\ T(r=\mathrm{0,}x=\mathrm{0...}\mathrm{0,175})=T_{\text{г}}\left(\text{τ}\right);\\ T(r,x\mathrm{,0})=T_{\text{в}}=\text{const.}\end{array}\right\}$                                                                                             (11)

Достоинства системы уравнений (11) заключаются в простоте группирования исходных данных и наглядности их распределения по уровням иерархии с постоянно нарастающей степенью детализации.

Тем не менее система уравнений (11) не может быть непосредственно применена для решения поставленной задачи. Дело в том, что аналитическое решение задач теплопроводности и теплоотдачи возможно для тел простых форм и при несложных граничных условиях [21; 22]. Аналитические методы решения задач теплопередачи подробно изложены в работе [23]. Сложная двумерная, нестационарная, нелинейная прикладная задача теплопередачи в области со сложной конфигурацией границ в такой постановке не имеет точного аналитического решения. Получение результата, близкого к вероятно существующему аналитическому решению, возможно лишь селекцией численных методов, суть которых сводится к конечно-разностным аппроксимациям системы уравнений (11).

Формирование разностной схемы и селекция методов расчета температурного поля ствола

Приближенное решение задачи теплообмена в стволе проведем применением метода конечных разностей как наиболее апробированного и приемлемого для тел с переменными по длине сечениями [21–24]. Область непрерывного изменения аргумента заменена разностной сеткой (далее – сетка) – дискретным множеством точек, пересечения которых образует узлы. Вместо функций непрерывных аргументов введены функции дискретных аргументов – сеточные функции, определяемые в узлах сетки. Частные производные, входящие в инегродифференциальные уравнения и краевые условия, заменены (аппроксимированы) разностными соотношениями. В результате подобных замен, систему дифференциальных уравнений в частных производных (11) требуется свести к системе конечно-разностных алгебраических (далее – разностных) уравнений – разностным схемам. Если решение системы разностных уравнений существует и при увеличении степени подробности сетки стремится к решению поставленной задачи, то это решение и будет являться искомым приближенным решением аналитической задачи. Несмотря на то, что число неизвестных в системе разностных уравнений будет весьма значительно, ее решение с точки зрения математических трудностей проще, чем решения исходной системы дифференциальных уравнений (11). Следовательно, при решении системы (11) необходимо: осуществить выбор конфигурации и размеров сетки; построить разностную схему; определить устойчивость разностной схемы; выполнить многошаговую аппроксимацию исходной дифференциальной задачи; выяснить сходимость разностной схемы.

Итак, заменим область Ω_T непрерывного изменения аргументов искомой величины T некоторым конечным множеством точек, лежащих в этой области. Точки сетки формирования конечной разности функции целочисленного аргумента T_kj по оси х будем обозначать через k, а такие же точки по оси r – через j. В соответствии со спецификой решаемой задачи, область Ω_T трансформируется в область W_Tkj  расчета температуры T_kj в k_j-точках сечений ствола при величине τ. Решение исходной задачи сводится к нахождению таблицы числовых значений функции T_kj в узлах сетки области W_Tkj.

В соответствии с выбранной системой координат (0, x, r,) разобьем ствол в направлении оси x на ϑ равных частей ϑ = l / Δx, а в направлении оси r на ν равных частей ν = r_y / Δr, где Δx, Δr – шаги сетки по соответствующим координатам; l – длина ствола; r_y – максимальная толщина ствола. Для этого проведем ϑ – 1 лучей в направлении, перпендикулярном оси x и ν – 1 лучей в направлении, перпендикулярном оси r, как представлено на рис. 2. В результате такого разбиения будем иметь сетку, состоящую из совокупности внутренних (на рис. 2 обозначены ·  ) и граничных (на рис. 2 обозначены  ) узлов. Так как в рассматриваемом случае Δx = l / ϑ = const и Δr = r_y / ν = const, то множество узлов x_k, определяемых точками с номерами k = 0, 1, 2, …,K_ϑ, и множество узлов r_j, определяемых точками с номерами j = 0, 1, 2, …, J_ν, представляет собой равномерную пространственную сетку в области W_Tkj.

 

Рис. 2. Сеточная схема ствола авиационной пушки ГШ-301

Fig. 2. Grid diagram of the GSh-301 aircraft gun barrel

 

Здесь следует остановиться на одной из особенностей построения сетки. Эта особенность заключается в том, что разбиение осей х и r на отрезки производится не от начала координат, а от поверхности ствола с учетом конфигурации его продольного сечения и нарезов. Такое разбиение делается для того, чтобы в процессе составления разностных уравнений граничные узлы максимально совпадали с положением поверхности области W_Tkj. При этом, конечно, не все граничные узлы окажутся на линии, определяющей границу поверхности ствола, а расстояние между r_j–1 и r_j узлами не на всех ν лучах может оказаться равным Δr. Однако эти обстоятельства менее существенны по сравнению с выполнением граничных условий III-го рода, так как они в основном и определяют характер описания процесса его нагрева и охлаждения.

При построении разностных аналогов дифференциальных операторов системы (11) воспользуемся методом формальной замены производных конечно-разностными отношениями. Наиболее естественный способ замены производной основывается на определении производной (например, по координате r 1-го уравнения системы (11)) как предела [17; 24]

$\frac{\partial T}{\partial r}=\underset{\Delta r\to 0}{\lim }\left[T(r+\Delta r)-T\left(r\right)\right]\frac{1}{\Delta r}.$                                                                               (12)

Если в равенстве (12) зафиксировать шаг Δr, то получим приближенную формулу для первой производной, выраженной через конечные разности.

Для так называемого правого разностного отношения или разности «вперед»

$\frac{\partial T}{\partial r}\approx \left[T(r+\Delta r)-T\left(r\right)\right]\frac{1}{\Delta r}.$                                                                                       (13)

Аналогично вводится левое разностное отношение (разность «назад»), записываемое в виде:

$\frac{\partial T}{\partial r}\approx \left[T\left(r\right)-T(r-\Delta r)\right]\frac{1}{\Delta r}.$                                                                                       (14)

При решении задач теплопроводности необходимо аппроксимировать и вторую производную. Для второй производной рассматривается линейная комбинация отношений (13) и (14)

$\frac{{\partial }^{2}T}{\partial r^2}\approx \left[T(r+\Delta r)-2T\left(r\right)+T(r-\Delta r)\right]\frac{1}{\Delta r^2}.$                                                                  (15)

Каждый переход на 1 шаг «вперед» условно обозначен через «+1», а «назад» – через «–1». Тогда для j-й точки сетки формирования конечной разности величины T_kj по оси r правое разностное соотношение (13) преобразуется к виду

$\frac{\partial T}{\partial r}=(T_{j+1}-T_j)\frac{1}{\Delta r}.$                                                                                             (16)

Аналогично преобразовано левое разностное отношения (14)

$\frac{\partial T}{\partial r}=(T_j-T_{j-1})\frac{1}{\Delta r}.$                                                                                             (17)

Разностный аналог второй производной, соответствующий формуле (15), представлен отношением:

$\frac{{\partial }^{2}T}{\partial r^2}=(T_{j+1}-2T_j+T_{j-1})\frac{1}{\Delta r^2}.$                                                                                (18)

Формулы (12) – (18) и их обоснования справедливы и при замене разностными отношениями производной ∂²T/∂x² во втором уравнении системы (11). В этом случае в уравнениях-аналогах вместо переменной r будет присутствовать переменная x, а индекс j будет замещен индексом k. Показанные аналогии будем иметь в виду далее, иногда не прибегая к прямой детализации задачи по второй пространственной переменной x.

Для построения соотношений, аппроксимирующих временнýю производную ∂T/∂τ в первом и втором уравнениях системы (11), в принципе можно использовать значения температур в kj-точках сечений ствола в различные моменты времени: T_{k, j, i}, T_{k, j, i–1}, T_{k, j, i–2}, … , где i – точка сетки формирования конечной разности величины температуры T_kj в kj-точках сечений ствола по времени τ. Однако в практике решения большинства задач теплопроводности в подавляющем большинстве случаев используются исключительно двухслойные (по времени τ) разностные схемы, аппроксимирующие значения искомых температур на текущей i-й и предыдущей (i – 1)-й временнóй точке. Значительно реже учитываются значения температур в (i – 2)-й момент времени получением трехслойных разностных схем [21–24].

При получении вариантов двухслойных разностных схем, производная по времени τ аппроксимирована временнóй разностью «назад»

$\frac{\partial T}{\partial \text{τ}}=(T_i-T_{i-1})\frac{1}{\Delta \text{τ}},$                                                                                             (19)

где Δτ – шаг сетки по времени τ.

Пространственные дифференциальные операторы при двухслойной разностной схеме аппроксимированы также на основе значений температуры T_kj в kj-точках сечений ствола в i-й и (i – 1)-й моменты времени τ. При этом возможны два предельных случая.

В первом случае в аппроксимации задействованы только значения температуры T_kj в kj-точках сечений ствола для текущего i-го момента времени τ. Так, для пространственной переменной r одномерная пространственно-временнáя аппроксимация соответствующего дифференциального оператора уравнения (11) будет иметь вид:

$\frac{{\partial }^{2}T}{\partial r^2}=(T_{j+\mathrm{1,}i}-2T_{j,i}+T_{j-\mathrm{1,}i})\frac{1}{\Delta r^2}.$                                                                             (20)

Во втором случае при аппроксимации используются только значения температуры T_kj в kj-точках сечений ствола для предыдущего момента времени (i – 1)

$\frac{{\partial }^{2}T}{\partial r^2}=(T_{j-\mathrm{1,}i-1}-2T_{j,i-1}+T_{j-\mathrm{1,}i-1})\frac{1}{\Delta r^2}\text{\hspace{0.17em}.}$                                                                  (21)

Принимая во внимание формулу (2), в соответствии с уравнениями (19) – (21) представлено два различных разностных уравнения, аппроксимирующих первое уравнение системы (11) в одномерной версии:

$\frac{T_{j,i}-T_{j,i-1}}{\Delta \text{τ}}=\frac{\text{λ}}{с\text{ρ}}(T_{j+\mathrm{1,}i}-2T_{j,i}+T_{j-\mathrm{1,}i})\frac{1}{\Delta r^2};$                                                               (22)

$\frac{T_{j,i}-T_{j,i-1}}{\Delta \text{τ}}=\frac{\text{λ}}{с\text{ρ}}(T_{j-\mathrm{1,}i-1}-2T_{j,i-1}+T_{j-\mathrm{1,}i-1})\frac{1}{\Delta r^2}.$                                                    (23)

Разностное уравнение вида (23) позволяет выразить решение задачи теплопроводности в стволе в явном виде на i-м временнóм слое через известные решения на предыдущем (i – 1)-м слое. В совокупности с аппроксимацией условий однозначности (4)–(10), разностное уравнение (23) образует явную разностную схему. Алгоритмы численного расчета по явной разностной схеме достаточно просты в программировании, но предъявляют требования к машинному времени.

Разностная схема, заданная разностным уравнением вида (22), сложнее, поскольку в каждое разностное уравнение вида (22) кроме неизвестного решения для j-й пространственной точки входят еще два искомых решения для соседних (j – 1)-й и (j + 1)-й пространственных точек. Все искомые решения оказываются «завязанными» друг с другом в общую не распадающуюся систему разностных уравнений. Таким образом, в данном случае на каждом i-м временнóм слое решения определяются не по явным формулам вида (23), а из решения системы разностных уравнений, поэтому разностная схема, заданная разностным уравнением вида (22) неявна. Эффективные алгоритмы решения системы уравнений (11) по неявной разностной схеме значительно сложнее алгоритмов численного расчета по явной разностной схеме, но время решения задачи может быть существенно уменьшено рациональным выбором шагов Δx, Δr и Δτ.

Для построения конечно-разностного аналога системы уравнений (11) в работе использована исключительно неявная разностная схема в силу ее безусловной устойчивости, т. е. способности обеспечивать точность решения при любой степени подробности сетки.

При решении двумерной в пространственных координатах задачи достаточно действенным оказался локально-одномерный метод, относящийся к группе методов расщепления [21; 23] Сложная краевая задача математической физики сведена к последовательному решению одномерных задач. При решении двумерной задачи по определению температурного поля ствола локально-одномерный метод позволил в два этапа вычислить температуру Т в множестве точек ствола в фиксированное время τ путем последовательного решения двух одномерных задач. Шаг сетки по времени t при этом на каждом этапе решения поделен пополам – 0,5 Δτ.

В развитие приведенных обоснований, с учетом выражения (22), система уравнений (11) примет вид:

$\begin{array}{l}\frac{T_{k,j,i}-T_{k,j,i-1/2}}{\mathrm{0,5}\Delta \text{τ}}=\frac{\text{λ}}{с\text{ρ}}\left[\frac{T_{k,j+\mathrm{1,}i}-2T_{k,j,i}+T_{k,j-\mathrm{1,}i}}{\Delta r^2}+\frac{1}{r}\left(\frac{T_{k,j,i}-T_{k,j,i-1}}{\Delta r}\right)\right];\\ \frac{T_{k,j,i}-T_{k,j,i-1/2}}{\mathrm{0,5}\Delta \text{τ}}=\frac{\text{λ}}{с\text{ρ}}\left(\frac{T_{k+\mathrm{1,}j,i}-2T_{k,j,i}+T_{k-\mathrm{1,}j,i}}{\Delta x^2}\right);\\ {\text{α}}_{\text{г}}\left[T_{\text{г}}\left(\text{τ}\right)-T_{k,j-\mathrm{1,}i}\left(x\text{, τ}\right)\right]=-\text{λ}\frac{T_{k,j,i}\left(x\text{, τ}\right)-T_{k,j-\mathrm{1,}i}\left(x\text{, τ}\right)}{\Delta r};\\ {\text{α}}_{\text{в}}\left[T_{\text{в}}-T_{k,j,i}\left(x\text{, τ}\right)\right]=-\text{λ}\frac{T_{k,j,i}\left(x\text{, τ}\right)-T_{k,j-\mathrm{1,}i}\left(x\text{, τ}\right)}{\Delta r};\\ {\text{α}}_{\text{г}}\left[T_{\text{д}}\left(\text{τ}\right)-T_{k-\mathrm{1,}j,i}\left(r\text{, τ}\right)\right]=-\text{λ}\frac{T_{k,j,i}\left(r\text{, τ}\right)-T_{k-\mathrm{1,}j,i}\left(r\text{, τ}\right)}{\Delta х};\\ {\text{α}}_{\text{в}}\left[T_{\text{в}}-T_{k-\mathrm{1,}j,i}\left(r\text{, τ}\right)\right]=-\text{λ}\frac{T_{k,j,i}\left(r\text{, τ}\right)-T_{k-\mathrm{1,}j,i}\left(r\text{, τ}\right)}{\Delta х};\\ T(r=\mathrm{0,}x=\mathrm{0...}\mathrm{0,175})=T_{\text{г}}\left(\text{τ}\right);\\ T(r,x\mathrm{,0})=T_{\text{в}}=\text{const.}\end{array}$                 (24)

Разрешением разностных уравнений системы (24) относительно T_k,j,i–1/2, T_k,j–1,i (x,τ), T_k,j,i (x,τ), T_k–1,j,i (r,τ), T_k,j,i (r,τ) получена система вида:

            $\left.\begin{array}{l}{T}_{k,j,i-1/2}=-\left(\frac{\mathrm{0,5}\text{\hspace{0.17em}λ\hspace{0.05em}}\Delta \text{τ}}{с\text{ρ\hspace{0.05em}}\Delta r^2}\right)T_{k,j+\mathrm{1,}i}+\left(\frac{\text{λ\hspace{0.05em}}\Delta \text{τ}}{с\text{ρ\hspace{0.05em}}\Delta r^2}-\frac{\mathrm{0,5}\text{\hspace{0.17em}\hspace{0.05em}λ}\Delta \text{τ}}{с\text{ρ(}{r}_{\text{0}}\Delta r+\Delta r^2{J}_{\text{ν}}\text{)}}+1\right)T_{k,j,i}-\\ -\left(\frac{\text{0,5λ\hspace{0.05em}}\Delta \text{τ}}{с\text{ρ\hspace{0.05em}}\Delta r^2}-\frac{\mathrm{0,5}\text{\hspace{0.17em}\hspace{0.05em}λ\hspace{0.05em}}\Delta \text{τ}}{с\text{ρ(}{r}_{\text{0}}\Delta r+\Delta r^2{J}_{\text{ν}}\text{)}}\right)T_{k,j-\mathrm{1,}i};\\ T_{k,j,i-1/2}=-\left(\frac{\mathrm{0,5}\text{\hspace{0.17em}λ\hspace{0.05em}}\Delta \text{τ}}{с\text{ρ\hspace{0.05em}}\Delta x^2}\right)T_{k+\mathrm{1,}j,i}+\left(1+\frac{\text{λ\hspace{0.05em}}\Delta \text{τ}}{с\text{ρ\hspace{0.05em}}\Delta x^2}\right)T_{k,j,i}-\left(\frac{\text{0,5λ\hspace{0.05em}}\Delta \text{τ}}{с\text{ρ\hspace{0.05em}}\Delta x^2}\right)T_{k-\mathrm{1,}j,i};\\ T_{k,j-\mathrm{1,}i}\left(x\text{, τ}\right)=T_{k,j,i}\left(x\text{, τ}\right)\left(\frac{\text{λ}}{\text{λ}+{\text{α}}_{\text{г}}\Delta r}\right)+\left(\frac{{\text{α}}_{\text{г}}{T}_{\text{г}}\left(\text{τ}\right)\Delta r}{\text{λ}+{\text{α}}_{\text{г}}\Delta r}\right);\\ T_{k,j,i}\left(x\text{, τ}\right)=T_{k,j-\mathrm{1,}i}\left(x\text{, τ}\right)\left(\frac{\text{λ}}{\text{λ}-{\text{α}}_{\text{в}}\Delta r}\right)-\left(\frac{{\text{α}}_{\text{в}}{T}_{\text{в}}\left(\text{τ}\right)\Delta r}{\text{λ}-{\text{α}}_{\text{в}}\Delta r}\right);\\ T_{k-\mathrm{1,}j,i}\left(r\text{, τ}\right)=T_{k,j,i}\left(x\text{, τ}\right)\left(\frac{\text{λ}}{\text{λ}+{\text{α}}_{\text{г}}\Delta x}\right)+\left(\frac{{\text{α}}_{\text{г}}{T}_{\text{д}}\left(\text{τ}\right)\Delta x}{\text{λ}+{\text{α}}_{\text{г}}\Delta x}\right);\\ T_{k,j,i}\left(r\text{, τ}\right)=T_{k-\mathrm{1,}j,i}\left(r\text{, τ}\right)\left(\frac{\text{λ}}{\text{λ}-{\text{α}}_{\text{в}}\Delta x}\right)-\left(\frac{{\text{α}}_{\text{в}}{T}_{\text{в}}\left(\text{τ}\right)\Delta x}{\text{λ}-{\text{α}}_{\text{в}}\Delta x}\right);\\ T(r=\mathrm{0,}x=\mathrm{0...}\mathrm{0,175})=T_{\text{г}}\left(\text{τ}\right);T(r,x\mathrm{,0})=T_{\text{в}}=\text{const,}\end{array}\right\}$                 (25)

где r₀ – расстояние от продольной оси ствола до его внутренней стенки.

Для решения системы разностных уравнений (25) можно применять обычные методы линейной алгебры или методы итераций. Но реализация локально-одномерного метода расщепления позволяет организовать вычисления методом прогонки или методом факторизации [21–23], наиболее выгодным и экономичным по объему вычислений.

Обозначением выражений в скобках через коэффициенты метода прогонки при искомой температуре T_kj на текущем i-м временнóм слое, система разностных уравнений (25) преобразована к компактной форме

$\left.\begin{array}{l}{T}_{k,j,i-1/2}=-A_j{T}_{k,j+\mathrm{1,}i}+C_j{T}_{k,j,i}-B_j{T}_{k,j-\mathrm{1,}i};\\ T_{k,j,i-1/2}=-A_k{T}_{k+\mathrm{1,}j,i}+C_k{T}_{k,j,i}-B_k{T}_{k-\mathrm{1,}j,i};\\ T_{k,j-\mathrm{1,}i}\left(x\text{, τ}\right)=T_{k,j,i}\left(x\text{, τ}\right){\text{ξ}}_{\text{г}}^r+{\phi }_{\text{г}}^r;\\ T_{k,j,i}\left(x\text{, τ}\right)=T_{k,j-\mathrm{1,}i}\left(x\text{, τ}\right){\text{ξ}}_{\text{в}}^r-{\phi }_{\text{в}}^r;\\ T_{k-\mathrm{1,}j,i}\left(r\text{, τ}\right)=T_{k,j,i}\left(x\text{, τ}\right){\text{ξ}}_{\text{г}}^x+{\phi }_{\text{г}}^x;\\ T_{k,j,i}\left(r\text{, τ}\right)=T_{k-\mathrm{1,}j,i}\left(r\text{, τ}\right){\xi }_{\text{в}}^x-{\phi }_{\text{в}}^x;\\ T(r=\mathrm{0,}x=\mathrm{0...0,175})=T_{\text{г}}\left(\text{τ}\right);\\ T(r,x\mathrm{,0})=T_{\text{в}}=\text{const.}\end{array}\right\}$                                                          (26)

Идентичность в системе разностных уравнений (26) группировок параметров при коэффициентах метода прогонки значительно упрощает составление машинной программы (далее – программа) расчета температурного поля ствола, учитывающей специальный вид матриц прогонки первых двух разностных уравнений системы (26) – их трехдиагональность. Так, например, для первого разностного уравнения системы (26) матрица прогонки имеет вид

$\left|\begin{array}{ccccccccccc}1& {\text{ξ}}_{\text{г}}^r& 0& \mathrm{...}& 0& 0& 0& \mathrm{...}& 0& 0& 0\\ A_1& C_1& B_1& \mathrm{...}& 0& 0& 0& \mathrm{...}& 0& 0& 0\\ \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}\\ 0& 0& 0& \mathrm{...}& A_j& C_j& B_j& \mathrm{...}& 0& 0& 0\\ \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}\\ 0& 0& 0& \mathrm{...}& 0& 0& 0& \mathrm{...}& A_{r_j-1}& C_{r_j-1}& B_{r_j-1}\\ 0& 0& 0& \mathrm{...}& 0& 0& 0& \mathrm{...}& 0& {\text{ξ}}_{\text{в}}^r& 1\end{array}\right|$                        (27)

При заданных граничных условиях однозначности III-го рода порядок матрицы (27) равен r_j – 1, причем отличны от нуля только коэффициенты прогонки, располагающиеся на трех диагоналях – главной и двух соседних. Трехдиагональная форма матрицы вида (27) позволяет организовать вычисления по методу Гаусса [17; 23; 24] так, чтобы не проводить операции с нулевыми элементами. Тем самым объем вычислений удастся значительно уменьшить. Достоинством схемы (26), (27) является то, что она позволяет реализовать прямоугольную сетку, наиболее апробированную при решении задач теплопроводности, а также; существенно упростить процессы получения объективных результатов. Кроме того, система уравнений (26) позволяет применить неявную форму представления конечно-разностных аналогов дифференциального уравнения (11), описывающего теплопередачу в стволе, что обеспечивает абсолютную устойчивость разностной схемы. Время решения задачи при этом существенно уменьшается путем выбора сколь угодно большой величины Δτ шага сетки по времени τ без риска нарушить устойчивость разностной схемы.

Однако выполнение требования устойчивости разностной схемы совершенно не означает выполнение условия аппроксимации исходной дифференциальной задачи. Поэтому существенно важным является решение вопроса выбора степени подробности сетки, влияющей на точность результатов. При этом необходимо учитывать, что чрезмерное уменьшение шагов Δx и Δr сетки по соответствующим осям x и r приводит к возрастанию числа неизвестных в разностной схеме (26), (27), а, следовательно, и к увеличению времени машинного счета. Укрупнение шагов Δx и Δr сетки по соответствующим осям x и r не способствует повышению точности требуемого результата. Сопряжение величины Δτ шага сетки по времени τ со значением аналогичной величины, задаваемой в процедурах решения основной задачи внутренней баллистики и периода последействия [1; 19; 20] определяет переменность отрезков времени τ. При отладке модели выбор значений шагов сетки Δx, Δr по соответствующим координатам x, r и величины шага по времени Δτ осуществлен, исходя из условий наибольшей точности решения, наименьшего числа неизвестных в системе разностных уравнений, учета реальных размеров области дискретного изменения аргументов W_Tkj величины T_kj температуры в kj-точках сечений ствола. Первоочередным в решении задачи представлялся учет конфигурации нарезов, поскольку их наличие приводит к неравномерности распределения температуры по периметру нарезной части канала ствола [8; 9; 11]. Исходное требование несравненной малости величина Δr шага сетки по оси r в отношении к высоте поля нареза очевидна.

Математический анализ степени подробности разностных схем представляет собой достаточно сложную задачу. Результаты исследований влияния значений шагов сетки Δx, Δr по соответствующим координатам x, r и величины шага по времени Δτ на точность решения задачи в данной работе приведены без доказательств. Наличие аппроксимации легко устанавливалось на практике путем численных экспериментов в процессе проведения пробных расчетов элементарных задач теплопроводности в цилиндрических стенках и сравнения полученных результатов с известными [25].

При измельчении пространственно-временнóй сетки точность решения задач возрастает, что свидетельствует о стремлении к нулю погрешности аппроксимации, но при этом, естественно, увеличивается время машинного счета. Приемлемыми следует считать шаги сетки по координате x – Δx = 0,001 м, по координате r – Δr = 4 · 10^–4 м, поскольку осредненная среднеквадратическая погрешность по всем точкам в этом случае не превышает 10 %, а время моделирования однократного нагружения – не более 17 с. Дальнейшее увеличение степени подробности сетки, например, в 2 раза, приводит к увеличению времени решения задач на порядок при увеличении точности счета всего на 3 %.

Так как скорость газов в канале ствола v во времени τ и по длине ствола l постепенно возрастает, то такая особенность не позволяет построить равномерную сетку по времени Δτ, поскольку по длине ствола l величина шага Δx сетки по оси x так же будет увеличиваться. Это, в свою очередь, может привести к тому, что точность результатов решения полученных в различных точках области дискретного изменения аргументов W_Tkj температуры в kj-точках сечений ствола T_kj, будет существенно отличаться между собой, что недопустимо. Учитывая также и то, что вычисления на каждом i-м временнóм слое производятся на основании значений предыдущего (i – 1)-го временнóго слоя, ошибка будет накапливаться достаточно быстро. В целях исключения этого события целесообразно использовать переменный шаг по времени, определяемый исходя из скорости газов в канале ствола v, получаемой решением основной задачи внутренней баллистики [1; 19; 20]

$\Delta \text{τ}=\frac{\Delta x}{v\left(l\right)}.$                                                                                                                   (28)

В отличие от пространственной сетки множество узлов τ_i, определяемых точками i = 0, 1, 2, …, I_ο, представляет собой неравномерную временнýю сетку в области W_ti

Формула (28) предъявляет жесткую связь величины Δτ шага сетки по времени τ со значениями Δx шага сетки по оси x, поскольку точность решения задачи непосредственно зависит от правильности выбора последнего.

Сходимость всех типов разностных схем при наличии условий их устойчивости и аппроксимации доказана в трудах [21–24].

Расчет температурного поля ствола при применении ААО сводится к многократному (по количеству выстрелов в очереди, количеству очередей выстрелов и интервалам времени между очередями выстрелов) решению системы уравнений (26) при начальном распределении температуры ствола T, которое устанавливается к началу очередного выстрела и определяется при решении той же системы уравнений (26) для предыдущего выстрела. Завершенной формой разработки модели явилась программа, позволяющая рассчитывать температурное поле ствола авиационной пушки ГШ-301.

Вследствие принятия некоторых допущений модель обладает определенным уровнем абстракции и, по причине неизбежной потери информации, не дает полную картину, характеризующую исследуемые физические процессы. Обоснование частных формулировок и дальнейшие обсуждения следствий представляются возможными после проверки адекватности модели реальным процессам теплообмена в стволе.

Проверка адекватности и ресурсоемкости модели

Достоверность разработанной модели установлена проверкой пространственного распределения нестационарных значений температуры ствола T при различных режимах стрельбы. Значения рассматриваемого основного теплофизического параметра для орудийной стали ОХН3МФА хорошо известны и позволяют сравнить результаты моделирования с экспериментальными данными, классифицированными, например, в [8; 10; 11] и результатами моделирования, наиболее близкими к опытным, зарегистрированными в работах [8; 9; 12; 13].

Анализ выходных данных показал, что расхождение результатов вычислений и эксперимента не превышает 10 % и отличие результатов моделирования от других теоретических расчетов не существенно. Сопоставление совместных расчетных и экспериментальных данных показывает незначительное (не более 1,3 %) увеличение точности получения совокупных нестационарных значений температуры ствола T. Вместе с тем оценка ресурсов, необходимых для проведения заявленных вычислений, обнаружила некоторые преимущества модели. В качестве оцениваемого ресурса рассмотрена временнáя эффективность [26], зависящая от списка, типа и структуры взаимодействия операторов программы и, естественно, от быстродействия вычислительной машины (тактовой частоты процессоров, объема оперативной памяти). Учитывалось то обстоятельство, что одна и та же программа при одних и тех же случаях на разных вычислителях выполняется за разное время.

С целью сравнения временнóй эффективности модели версии программ, ранее применяемые в работах [9; 12; 13], импортировались на один и тот же современный вычислитель. По итогам прямого хронометража и абсолютных количественных сравнений наблюдалось сокращение времени машинного счета с использованием модели в 1,6 раза, что считается вполне приемлемым для задач такого класса. Указанные выводы обусловливают возможность моделирования рабочих процессов выстрела и нагружения ствола с точностью, как минимум не снижающей прежнюю, но с немаловажным выигрышем во времени.

Экономии временных затрат на моделирование теплового состояния ствола способствовало, по-видимому, проведение следующих мероприятий:

Рекомендации по распространению решений при исследовании разнотипных образцов артиллерийских орудий

В качестве опорного образца в работе выбрано ААО типа ГШ-301, имеющее базовую одноствольную схему автоматики, стоящее на вооружении большей части современных, а также планирующееся на оснащение перспективных комплексов авиационного вооружения. Однако направленность исследований не исключает факта трансформирования разработанной модели к исследованиям других моделей ААО, прочих типов сухопутных и морских скорострельных артиллерийских орудий, а также артиллерийских орудий крупных калибров при выполнении ряда условий:

1. Решение основной задачи внутренней баллистики применительно к артиллерийскому орудию конкретной схемы действия и конструктивного оформления. Так, известны решения для безоткатных артиллерийских орудий, гаубиц, минометов и т. д. [19; 20].
2. Построение схемы осевой симметрии и сеточной схемы для области с конкретной геометрией. Отдание предпочтения тому или иному варианту в настоящее время является вопросом разрешаемым, но не закрытым [21–25]. Во многом предпочтения диктуются особенностями решаемой задачи, во многом зависят от вкуса исследователя.
3. Многофакторная композиция дифференциальных и разностных операторов, представляющая собой две стороны единого процесса моделирования. К примеру, несравненно бóльшая длина и относительно малая толщина стенок в аналогичных агрегатах танковых и противотанковых пушек заставляют дополнительно учитывать тепло от солнечной радиации. В части касаемой аппроксимации системы уравнений (11), отметим, что обычно после проведения сопоставления явной и неявной разностных схем, делается вывод о нецелесообразности применения первой. Однако практика решения реальных задач не подтверждает безусловную правильность такой рекомендации. В пользу явной разностной схемы можно привести следующие соображения. Во-первых, при анализе быстропротекающих процессов преимущество неявной разностной схемы, заключающееся в более свободном выборе величины временнóго шага, может не проявиться. Во-вторых, явные разностные схемы удобны при реализации на вычислителях с несколькими параллельными процессорами, которые получили широкое распространение в настоящее время.
4. Нахождение решения сложнейшей задачи о теплообмене между стволом и окружающими газами. Полная система дифференциальных уравнений теплоотдачи включает уравнения переноса теплоты, теплоотдачи, движения, неразрывности [16], а решение должно быть относительно простым и приемлемым для инженерной практики.
5. Непременный учет взаимного влияния нагретых стволов в ААО, высокотемпных сухопутных и морских артиллерийских орудиях многоствольной схемы служит критерием полезности выводов и повышает практическую значимость модели.

Заключение

Способами и приемами дифференциально-разностного описания быстропротекающих тепловых процессов образован один из возможных вариантов модели, позволяющей достаточно объективно, но с меньшими затратами времени производить расчет нестационарного температурного поля разностенного составного цилиндра. Определены направления модификации модели при включении в инструменты исследований многообразной типологии артиллерийских орудий.

Таким образом, предлагаемая модель может быть использована в научно-исследователь-ских организациях Воздушно-космических сил России для определения оптимальных условий применения ААО в составе комплексов авиационного вооружения ударных носителей; в проектно-конструкторских организациях при проведении расчетов, связанных с определением теплового состояния ствола; в образовательной деятельности для развития научных основ учебных дисциплин.

作者简介

Aleksandr Podkopaev

Air Force Military educational and scientific center “Air Force academy named after professor N. E. Zhukovsky and Y. A. Gagarin”

编辑信件的主要联系方式.
Email: aleksanpodkopaev@mail.ru

Cand. Sc., associate professor, professor of the department operation of aircraft weapon systems (and sighting systems)

俄罗斯联邦, 54а, Starykh bol’shevikov St., Voronezh, 394064

Azizullo Babadzanov

Air Force Military educational and scientific center “Air Force academy named after professor N. E. Zhukovsky and Y. A. Gagarin”

Email: azizullobabadzanov@mail.ru

adjunct

俄罗斯联邦, 54а, Starykh bol’shevikov St., Voronezh, 394064

Ilya Podkopaev

Air Force Military educational and scientific center “Air Force academy named after professor N. E. Zhukovsky and Y. A. Gagarin”

Email: ilya.podkopaev.96@bk.ru

adjunct

俄罗斯联邦, 54а, Starykh bol’shevikov St., Voronezh, 394064

Vasily Dolzhikov

Air Force Military educational and scientific center “Air Force academy named after professor N. E. Zhukovsky and Y. A. Gagarin”

Email: Link707@mail.ru

Cand. Sc., associate professor, head of the department operation of aircraft weapon systems (and sighting systems)

俄罗斯联邦, 54а, Starykh bol’shevikov St., Voronezh, 394064

参考

补充文件