Mathematical model of thermophysical loading of a small-caliber artillery barrel with variant discretization of half-integer layers of the computational domain

Full Text

Introduction

An analysis of existing trends in the development of artillery convincingly shows that at present the main attention of specialists is not so much the creation of new models, but rather the optimization of the tactical and technical characteristics of serial types of AAO [1]. An important obstacle when searching for the reserve functionality of AAO is manifested in the phenomenon of heating the barrel, which is cyclically subjected to high thermomechanical loads created by firing modes. The barrel largely determines the combat properties of the AAO, since it is in the barrel that the ballistic characteristics are realized and the design of all elements of the “cartridge-barrel” system largely depends on its design. As a result, the scientific and technical task of formalizing the temperature field of the barrel seems to be a priority task of AAO research.

The physical meanings of the automatic firing process indicate the need for an indispensable description of the non-stationary heating and cooling of the barrel by solving the differential equation of thermal conductivity and the conditions of uniqueness with variable, continuous and discontinuous coefficients [2]. However, the exact solution of the thermophysics equation is limited for a certain range of problems. Such problems include the multidimensional, nonstationary, nonlinear problem of heat transfer in a cylindrical wall with a cross section varying along its length. Without dwelling on the diverse variations of approximation schemes for the differential heat equation and uniqueness conditions in various subject areas, we note the most successful approaches developed by domestic and foreign scientists. Thus, in the articles [3–5], experimental research schemes and methods for processing output data are proposed that provide increased accuracy in determining body temperature and expanded the measurement range; the article [6] presents a unique thermal model developed based on the apparatus of probability theory; in the article [7], the temperature fields of finned walls of various configurations were determined by numerical solutions of the multidimensional heat conduction problem; the article [8] proposed tools for mathematical modeling (hereinafter referred to as modeling) of the temperature field in gas turbine units, taking into account as much as possible the set of parameters in multifactor boundary conditions of the boundary layer; in the article [9], correlation regression dependencies of the optimal extrema of loading barrels of small arms and cannon artillery weapons were obtained. Examples of works on similar topics in the field of aviation artillery science include the articles [10–13].

Despite the fact that in the analyzed works almost all of the presentation of the material, naturally, is of a purely specific nature, some ideas of colleagues turned out to be useful in achieving the goal of this work.

Formation of a model scheme for studying the temperature field of the barrel

Obtaining the desired solution to the problem posed in a non-stationary formulation, with thermophysical coefficients depending on temperature, is carried out in a sequence that ensures step-by-step specification of dependent actions.

Since the barrel has the shape of a limited cylinder of finite length, with the structural absence of heat sources in the internal sections of the barrel, the basic equation of thermal conductivity is presented in a cylindrical coordinate system in the form [2; 14; 15]:

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

where T – barrel temperature; t – time; a – thermal diffusivity coefficient of barrel steel; z, r, θ, – radius vector, applicate and polar angle, respectively, of the cylindrical coordinate system.

Coefficient a in the equation (1) is significant for non-stationary thermal processes and characterizes the rate of change in body temperature:

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

where λ, с, ρ – coefficients of thermal conductivity, specific heat capacity and density, respectively, of barrel steel.

If the thermal conductivity coefficient of barrel steel λ characterizes the ability of a material to conduct heat, then the thermal conductivity coefficient of barrel steel a is a measure of the thermal inertia properties of the body under study. The rate of temperature change at any point in the barrel will be greater, the greater the value of the coefficient a, which is revealed by the test condition when operating with formula (2), formed as a table of dependences of the thermal conductivity coefficients λ and specific heat c of the barrel steel on the barrel temperature T [16].

The most complete mathematical models of heat exchange processes occurring in various products with various configurations take into account the presence of uneven space-time fields in the desired quantities: temperatures of solids, liquids, gases, heat flows, radiation intensities, etc. [6–9]. Such mathematical models are systems of partial differential equations, integral and integrodifferential equations. However, the solution to the problem under consideration is limited to the construction of a model based on specific assumptions, which is explained by the following reasons:

– direct implementation of complete thermal mathematical models is possible exclusively for elementary volumes under simple boundary conditions;

– the use of an absolute mathematical model of the functioning of a pulsed heat engine is complicated by the difference in the boundaries of the AAO elements and a large number of not always deterministic initial data;

– the issue of harmonizing the accuracy characteristics of physical and mathematical methods with the available characteristics of computer time, memory and bit grid involves the consistent use of more simplified, compared to the full, mathematical models that describe the thermophysical loading of the barrel with varying degrees of detail.

When solving the problem of the most complete and objective determination of the temperature field of a barrel heated by firing, the following assumptions are made that relate to the basic assumptions of the subject area of knowledge:

– the initial temperature of the barrel is approximately equal to the ambient temperature (T₀ = T₂) or corresponds to its distribution over the surface of the barrel; subsequent loading with shots is characterized by the presence of a very specific temperature field of the barrel before each shot;

– the material of the barrel steel ОХН2МФА is considered isotropic and homogeneous, that is, the coefficients of thermal conductivity λ and specific heat capacity c of the barrel steel do not depend on spatial coordinates;

– the contact of the cartridge case with the chamber wall is assumed to be ideal, due to the tight pressing of the cartridge case under the influence of the pressure of the powder gases (hereinafter referred to as gases) when fired;

– the cartridge is represented as a model temperature concentrator and is simulated by a concentrated heat capacity with constant thermophysical characteristics.

The first and second assumptions about the mechanism of heat transfer in the barrel allow us to assume that there are no temperature fluctuations T on the outer and inner surfaces of the barrel sections after the shot. Then the isothermal surfaces remain cylindrical, having a common axis with the pipe, and the barrel temperature T will change only in the radial and longitudinal directions, that is, ∂T/∂θ = 0 and ∂²T/∂θ² = 0 [2; 14]. Of the three coordinates written in equation (1) for the three-dimensional case, when considering the applied axisymmetric problem of determining the temperature field of the barrel, two coordinates z and r will remain. In addition, since the barrel is a body of rotation and is symmetrical about the longitudinal axis, after some transformations carried out for the convenience of data grouping, formula (1) is reduced to the equation for finding a two-dimensional temperature field of the barrel on the plane (0, z, r):

$\frac{1}{a}\frac{\partial T}{\partial t}=\frac{\partial }{\partial z}\left(\frac{\partial T}{\partial z}\right)+\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial T}{\partial r}\right).$   (3)

It should also be noted here that the accepting of the extreme two assumptions determines the need, discussed above, to take into account in equation (3) the dependence of the thermal conductivity coefficients λ and specific heat capacity c of barrel steel, included in formula (2), on the barrel temperature T when studying applied issues of safe placement of the next cartridge in the barrel heated by shooting.

The basic differential equation of thermophysics (3) establishes a connection between temporal and spatial changes in temperature at any point in the barrel at which the phenomenon of thermal conductivity occurs. A differential equation of the form (3) can have an infinite number of solutions. Isolating from this set a solution that reflects the conditions of thermal interaction in the barrel and specifies the problem posed was carried out by adding geometric, boundary and physical conditions of uniqueness to equation (3). The boundary conditions of uniqueness are further understood as a set of initial and boundary conditions.

When arguing the geometric conditions of unambiguity, the world's lightest 30-millimeter aircraft gun GSh-301 with a unique single-barrel automation circuit, which is in service with most modern aircraft and is planned to equip future aircraft weapons systems, was chosen. Since the barrel is a symmetrical body of rotation relative to the longitudinal axis, the introduction into consideration of a truncated region consisting of internal Γ₁, external Γ₂ and vertical boundaries Γ₃, Γ₄ located on one side of the longitudinal axis of the trunk is quite sufficient. Fig. 1 shows a diagram of the axial symmetry of the AAO type GSh-301 barrel in a cylindrical coordinate system (0, z, r,), specifying the diagram presented in article [17] by including the boundary designations Γ₁ – Γ₄, required for further clarifications. As before, the z axis coincides with the longitudinal axis of the barrel, and the temperature distribution in each calculated cross section of the barrel is symmetrical relative to the channel axis T = T(r).

 

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

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

 

In the process of applying AAO, the flight of an aircraft, as a rule, is carried out in a quasi-steady mode v₂ ≈ const and, based on the first assumption, the initial conditions of the problem are written in the form:

T(z, r, 0) = T₂ = const.   (4)

The boundary conditions for the simulated process must reflect the conditions of thermal interaction between the environment and the surface of the body. In general, boundary conditions can be specified in several ways. In the theory of heat transfer, boundary conditions of four types are distinguished [2; 14]. First type boundary conditions are specified in the form of temperature distribution on the surface of bodies. A mathematical description of heat transfer by first type boundary conditions is used for given temperature changes at the boundaries of bodies or very intense thermal conductivity on surfaces, when the temperatures of the surfaces are close to each other. The range of such practical problems is limited, and first type boundary conditions are used mainly in estimation calculations. Boundary conditions of the second kind are specified by the distribution of heat flux density on the surface of the body. The physical essence of the heat exchange conditions corresponding to second type boundary conditions reflects the heating and cooling of bodies through radiation, when heat exchange occurs mainly according to the Lambert-Beer law with uniform heating of the surface of the body. Third type boundary conditions are specified in the form of a dependence of the heat flux density due to thermal conductivity from the body on the temperatures of the body surface and the environment. The mathematical description of the processes of heating and cooling a body is carried out by Newton's law. Analytical expressions for boundary conditions of the third kind have found wide application in studies of heat transfer at the boundaries of materials and substances. Fourth type boundary conditions (conjugation conditions) are specified as conditions for the continuity of the temperature field and conservation of energy on the contact surfaces of multilayer structures.

In research practices heat transfer in solid bodies flown around by gas flows, setting third type boundary conditions at the boundary between the body and the flow has found wide application. Also taking into account the fact that the barrels are not thermally insulated, when solving the problem of determining the temperature field of the barrel of the GSh-301 aircraft gun, we will set the boundary conditions in the form of ambient temperatures and the laws of heat exchange between this environment and the surface of the barrel, depending on the design characteristics and conditions functioning.

At the inner Γ₁ and outer Γ₂ boundaries of the barrel, we will set the dependence of the thermal conductivity coefficient of the barrel steel λ on the gas temperature T₁ and air temperature T₂, respectively.

At the inner boundary Γ₁ of the barrel, convective heat exchange will take place between hot gases and the barrel channel:  

${\left.-{\text{λ}}_{\text{с}}\frac{\partial T}{\partial r}\right|}_{r_0}={\text{α}}_1\left(T_1-T\right),$   (5)

where r₀ – inner barrel radius; α₁ – heat transfer coefficient from gases to the bore.

Here and below, the dependence of the quantities under consideration on the current time t is obvious.

We note that to calculate the boundary conditions of heat transfer in the barrel channel, it is necessary to determine the intra-ballistic parameters of gases from the solution to the main problem of internal ballistics, set out in article [18].

At the outer boundary Γ₂ of the barrel, convective heat exchange occurs between the incoming air and the outer surface of the barrel:

$-{\left.{\text{λ}}_{\text{с}}\frac{\partial T}{\partial r}\right|}_{r_{y_j}}={\text{α}}_2\left(T-T_2\right),$   (6)

where $r_{y_j}$ – thickness (outer radii) of barrel elements; α₁ – heat transfer coefficient from the outer surface of the barrel to the air.

The development of a mathematical model of heat exchange inside and in the vicinity of the barrel during near-wall flows of coolants, which makes it possible to determine the heat transfer coefficients from gases to the barrel channel α₁ and from the outer surface of the barrel to air α₂, present in formulae (5) and (6), respectively, is the subject of the article [13].

In accordance with the second assumption, the vertical boundaries of the barrel Γ₃ and Γ₄ are considered adiabatic, that is, the heat flow through these boundaries can be neglected:

${\left.\frac{\partial T}{\partial z}\right|}_{z=0}=\mathrm{0,}{\left.\frac{\partial T}{\partial z}\right|}_{z=l}=0,$   (7)

where l – barrel length.

During bursts of shots, the channel and the outer surface of the barrel have quite high temperatures, so it is necessary to take into account the design features of the AAO reference sample. Modeling of the process of functioning of the standard cooling system of the GSh-301 aircraft gun is realized by introducing a local heat transfer coefficient.

In order to increase the accuracy of modeling the temperature field of the barrel, the influence of the cartridge case located in the chamber during the shot was taken into account. Based on the third assumption, it is possible to schematize heat transfer by describing the phenomenon of thermal conductivity. Since the thickness of the shell wall is relatively small, it is assumed that it will instantly warm up to the gas temperature when firing T₁. The boundary condition on the surface of the chamber at characteristic points of the barrel, where direct contact of the cartridge case with the wall occurs, is formulated as first type boundary condition [2; 14]:

$T(z=\mathrm{0...0,175};r=0)=T_1.$   (8)

The nonstationary temperature field of the barrel is definable with the known differential equation of the process (3) and given additional conditions (4) – (8), which completely determine the boundary value problem:

$\left.\begin{array}{l}\frac{1}{a}\frac{\partial T}{\partial t}=\frac{\partial }{\partial z}\left(\frac{\partial T}{\partial z}\right)+\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial T}{\partial r}\right);\\ {\left.-{\text{λ}}_{\text{с}}\frac{\partial T}{\partial r}\right|}_{r_0}={\text{α}}_1\left(T_1-T\right);\\ -{\left.{\text{λ}}_{\text{с}}\frac{\partial T}{\partial r}\right|}_{r_{y_j}}={\text{α}}_2\left(T-T_2\right);\\ {\left.\frac{\partial T}{\partial z}\right|}_{z=0}=0;{\left.\frac{\partial T}{\partial z}\right|}_{z=l}=0;\\ T(z=\mathrm{0...0,175};r=0)=T_1;\\ T(z,r\mathrm{,0})=T_2=\text{const.}\end{array}\right\}$.   (9)

Thus, with a number of simplifying assumptions, the problem of loading the barrel is formulated in a complete form. However, as noted in the papers [2; 14; 19–22], the objective lack of an exact analytical solution to direct, multidimensional, unsteady, nonlinear heat transfer problems in areas with a complex boundary configuration leads to the need to use numerical methods.

Synthesis of a finite-difference scheme for calculating the temperature field of the barrel

For most structures of complex shape, which also includes the shaft design, the system of eigenfunctions and the spectrum of eigenvalues of the corresponding homogeneous problem are not known and not tabulated [19]. Therefore, for such bodies, in this case, it is convenient to use the finite difference method as the most universal [19–22].

The area of continuous change of the argument is replaced by a discrete set of points, the intersections of which form nodes, that is, the construction of a difference grid (hereinafter referred to as the grid), as well as the reduction of the system of partial differential equations (9) to a finite-difference scheme, that is, the composition of a system of finite-difference algebraic (hereinafter referred to as difference) equations are performed by analogy with the techniques described in the publication [17]. Some of the author's duplication of information is mediated by the concentration of classical physical and mathematical meanings of the question of heating and cooling the barrel.

The area Ω_T of continuous change in the arguments of the desired value T is replaced by a certain finite set of points lying in this region. The grid points for forming the finite difference of the function of the integer argument T_kj along the z axis are designated by k, and similar points along the r axis are designated by j. In accordance with the specifics of the problem being solved, the region Ω_T is transformed into the area for calculating the temperature T_kj at kj -points of the barrel sections.  In accordance with the selected coordinate system (0, z r,) in the direction of the z axis, the barrel is divided into ϑ equal parts ϑ = l / Δz, and in the direction of the r axis into ν equal parts ν = r_y / Δr, where Δz, Δr are grid steps at the corresponding coordinates; r_y – maximum barrel thickness. To do this, ϑ – 1 rays are drawn in the direction perpendicular to the z axis and ν – 1 rays are directed in the direction perpendicular to the r axis, as shown in Fig. 2. As a result of this partition, we have a grid consisting of a set of internal (in Fig. 2 indicated by ) and boundary (in Fig. 2 indicated by ) nodes. Since, in the case under consideration, Δz = l / ϑ = const and Δr = r_y / ν = const, then the set of nodes z_k, defined by points with numbers k = 0, 1, 2, ..., K_ϑ and the set of nodes r_j, defined by points with numbers j = 0, 1, 2, …, J_ν, is a uniform spatial grid in the area ${\Omega }_{{Т}_{kj}}\text{.}$

 

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

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

 

Unlike the previous version [17], here we consider two possible approaches to setting geometric conditions for uniqueness when the boundary nodes of the grid do not coincide with the boundaries of the barrel. One of them is the introduction of additional nodes at points where the grid lines do not coincide with the elements of the trunk geometry. The second approach is that the geometry of the trunk is approximated by lines passing through the boundary nodes of the grid, and the geometric conditions of uniqueness are transferred to these lines. Due to the inexpediency of introducing additional nodes, which leads to a significant complication of the problem of constructing a difference scheme, the second approach turned out to be more preferable, since it does not introduce additional difficulties in writing difference equations. The approximation of the trunk geometry is realized by conditionally dividing it into a finite number of sections, each of which is characterized by length and thickness, which are reduced to spatial grid steps Δz and Δr along the z and r axes, respectively.

By analogy with the grid for the spatial domain ${\Omega }_{{Т}_{kj}}$, a temporary grid of the domain ${\Omega }_{{Т}^i}$ for calculating the value T ⁱ is introduced in the set of nodes τ_i,, defined by points i = 0, 1, 2, …, I_ο, where i and I_ο are the current and boundary, respectively, grid points for the formation of the finite difference barrel temperature T over time t. The time grid step t is designated Δτ.

The solution to the non-stationary problem of thermal conductivity in the barrel predetermines the unconditional intersection of one-dimensional spatial grids in each direction with a time grid in the following form:

${\Omega }_{{Т}_k{Т}^i}={\Omega }_{{Т}_k}\times {\Omega }_{{Т}^i}=\left[\begin{array}{l}(z_k,{\text{τ}}_i),z_{k+1}\text{}=z_k+\Delta z,{\text{\hspace{0.17em}τ}}_{i+1}={\text{τ}}_i+\Delta \text{τ};\\ k=\mathrm{0,}\mathrm{1,}\mathrm{2,...,}{K}_{\vartheta };i=\mathrm{0,}\mathrm{1,}\mathrm{2,...,}{I}_{\text{ο}};\\ z_0=\mathrm{0,}{z}_{K_{\vartheta }}\text{}=\mathrm{1,5}\text{\hspace{0.17em}м},{\text{\hspace{0.17em}τ}}_0=\mathrm{0,}{\text{\hspace{0.17em}τ}}_{I_{\text{ο}}}=t.\end{array}\right]$

${\Omega }_{{Т}_j{Т}^i}={\Omega }_{{Т}_j}\times {\Omega }_{{Т}^i}=\left[\begin{array}{l}(r_j,{\text{τ}}_i),r_{j+1}\text{}=r_j+\Delta r,{\text{\hspace{0.17em}τ}}_{i+1}={\text{τ}}_i+\Delta \text{τ};\\ j=\mathrm{0,}\mathrm{1,}\mathrm{2,...,}{J}_{\text{ν}};i=\mathrm{0,}\mathrm{1,}\mathrm{2,...,}{I}_{\text{ο}};\\ r_0=15\cdot {10}^{-3}\text{м},r_{J_{\text{ν}}}\text{}=42\cdot {10}^{-3}\text{м},{\text{\hspace{0.17em}τ}}_0=\mathrm{0,}{\text{\hspace{0.17em}τ}}_{I_{\text{ο}}}=t.\end{array}\right],$   (10)

Expression (10) forms a stencil of a space-time grid, the diagram of which along the longitudinal coordinate z is shown in Fig. 3.

 

Рис. 3. К выбору пространственно-временной сетки (на примере пространственной координаты z и времени t)

Fig. 3. On the choice of a space-time grid (on the example of the spatial coordinate z and time t)

 

To construct difference analogues of differential operators of the system of equations (9), the method of formally replacing derivatives with finite-difference relations was used. This method is the most justified and applicable in problems of this class and is based on the Taylor series expansion of fairly smooth functions, which, as a rule, allows one to preserve the local properties of differential equations [15]. In addition, the method of approximating derivatives by Taylor series has two main advantages:

– when the size of the unit cell tends to zero, the difference equation is reduced to a differential equation, that is, the compatibility of the equations is ensured, which is an important criterion for stability;

– Difference equations of any degree of accuracy can be obtained by adding or removing the required number of terms in the approximating series, and if mathematical verification is necessary, the accuracy of the approximation is estimated from the discarded terms of the series.

The most natural way to replace the derivative is based on defining the derivative (for example, with respect to the z coordinate) as a limit [15; 19]:

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

If we fix the step Δz in equality (11), we obtain an approximate formula for the first derivative expressed in terms of finite differences.

For the so-called right difference relation or “forward” difference:

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

Similarly, the left difference relation (the “backward” difference) is introduced, written in the form:

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

When solving heat conduction problems, it is necessary to approximate the second derivative. For the second derivative, a linear combination of relations (12) and (13) is considered:

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

Each transition to one step “forward” is conventionally designated by “+1”, and “backward” by “-1”. Then, for the k grid point of the formation of a finite difference in the value of T_kj along the z axis, the right difference relation (12) is transformed to the form:

$\frac{\partial T}{\partial z}=(T_{k+1}-T_k)\frac{1}{\Delta z}.$   (15)

The left difference relation is transformed similarly (13):

$\frac{\partial T}{\partial z}=(T_{k+1}-T_k)\frac{1}{\Delta z}.$   (16)

The difference analogue of the second derivative, corresponding to formula (14), is represented by the relation

$\frac{{\partial }^{2}T}{\partial z^2}=(T_{k+1}-2T_k+T_{k-1})\frac{1}{\Delta z^2}.$   (17)

The formulas (11) – (17) and their justifications are also valid when replacing the derivative with respect to coordinate r in the system of equations (9) by difference relations. In this case, in analogue equations, instead of the variable z, the variable r will be present, and the index k will be replaced by the index j. We will keep in mind the discovered analogies further, sometimes without resorting to direct detailing of the difference scheme for the spatial variable r.

When constructing relations that approximate the time derivative ∂T/∂t in the system of equations (9), it is permissible to use temperature values at kj-points of the barrel sections at different times: T_k,j,i,  T _k,j,i–1,  T _k,j,i–2, … . However, in the practice of solving most applied problems of thermal conductivity, in the vast majority of cases, exclusively two-layer (in time t) difference schemes are used, approximating the values of the desired temperatures at the current i-th and previous (i – 1) time point. Much less frequently, the temperature values at the (i – 2)^nd moment of time are taken into ^st account by obtaining three-layer difference schemes [19–22].

When obtaining variants of two-layer difference schemes, the time derivative is approximated by the “backwards” time difference:

$\frac{\partial T}{\partial t}=(T^i-T^{i-1})\frac{1}{\Delta \text{τ}}.$   (18)

Spatial differential operators in a two-layer difference scheme are also approximated based on the temperature values T_kj at kj -points of the barrel sections at the i-th and (i – 1) ^st moments of time. In this case, two limiting cases are possible.

In the first case, only the temperature values T_kj at kj -points of the barrel sections for the current i moment of time are involved in the approximation. Thus, for the spatial variable z, the one-dimensional space-time approximation of the first differential operator to the system of equations (9) will have the form:

$\frac{{\partial }^{2}T}{\partial z^2}=(T_{k+1}^i-2T_k^i+T_{k-1}^i)\frac{1}{\Delta z^2}.$   (19)

In the second case, during approximation, only the temperature values T_kj at kj -points of the barrel sections for the previous time point (i – 1) ^st are used:

$\frac{{\partial }^{2}T}{\partial z^2}=(T_{k+1}^{i-1}-2T_k^{i-1}+T_{k-1}^{i-1})\frac{1}{\Delta z^2}.$   (20)

In accordance with options (18) – (20), we present two types of difference equations that approximate the first equation of system (9) in a one-dimensional version:

$\frac{\text{1}}{a}\frac{T_k^i-T_k^{i-1}}{\Delta \text{τ}}=(T_{k+1}^i-2T_k^i+T_{k-1}^i)\frac{1}{\Delta z^2};$   (21)

$\frac{\text{1}}{a}\frac{T_k^i-T_k^{i-1}}{\Delta \text{τ}}=(T_{k+1}^{i-1}-2T_k^{i-1}+T_{k-1}^{i-1})\frac{1}{\Delta z^2}.$   (22)

A difference equation of the form (22) makes it possible to express the solution to the problem of thermal conductivity in the wellbore in explicit form on the i time layer through the known solutions on the previous (i – 1)1^st layer. Difference equation (22) forms an explicit difference scheme. Algorithms for the numerical calculation of the system of equations (9) using an explicit difference scheme are quite compact when programming, but they impose requirements on computer time.

The difference scheme specified by a difference equation of the form (21) is more complicated, since each difference equation of the form (21), in addition to the unknown solution for the k-th spatial point, includes two more sought-after solutions for the neighboring (k – 1)-th and (k + 1)-th spatial points. All the sought solutions turn out to be “tied” with each other into a common non-degenerate system of difference equations. Thus, in this case, at each i-th time layer, the solutions are determined not by explicit formulas of the form (22), but from the solution of the system (K_ϑ – 1) of difference equations, as a result of which the difference scheme specified by the difference equation of the form (21) is implicit. Effective algorithms for solving the system of equations (9) using an implicit difference scheme are much more complicated than numerical algorithms using an explicit difference scheme, but the time for solving the problem can be significantly reduced by a rational choice of steps Δz, Δr and Δτ.

The obvious difference in the behavior of the solutions obtained in cases of implementation of the template in Fig. 3 using explicit (22) and implicit (21) difference schemes, a proper physical and mathematical explanation can be given. The value of the time derivative with an explicit difference scheme (22) is calculated from the values of the desired function at the beginning of the time interval, therefore the increment ($T_k^i-T_k^{i-1}$) does not depend on the obtained values, and the absolute value of this increment is proportional to the step. As a result, at some critical step Δτ, new values $T_k^i$, пcan be obtained that contradict the physical meaning of the problem (for example, a drop in the barrel temperature T on the i-th time layer compared to the (i – 1)-th time layer with continued exposure to gas temperature T₁). In the implicit difference scheme (21), the increment  ($T_k^i-T_k^{i-1}$) depends on all values $T_k^i$ on the new time layer, that is, there is a kind of “feedback” that does not allow obtaining absurd increments of the grid function. However, the practice of solving real problems does not at all exclude the advisability of including an explicit difference scheme in the stencil shown in Fig. 3. Firstly, when describing the fast processes under study, the advantage of the implicit scheme, which consists in a more free choice of the value of the time step Δτ, may not appear. Secondly, explicit schemes are more resource-intensive, especially when calculating on computers with several parallel processors, which are widely used nowadays.

Due to the uniformity of the grid over all spatial coordinates, the fact of difference approximation of the differential operator for the variable r for each value of z at any local point, both along an isolated coordinate r and when solving a problem with a time variable t simultaneously, can be shown in a similar way.

One of the most important achievements of computational mathematics is the development of various difference schemes for solving multidimensional partial differential equations of thermophysics [19–22]. The desire to obtain a close approximation of the problem of temperature loading of the barrel to its possibly existing exact analytical solution was facilitated by selection and some techniques for transforming the longitudinal-transverse difference scheme of the Peaceman-Rackford two-dimensional sweep method. The main advantages of the preferred explicit-implicit difference scheme include: a combination of the strengths of explicit difference schemes (low computer time consumption at the time step Δτ) and implicit difference schemes (unconditional stability, that is, the ability to ensure the accuracy of the solution at any degree of mesh detail); possibility of application to multidimensional areas and co-occurring processes; adaptability to compiling efficient machine codes on high-speed computers with a sufficiently large amount of RAM.

The course of the two-dimensional physical process of heating and cooling the barrel at each time step Δτ in spatial steps Δz and Δr is delivered as a result of the sequential implementation of one-dimensional processes, each of which begins from the distribution of the temperature field of the barrel that arose after the end of the previous one-dimensional process. Based on this representation, called splitting [20; 22] modeling of one-dimensional processes is carried out implicitly, and the sequential action of processes is taken into account in an essentially explicit way. Given the given boundary conditions and the same initial temperature T₀ at all points in the region of a complex-shaped barrel, the optimal solution is achieved by reducing the multidimensional problem at each time step Δτ to a set of one-dimensional problems solved by the sweep method.

The specificity of the stability of the implicit approximation of locally one-dimensional problems with any division of the time step Δτ determined the method for increasing the accuracy of the formation of an array of barrel temperatures T. The essence of the method is to select a template on the time grid containing a half-integer layer:

${\text{τ}}_{i+1/\text{2}}=\frac{{\text{τ}}_{i+1}-{\text{τ}}_i}{2}=\mathrm{0,5}\cdot \Delta \text{τ\hspace{0.17em},\hspace{0.33em}}$   (23)

as shown in fig. 4.

 

Рис. 4. К выбору временнóго шаблона продольно-поперечной разностной схемы методом двумерной прогонки Писмена – Рэкфорда

Fig. 4. To the choice of the time template of the longitudinal-transverse difference scheme by the method of two directions of Peaceman-Rackford

 

Then, taking into account difference equations (21) and (22), difference relation (23), as well as the discussed spatial analogies, the finite-difference approximation of the first equation of system (9) according to the longitudinal-transverse difference scheme for the direction z for any value of r will be look like:

$\frac{1}{a}\frac{T_{kj}^{i+1/2}-T_{kj}^i}{\mathrm{0,5}\cdot \Delta \text{τ}}=\frac{T_{(k+1)j}^{i+1/2}-2T_{kj}^{i+1/2}+T_{(k-1)j}^{i+1/2}}{{\left(\Delta z\right)}^2}+\frac{r_j{T}_{k(j+1)}^i-(r_{j-1}+r_j)T_{kj}^i+r_{j-1}{T}_{k(j-1)}^i}{r_j{\left(\Delta r\right)}^2}.$   (24)

The boundary and initial conditions along the z coordinate for each fixed value of r are approximated as follows

– initial condition:

$\left.\begin{array}{l}i=0T_{kj}^0=T_2,\\ i>0T_{kj}^i=T_{kj}^{i+1}.\end{array}\right\};$   (25)

– boundary conditions:

$\left.\begin{array}{l}{T}_{1j}^{i+1/2}=T_{2j}^{i+1/2},\\ T_{K_{\vartheta }j}^{i+1/2}=T_{(K_{\vartheta }-1)j}^{i+1/2}.\end{array}\right\}.$   (26)

When synthesizing a modified two-layer difference scheme, the solution to the non-stationary heat conduction problem on a separate layer can be considered as the initial condition for subsequent layers. Consequently, we write the finite-difference approximation of the first equation of system (9) for the direction r for any value of z in the following form:

$\frac{1}{a}\frac{T_{kj}^{i+1/2}-T_{kj}^i}{\mathrm{0,5}\cdot \Delta \text{τ}}=\frac{r_j{T}_{k(j+1)}^{i+1}-(r_j+r_{j-1})T_{kj}^{i+1}+r_{j-1}^{}{T}_{k(j-1)}^{i+1}}{r_j{\left(\Delta r\right)}^2}+\frac{T_{(k+1)j}^{i+1/2}-2T_{kj}^{i+1/2}+T_{(k-1)j}^{i+1/2}}{{\left(\Delta z\right)}^2}.$   (27)

The boundary and initial conditions along the z coordinate for each fixed value of r are approximated as follows:

– initial condition:

$T_{kj}^{i+1}=T_{kj}^{i+1/2}+T_{kj}^{i+1/2}$– solution to the equation (24);   (28)

– boundary conditions:

$\left.\begin{array}{l}-{\text{λ}}_{с}\frac{T_{k1}^{i+1}-T_{k2}^{i+1}}{\Delta r}={\text{α}}_1\left(T_1{,}_k^{i+1}-T_{k1}^{i+1}\right),\\ -{\text{λ}}_{с}\frac{T_{k(J_{\text{ν}}-1)}^{i+1}-T_{k{J}_{\text{ν}}}^{i+1}}{\Delta r}={\text{α}}_2\left(T_{k{J}_{\text{ν}}}^{i+1}-T_2{,}_k^{i+1}\right).\end{array}\right\}.$   (29)

Eliminating possible discrepancies, we note that in difference relation $T_1{,}_k^{i+1}$ and $T_2{,}_k^{i+1}$ denote, respectively, the temperatures of gases and air at the k-th spatial grid point on the (i + 1)-th time layer.

From equations (24), (27) it is clear that in the constructed difference scheme the transition from the i-th to the (i + 1)-th time layer occurs in two stages with steps of 0.5 Δτ = 0.5 (τ_i+1 – τ_i). Along with the main values of the grid function $T_{kj}^i$ and $T_{kj}^{i+1}$, intermediate values $T_{kj}^{i+1/2},$ are introduced which are formally considered as the values of T_kj at (τ_i+1 – 2τ_i+1/2). Relation (24) contains three unknown quantities $T_{(k+1)j}^{i+1/2},T_{kj}^{i+1/2},T_{(k-1)j}^{i+1/2},$ values $T_{k(j+1)}^i,T_{kj}^i,T_{k(j-1)}^i$ can be determined on the initial layer by integrating systems of equations of internal and intermediate ballistics [18]. That is, by relation (24) the difference scheme is classified as implicit in the z coordinate and explicit in the r coordinate. For any value of r, the numerical solution can be found by sweeping in the z direction. The desired temperature values T_kj at kj -points of the barrel sections are related to each other “horizontally” and “vertically”. Moreover, the unknowns of any internal horizontal straight line “interact” on the time half-layer exclusively with the unknowns of two adjacent straight lines – the upper and lower ones. Next, using relation (27), which contains three unknown quantities  $T_{k(j+1)}^{i+1},T_{kj}^{i+1},T_{k(j-1)}^{i+1},$ (the values $T_{(k+1)j}^{i+1/2},T_{kj}^{i+1/2},T_{(k-1)j}^{i+1/2}$ are recorded by sweeping in the z direction at values of r), the difference scheme is translated into a form that is implicit in the r coordinate and explicit in the z coordinate. Therefore, the final distribution of temperature T_kj at kj-points of the barrel sections is found by sweeping in the direction r at any value of z, where the transition between time layers is also performed in half-steps in the longitudinal and transverse directions, respectively, along the rows and columns on the grid.

The problem of optimal selection of grid steps Δz, Δr, Δτ and thus the number of its nodes is not easy. On the one hand, the greater the accuracy required, the finer the step is desirable. On the other hand, too small a step significantly increases the requirements for the speed and memory capacity of computers. Obviously, there must be some meshes with an optimal number of nodes. We will optimize the grid based on the conditions for the best convergence of the results of the numerical calculation with the likely existing true analytical solution and borrowed experimental data. 

First of all, in order to most accurately determine the temperature field of the barrel, it is advisable to solve the problem taking into account the configuration of the rifling, since their presence leads to uneven temperature distribution along the perimeter of the rifled part of the barrel bore [9; 10]. The initial requirement of incomparably small size Δr of the grid pitch along the r axis in relation to the height of the rifling field is obvious. In general, the value Δr of the grid step along the r axis is assigned according to the approximate dependence of the stationary and linear components [23]: 

$\Delta r\approx \frac{\text{λ}\cdot \Delta T}{{\text{α}}_1\left({Т}_1-{Т}_0-\Delta T\right)},$

where ΔT – temperature gradient on the heat exchange surface (for AAO ΔT ≤ 323 K).

Since the velocity of the projectile (gases) when fired v₁ in time t and along the barrel length l gradually increases, reaching the value v_д at the muzzle of the barrel, this feature does not allow constructing a uniform grid in time Δτ, since along the barrel length l the grid step size Δz along the z axis will also increase. This, in turn, can lead to the fact that the accuracy of the solution results obtained at different points in the area of discrete changes in the arguments  ${\Omega }_{{Т}_{kj}}$ of the value T_kj, will differ significantly from each other, which is unacceptable. Taking into account also the fact that calculations at each i-th time layer are performed both on the basis of the value of the previous (i – 1)-th and the previous (i – 0.5)-th time layer, the error will accumulate quite quickly. In order to eliminate this event, when calculating the heating of the barrel during the time of movement of the projectile (gases) along the barrel bore t_д, it is advisable to use a variable time step Δτ ≠ const, assigned when solving the main problem of internal ballistics [18]. Further, during the aftereffect period t_п and in the time intervals between bursts of shots Δt, a constant time step Δτ = const is established, assigned, in turn, for the period of intermediate ballistics:  

$\Delta \text{τ}=\left\{\begin{array}{l}\mathrm{0,0002}\frac{l}{v_1\text{(}l)},\text{\hspace{0.33em}\hspace{0.33em}если\hspace{0.33em}}t\le t_{\text{д}};\\ \mathrm{0,0002}\frac{l}{v_{\text{д}}},\text{\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}если\hspace{0.33em}}t>t_{\text{д}}.\end{array}\right.$   (30)

In contrast to the spatial grid, the set of nodes τ_i, defined by points i = 0, 1, 2, ..., I_ο is a non-uniform temporary grid in the area ${\Omega }_{{Т}^i}.$

The justification for the value Δz of the grid step along the z axis, providing the desired accuracy, of the solution was made using the stability condition for the explicit components of the difference scheme (24)–(29) [20; 22], including, among other things, the constancy of the time step Δτ of the lower part of formula (30):

$\frac{{\left(\Delta z\right)}^2}{\Delta \text{τ}}\ge 2a,\text{\hspace{0.33em}\hspace{0.33em}при\hspace{0.33em}}\Delta \text{τ}=\text{const,\hspace{0.33em}}a=\text{const\hspace{0.17em}}.$   (31)

Formula (31) shows a strict connection between the value Δz of the grid step along the z axis and the values Δτ of the grid step over time t, since the accuracy of solving the problem directly depends on the correct choice of the latter. From the stability condition (31) follows a guide to action - the refinement of the spatial grid must be accompanied by the refinement of the time grid. For example, when the number of spatial nodes z_k increases by 4 times, it is necessary to increase the number of time steps t of the difference grid Δτ by 16 times. Previously, the need to comply with condition (31) led to the fact that when determining the step size Δτ in solving real non-stationary problems of thermophysics, it was not possible to proceed only from the nature of the physical process being studied. This in some cases led to unacceptable costs of machining time. In addition, with an unreasonably large number of time nodes τ_i, a rounding error was observed, which occurs during numerical calculations in calculating machines of early generations.

The stability property of the explicit part of the difference scheme (24) – (29) has also been established in practice, by ascertaining the absence of “divergent mode” of the numerical solution in the process of trial calculations.

When considering the approximation property of the formed difference scheme (24)–(29), a special concept of the so-called total approximation was introduced [20; 22] of locally one-dimensional difference schemes, which is as follows. Each of the intermediate difference equations (24) or (27) separately may not have the approximation property. However, the discrepancy arising at the first time half-step, as a rule, is compensated at the second time half-step with the correct combination of spatial steps Δz, Δr and time step Δτ, so that in general the approximation error is obtained, tending to zero at the given degree of detail of the space-time grid.

Such a way of discretizing the computational domain should be recognized, although labor-intensive, but also the most acceptable for solving an applied problem of thermophysics.

Thus, a discrete set of grid points is characterized by coordinates and parameters:

$\begin{array}{l}{z}_k=(k-1)\cdot \Delta z,\Delta z=\mathrm{0,001}\text{\hspace{0.17em}м},k=\overline{\mathrm{1,}151};\\ r_j=r_0+(j-1)\cdot \Delta r,\Delta r=\mathrm{0,125}\cdot {10}^{-3}\text{\hspace{0.17em}м},j=\overline{\mathrm{1,}161};\\ {\text{τ}}_i=(i-1)\cdot \Delta \text{τ},\\ {\text{τ}}_{i+1/\text{2}}=\mathrm{0,5}\cdot \Delta \text{τ},i=\mathrm{0,}t.\end{array}$

Thus, the resulting expressions (24)–(29) constitute a method for numerically solving the boundary value problem (9) for determining the grid temperatures of the barrel $T_{jk}^i$. Taking into account nonlinearity of the first kind in the numerical solution of the system of equations (9) is organized by an iterative process, in which the determination of the next approximation is carried out by including a linear solution, in which the coefficients of thermal conductivity λ and specific heat c of barrel steel are calculated from the values of barrel temperatures T found on previous iteration. The formation of a non-stationary temperature field of the GSh-301 aircraft gun barrel is generally feasible by software organization of a matrix of values $T_{jk}^i$ when studying the application modes of AAO.

In the papers [19–22] it was proven that in the presence of approximation and stability, convergence of all types of difference schemes will always take place. However, this fact does not exclude the scientific and methodological significance of the procedure for checking the improvement in the convergence of the results of modeling the thermophysical loading of the system in comparison with the known results.

Checking the reliability of the thermophysical model of barrel loading

Applied research into the quality of AAO is preceded by checking the developed model for the adequacy of reflecting the simulated thermophysical processes occurring in a gas-dynamic pulse machine. Establishing the set of properties of the model that determine its suitability for conducting diverse numerical experiments is possible by comparing the modeling results with experimental data, as well as with the known averaged results of some theoretical works that are closest to the experimental data. This approach can significantly increase the reliability of the conclusions.

Based on these positions, the verification of the degree of objective representation by the calculation results of the actual values of the main parameter was carried out by numerical modeling of the process of heating and cooling of the barrel surface worked out at the test site during and after shooting a combat set of 75 rounds. The content of firing modes and conditions for the use of AAO are limited to the type of information that does not require further specifications. The combined results of the full-scale experiment and calculation are shown in the form of graphs of the dependence of the barrel temperature T on time t in Fig. 5.

 

Рис. 5. Зависимость температуры ствола авиационной пушки ГШ-301в районе компенсатора от времени при отстреле боевого комплекта в 75 патронов

Fig. 5. Dependence of the barrel temperature of the GSh-301 aircraft gun in the area of the compensator on time when firing a combat set of 75 rounds of ammunition

 

Analysis of the results obtained shows that there is a fairly good correlation between experiments and calculated data. Satisfactory agreement between the modeling results and the experimental data is confirmed by the fact that the averaged relative error in determining the barrel temperature T in the reference section does not exceed 0.6 %. In most works in the field of aviation artillery science, including the works of the co-author of the article, the discrepancy between this value in numerical and full-scale experiments is about 10 % [9; 11].

Thus, an increase in the accuracy of simulating thermal loading of the barrel by 9.4 % was achieved:

– taking into account the nonlinearity of the thermophysical properties of the barrel steel materialλ (T), c(T);

– choosing the values (probably close to optimal) of the grid steps Δz, Δr in the corresponding coordinates z, r, as well as the size of the step Δτ in time t in the thermophysical model of barrel loading;

–  an effective combination of the advantages of explicit and implicit difference schemes in the constructed explicit-implicit difference scheme for finite-difference approximation of the heat transfer problem in a body with a complex geometric shape.

The formalization of thermophysical processes of heat propagation in a thermally loaded AAO element is logically completed by a package of application programs designed to calculate the thermal state of the barrel during firing and determine safe firing modes for a range of flight conditions of the carrier aircraft [24]. Algorithms for the numerical calculation of the system of equations (9) using the corresponding finite-difference scheme (24) – (29) were debugged using the Microsoft Developer Studio software product, the Fortran Power Station 4.0 environment and the FORTRAN 90 algorithmic language.

Prospects for further improvement of the model

The software organization for calculating the temperature field of the barrel when using AAO comes down to multiple (according to firing modes) solution of the system of equations (24)–(29) with the initial distribution of the barrel temperature T, which is established at the beginning of the next shot and is determined by solving the same system of equations (24)–(29) for the previous shot. The proposed tools make it possible to adequately simulate the temperature field of the barrel under various firing conditions and create a basis for the composition of the maximum effective firing modes.

At the same time, there are practical applications of medium-sized special mechanical engineering, for which some of the assumptions adopted in the paper have to be removed. Thus, when analyzing the heating of a barrel in the area of gas outlet openings of gas automatics or muzzle devices, it is necessary to take into account local heat flows in the elements connected to the barrel. Then we should move on to a much more complex three-dimensional formulation in coordinates (0, z, r, θ). The need to solve three-dimensional heat transfer problems is not excluded when analyzing the effectiveness of cooling fins or grooves, the thermal state of the rifling, and taking into account the technological variation in the thickness of the barrel. In addition, when studying the mechanism of barrel wear when analyzing the thermally stressed state of a thin surface layer of metal adjacent to the channel surface, it is inevitable to take into account the dependence of the thermophysical characteristics of the barrel steel on not only temperature, but also on spatial coordinates. In practical calculations, it is increasingly necessary to abandon the assumption of the constancy of the thermophysical characteristics of the cartridge. Calculations based on the so-called “instantaneous” values of the thermophysical characteristics of the elements of ammunition located in the barrel during breaks between automatic firing help to clarify the thermodynamic state of the “cartridge-barrel” system and more closely link it with the combat properties of the AAO.

To obtain more complete information about the accuracy characteristics of the model, it is advisable to additionally conduct a series of flight experiments that provide natural conditions for thermal loading of the barrel. Then the assessment of the averaged relative error in modeling the heating and cooling of the barrel will undoubtedly be more objective.

Conclusion

By matching the accuracy characteristics of physical and mathematical methods for solving heat transfer problems and related problems with the colossal characteristics of speed, memory and bit grid of modern computer machines, a model of increased accuracy was synthesized, which differs from the known ones by the variable selection of the pitch of the template-grid of the barrel of the GSh-301 aircraft gun. The applied significance of the model is demonstrated by the availability of methods for its adaptation to solving other problems of thermodynamics and mechanics of the strength of barrels.

About the authors

Ilya A. 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

Russian Federation, 54а, Starykh bol’shevikov St., Voronezh, 394064

Aleksandr V. Podkopaev

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

Author for correspondence.
Email: aleksanpodkopaev@mail.ru

Cand. Sc., Associate Professor, Professor of the Department Operation of Aircraft Weapon Systems (and Sighting Systems)

Russian Federation, 54а, Starykh bol’shevikov St., Voronezh, 394064

Vasily I. 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)

Russian Federation, 54а, Starykh bol’shevikov St., Voronezh, 394064

References

Supplementary files