Thermoelasticity Coupling Transient Response of Composite-Material Direction Pipe

Abstract

Based on the positive scheme method, the thermal load of the jet flow in an inner composite-material direction pipe is obtained, and the thermoelasticity coupling transient response is investigated. The positive scheme method with second-order accuracy is extended for solving the axisymmetric Euler equations, and the supersonic axisymmetric jet flow over a missile afterbody containing jet exhaust is simulated. The correctness of the development for the positive scheme method is verified. With the developed positive scheme method used to simulate the jet flow in the inner direction pipe, the thermal load is obtained. The thermoelasticity coupling finite element model of the composite-material direction pipe is established, and the stress response under dynamic pressure, unsteady temperature, and coupling state is obtained. Results show that, at the beginning of engine ignition, the effect of dynamic pressure and temperature field on the coupling stress is basically the same, and after that, the contribution of the temperature field to the coupling stress increases, and the thermal stress is the main factor affecting the strength of the composite-material direction pipe.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

The composite-material direction pipe is an important component of a multiple launch rocket system. During the launch, the transient response of the direction pipe is one of the key factors affecting the initial disturbance and firing intensity of the launch system [1]. When the rocket is launched, the inner pipe is subjected to dynamic pressure, an unsteady temperature field, and boundary constraint. Under the thermal load of the jet flow, the transient response of the direction pipe shows a significant thermoelasticity coupling effect. Obtaining the thermal load of the inner jet flow of the direction pipe is the precondition for the study of its thermoelasticity coupling transient response. However, it is difficult to test the obtainment of the thermal load because the time of the formation and development of the rocket gas jet flow is very short. The high-resolution numerical scheme has become an important means to calculate the structure of wave systems, pressure fields, and temperature fields [2]. Demkowicz and Gopalakrishnan [3] proposed a Runge–Kutta discontinuous Galerkin method that had the same high-order accuracy and high resolution in the spatial and time dimensions. Chang [4–8] carried out a numerical calculation for the unsteady-state underexpanded jet flow with the CE/SE method and captured the complex wave structures such as shear layer, expansion waves, incident shock, and mach shock disk. Liu and Lax [9–12] proposed a positive scheme method for computational fluid mechanics which was used to solve hyperbolic conservation laws. At present, there are few reports regarding the positive scheme method in the calculation of gas jet flow. For the thermal-elasticity coupling problem of composite material, Yang et al. [13, 14] performed a study for the three-dimensional thermal-elasticity coupling problem of fiber-reinforced composite material under thermal-mechanical load with the finite element method. Palumbo et al. and Kylili et al. [15, 16] analyzed the thermal-mechanical stress for the failure mechanism of composite material and provided an analysis method for the thermal-mechanical stress of glass fiber-reinforced composite material. So far, the effects of the pressure field and the temperature field on the coupling state have seldom been studied. The steady-state response of composite materials has been the primary focus, and the study of the transient response at the millisecond level has been even less frequent.

In this study, the thermoelasticity coupling transient response of a composite-material direction pipe under the thermal load is researched. On the basis of the work of Fazio and Jannelli [10], the positive scheme method with second-order accuracy is extended for solving the axisymmetric Euler equations, the calculation code is programmed and applied to simulate the supersonic axisymmetric jet flow over a missile afterbody containing jet exhaust, and the correctness of the development for positive scheme method is verified. With the developed positive scheme method used in the numerical simulation for the gas jet in the direction pipe, the unsteady process of wave structure, pressure field, and temperature field is revealed. The thermoelasticity finite element model of the composite-material direction pipe is established, and the thermoelasticity coupling transient response of the direction pipe is calculated, which provides reference for the structural design.

2. Development of Positive Scheme Method

The numerical calculation of the oblique shock regular reflection problem indicates that the positive scheme method is of high calculation accuracy and perfect resolution [9]. To apply the positive scheme method to simulate the jet flow of a composite-material direction pipe, the method is developed to solve the axis-symmetric Euler equations.

2.1. Governing Equations

For supersonic jet flow, the governing equations can be expressed as inviscid, axisymmetric Euler equations in the cylindrical coordinate system [9] as$$\begin{array}{}\text{(1)}& \frac{\partial U}{\partial t}+\frac{\partial F\left(U\right)}{\partial x}+\frac{\partial G\left(U\right)}{\partial y}=S\left(U\right),\hspace{0.17em}\hspace{0.17em}\Omega \times \left(0,T\right),\end{array}$$where $U={\left[\rho ,\rho u,\rho v,E\right]}^T$, $F\left(U\right)={\left[\rho u,\rho u^2+p,\rho uv,\left(E+p\right)u\right]}^T$, $G\left(U\right)={\left[\rho v,\rho uv,\rho v^2+p,\left(E+p\right)v\right]}^T$, $S\left(U\right)=\left(-\rho v/y\right){\left[1,u,v,\left(E+p\right)/\rho \right]}^T$, $\Omega \in R^2$, $T$ is the time variable, $\rho$ is the density, $u$ and $v$ are the velocity components in the $x$ and $y$ directions, respectively, $E=\rho e+\rho \left(u^2+v^2\right)/2$ represents the total energy in unit volume, $e$ is the internal energy in unit mass, $p$ is the pressure for the ideal gas, the state equation is $p=\left(\gamma -1\right)\rho e$, and $\gamma$ is the specific heat capacity.

2.2. Space Discretization

The finite difference equations of the axisymmetric governing Equation (1) discretized using the positive schemes are given in the semidiscrete form as$$\begin{array}{}\text{(2)}& {\left(\frac{\partial U}{\partial t}\right)}_{i,j}=-\left[\frac{1}{\Delta x}\left(F_{i+1/2,j}-F_{i-1/2,j}\right)+\frac{1}{\Delta y}\left(G_{i,j+1/2}-G_{i,j-1/2}\right)\right]+S_{i,j},\end{array}$$where $\Delta x$ and $\Delta y$, respectively, are the spatial steps in the $x$ and $y$ axes. The construction process of positive scheme numerical flux, taking $F_{i+1/2,j}$, for example, is briefly described as follows. Other details of the calculation technique can be found in Reference [9].

The numerical flux $F_{i+1/2,j}$ in Equation (2) was constructed by mixing a centered difference flux $F_{i+1/2,j}^c$ and an upwind flux $F_{i+1/2,j}^{\text{up}}$. Introducing the limiter $L^0=R\mathrm{diag}\left({\varphi }^0\left({\theta }^k\right)\right)R^{-1}$, ${\varphi }^0\left(\theta \right)$ satisfies [17, 18]$$\begin{array}{}\text{(3)}& \left\{\begin{array}{c}0\le {\varphi }^0\left(\theta \right),\hspace{1em}\frac{{\varphi }^0\left(\theta \right)}{\theta }\le 2,\\ {\varphi }^0\left(1\right)=1.\end{array}\right.\end{array}$$

The numerical flux could be constructed as$$\begin{array}{}\text{(4)}& F_{i+1/2,j}^0=F_{i+1/2,j}^{0,\text{up}}+L^0\left(F_{i+1/2,j}^{0,c}-F_{i+1/2,j}^{0,\text{up}}\right),\end{array}$$where $F_{i+1/2,j}^{0,c}=\left(1/2\right)\left[F\left(U_{i,j}\right)+F\left(U_{i+1,j}\right)\right]$, $F_{i+1/2,j}^{0,\text{up}}=\left(1/2\right)\left[F\left(U_{i,j}\right)+F\left(U_{i+1,j}\right)\right]-\left(1/2\right)\left|\stackrel{\text{o}}{A}\right|\left(U_{i+1,j}-U_{i,j}\right)$, $\left|\stackrel{\text{o}}{A}\right|=R\left|\Lambda \right|R^{-1}$ is the absolute value of $A$, $A=\left(\partial F\left(U\right)/\partial U\right)$, and $\left|\Lambda \right|=\mathrm{diag}\left(\left|{\lambda }^k\right|\right)$. The numerical flux of the least dissipative scheme is given by$$\begin{array}{}\text{(5)}& F_{i+1/2,j}^0=\frac{1}{2}\left[F\left(U_{i,j}\right)+F\left(U_{i+1,j}\right)\right]-\frac{1}{2}\left|\stackrel{\text{o}}{A}\right|\left(U_{i+1,j}-U_{i,j}\right)+\frac{1}{2}{L}^0\left|\stackrel{\text{o}}{A}\right|\left(U_{i+1,j}-U_{i,j}\right).\end{array}$$

If $\left|\dot{A}\right|=R\mathrm{diag}\left({\mu }^k\right)R^{-1}$, where the diagonal matrix $\mathrm{diag}\left({\mu }^k\right)\mathrm{\ge }\left|\Lambda \right|$. Adopting the limiter $L^1=R\mathrm{diag}\left({\varphi }^1\left({\theta }^k\right)\right)R^{-1}$, ${\varphi }^1\left(\theta \right)$ is the minmod flux limiter that satisfies [17, 19]$$\begin{array}{}\text{(6)}& {\varphi }^1\left(\theta \right)=\left\{\begin{array}{c}0,\hspace{1em}\theta \le 0,\\ \theta ,\hspace{1em}0\mathrm{<}\theta \mathrm{\le }1,\\ 1,\hspace{1em}\theta \mathrm{>}1.\end{array}\right.\end{array}$$

The numerical flux can be constructed as$$\begin{array}{}\text{(7)}& F_{i+1/2,j}^1=F_{i+1/2,j}^{1,\text{up}}+L^1\left(F_{i+1/2,j}^{1,c}-F_{i+1/2,j}^{1,\text{up}}\right),\end{array}$$where $F_{i+1/2,j}^{1,c}=\left(1/2\right)\left[F\left(U_{i,j}\right)+F\left(U_{i+1,j}\right)\right]$ and $F_{i+1/2,j}^{1,\text{up}}=\left(1/2\right)\left[F\left(U_{i,j}\right)+F\left(U_{i+1,j}\right)\right]-\left(1/2\right)\left|\dot{A}\right|\left(U_{i+1,j}-U_{i,j}\right)$.

Similarly, the numerical flux of the more dissipative scheme is expressed as$$\begin{array}{}\text{(8)}& F_{i+1/2,j}^1=\frac{1}{2}\left[F\left(U_{i,j}\right)+F\left(U_{i+1,j}\right)\right]-\frac{1}{2}\left|\dot{A}\right|\left(U_{i+1,j}-U_{i,j}\right)+\frac{1}{2}{L}^1\left|\dot{A}\right|\left(U_{i+1,j}-U_{i,j}\right).\end{array}$$

Finally, by combining Equation (5) the least dissipative scheme and Equation (8) the more dissipative scheme, the numerical flux of positive schemes can be written as$$\begin{array}{}\text{(9)}& F_{i+1/2,j}^{\alpha ,\beta }=\frac{1}{2}\left[F\left(U_{i,j}\right)+F\left(U_{i+1,j}\right)\right]-\frac{1}{2}\left[\alpha \left|\stackrel{\text{o}}{A}\right|\left(I-L^0\right)+\beta \left|\dot{A}\right|\left(I-L^1\right)\right]\left(U_{i+1,j}-U_{i,j}\right).\end{array}$$under the condition $\left(\Delta t/\Delta x\right)\left(\alpha \underset{1\le k\le n,U}{\max }\left|{\lambda }^k\right|+\beta \underset{1\le k\le n,U}{\max }{\mu }^k\right)\le \left(1/2\right)$, where $\Delta t$ is the time step. The constants $\alpha$ and $\beta$ satisfy $0\mathrm{\le }\alpha \mathrm{\le }1$ and $\alpha +\beta \ge 1$.

2.3. Runge–Kutta Time Discretization

To match the accuracy with the space variables, we constructed the time discretization with the accurate second-order Runge–Kutta method [20] as$$\begin{array}{}\text{(10a)}& U_{i,j}^{\ast }=U_{i,j}^m-\left[\frac{\Delta t}{\Delta x}\left(F_{i+1/2,j}^{{\alpha }_x,{\beta }_x}-F_{i-1/2,j}^{{\alpha }_x,{\beta }_x}\right)+\frac{\Delta t}{\Delta y}\left(G_{i,j+1/2}^{{\alpha }_y,{\beta }_y}-G_{i,j-1/2}^{{\alpha }_y,{\beta }_y}\right)\right]+\Delta t{S}_{i,j},\text{(10b)}& U_{i,j}^{\ast \ast }=U_{i,j}^{\ast }-\left[\frac{\Delta t}{\Delta x}\left(F_{i+1/2,j}^{{\alpha }_x,{\beta }_x\ast }-F_{i-1/2,j}^{{\alpha }_x,{\beta }_x\ast }\right)+\frac{\Delta t}{\Delta y}\left(G_{i,j+1/2}^{{\alpha }_y,{\beta }_y\ast }-G_{i,j-1/2}^{{\alpha }_y,{\beta }_y\ast }\right)\right]+\Delta t{S}_{i,j}^{\ast },\text{(10c)}& U_{i,j}^{m+1}=\frac{1}{2}{U}_{i,j}^m+\frac{1}{2}{U}_{i,j}^{\ast \ast }.\end{array}$$

2.4. Verification of Supersonic Jet Flow

The positive scheme method was used to solve the axis-symmetric Euler equations, and the supersonic axis-symmetric flow over a missile afterbody was simulated. Half of the calculation domain cut by any meridian plane through the central axis of the jet is shown in Figure 1. The radial direction was 1.5, and the axial direction was 3.5. The calculation domain was dispersed by a uniform grid with $\Delta x=\Delta y=1/40$, the nozzle surfaces EF and FG used a mirror reflection boundary condition, the upper boundary DE used a supersonic flow condition, the outer boundary CD used a simple wave condition, the central axis of the flow AB used a symmetric condition, the lower boundary BC used an extrapolation boundary condition, and the other calculation domains used the initial condition of supersonic free flow.

[figure omitted; refer to PDF]

Figure 2(a) shows the density contours, Figure 2(b) the pressure contours, and Figure 2(c) the experimental picture under the same flow condition of a supersonic gas jet field.

[figures omitted; refer to PDF]

As shown in Figure 2, the exhaust jet was compressed by the supersonic freestream around the lip of the nozzle exit, also the Mach disc structure almost disappeared near the jet upstream. Meanwhile, two oblique shock waves, a contact discontinuity jet shock wave, and a bottle shock wave were formed. Behind the first oblique shock wave are the expansion waves around the lip of the nozzle exit, and under the second oblique shock wave is a jet shock wave. These expansion waves propagate across the jet flow and are reflected back into the jet as compression waves from the jet boundary, thus forming an incident shock. Later, a reflected shock wave appears when the incident shock encounters the jet centerline and intersects with jet shock wave. Overall, the several primary flow characteristics between the numerical results and the experimental picture [21] obtained under the same flow conditions as shown in Figure 2(c) agreed reasonably well, indicating that the numerical results obtained by the positive scheme method was reliable.

3. Positive Scheme Numerical Study for Gas Jet in a Direction Pipe

The nozzle outlet pressure of the rocket engine $p_{\text{c}}$ was 5.058 atm, the nozzle outlet density ${\rho }_{\text{c}}$ was 0.81 kg/m³, the Mach number of the nozzle exit ${\text{Ma}}_{\text{c}}$ was 2.63, the nozzle exit temperature $T_{\text{c}}$ was 1,780 K, the nozzle outlet diameter $d_{\text{c}}$ was 0.0613 m, and the inner diameter of the direction pipe $D$ was 0.07 m. Setting the characteristic length $L_0=d_{\text{c}}$, the characteristic density ${\rho }_0={\rho }_{\text{c}}$, and the characteristic speed $U_0=a_{\text{c}}=881\text{.}57$ m/s as three independent reference quantities, the dimensionless treatment of the above parameters was conducted.

Half of the calculation domain cut by any meridian plane that passes through the nozzle jet is shown in Figure 3, and the nondimensional calculation domain was [0 : 53.8] × [0 : 1.14]. The central axis AB of the direction pipe, the inner surface CD of the direction pipe, and the surface EF of the nozzle solid wall all used the mirror reflection boundary condition, the exit cone boundary AF of the nozzle used the cone inflow condition, the lower boundary BC used the extrapolation boundary condition, the upper boundary DE used the static atmosphere condition, and the static atmosphere condition was assigned to the initial jet flow.

[figure omitted; refer to PDF]

The density contour of the inner jet flow at different times (dimensionless time) is shown in Figure 4, which involves 1/4 the length of the pipe. As indicated from Figures 4(a)–4(f), after the unexpanded gas flowed through the nozzle into the atmospheric environment, the initial shock wave was formed by mutual compression together with the air inside the pipe. Owing to the intersection of the boundary of the jet flow with the direction pipe’s inner surface, when the jet flow collided with the direction pipe’s wall surface during its propagation process, a back jet flow together with two reflection shock waves were generated. As time developed, the two reflection shock waves intersected at the center axial of the direction pipe, while new reflection shock waves were penetrated and formed on the corresponding directional surface. Along with repeated intersecting and penetrating, the first intersection at the center axial was about 5.5 times the pipe’s diameter away from the outer pipe surface (shown in Figure 4(f)). The flow characteristics of the pipe were consistent with the calculation results of Reference [22], which demonstrates that the results of the positive scheme method were effective and reliable.

[figures omitted; refer to PDF]

The pressure and temperature curves of the inner surface of the direction pipe at t = 7 ms are shown in Figure 4. As shown in Figure 5(a), if the direction pipe is not long enough, the impinging shock intersects and reflects periodically inside the direction pipe, and the distance of the two pressure spikes caused by the impinging shock is approximately the same. At a certain length of the direction pipe, the structure of the impinging shock had nothing to do with the length of the direction pipe. The strength of the impinging shock decreased along the direction pipe, and the variation of pressure curves matched well with those of the relevant literature [23].

[figures omitted; refer to PDF]

4. Finite Element Model of the Composite-Material Direction Pipe

4.1. Calculation Model of the Direction Pipe

The structure of the composite-material direction pipe using a cylindrical and spiral design is shown in Figure 6. The front, middle, and back bourrelets were connected with the three clapboards of the launcher. The direction pipe plays a role in fixing the rockets in the state of marching, giving the rocket initial direction and a reasonable muzzle speed at launch time. There exists spiral guide groove in the inner surface of the direction pipe; by acting with the directional knob of the rocket, the rocket produces a rotational motion around its axis, which provides a suitable off-pipe rotational speed for the rocket.

[figure omitted; refer to PDF]

The finite element model of the composite-material direction pipe was established using commercial software ABAQUS [24]. In the PART module, the geometric model of the direction pipe was established according to the actual structural parameters. In the PROPERTY module, the composite shell section property was assigned to the geometric model; the composite material used the epoxy resin as the matrix and glass fiber as the reinforced material; the composite modes adopted the circumferential stacking, cross stacking, and axial stacking; the ply order was $\left[{\left({90}^{\circ }/\pm 45/0^{\circ }\right)}_2\right]$; the single layer thickness was 0.25 mm; the density was $\rho =1.8\times {10}^3{\text{kg/m}}^3$; the elasticity modulus was $E_1=1.3\times {10}^4\text{\hspace{0.17em}\hspace{0.17em}MPa}$, $E_2=1.3\times {10}^4\text{\hspace{0.17em}\hspace{0.17em}MPa}$, and $E_3=0.27\times {10}^4\text{\hspace{0.17em}\hspace{0.17em}MPa}$; Poisson’s ratio was ${\mu }_{12}=0.25$, ${\mu }_{23}=0.35$, and ${\mu }_{13}=0.35$; and the thermal expansion coefficient was ${\alpha }_1=1.1\times {10}^{-5}$, ${\alpha }_2=1.1\times {10}^{-5}$, and ${\alpha }_3=1.4\times {10}^{-5}$. In the STEP module, the thermoelasticity coupling algorithm was adopted. In the MESH module, in order to reflect the structure characteristic of the composite-material direction pipe, the 8-node reduced integration brick element was adopted; there were 1762 elements in all. In the LOAD module, the front bourrelet, middle bourrelet, and back bourrelet were constrained. The finite element model of the composite-material direction pipe is shown in Figure 7.

[figure omitted; refer to PDF]

4.2. Boundary Condition

As indicated in Figure 7, three bourrelets were connected with clapboards; therefore, the bourrelets of the direction pipe should meet the displacement boundary conditions$$\begin{array}{}\text{(11)}& \left\{\begin{array}{c}u\left(x,y,z,t\right){|}_{x=x_{0i},y=y_{0i},z=z_{0i},t=0}=0,\\ v\left(x,y,z,t\right){|}_{x=x_{0i},y=y_{0i},z=z_{0i},t=0}=0,\\ w\left(x,y,z,t\right){|}_{x=x_{0i},y=y_{0i},z=z_{0i},t=0}=0,\end{array}\right.\end{array}$$where $x_{0i}$, $y_{0i}$, and $z_{0i}\left(i=\mathrm{1,2,3}\right)$ represent, respectively, the three coordinate values of the elements of the bourrelets.

The thermal stress effect is generated when the composite material is forced by the temperature field of the gas jet; therefore, the finite element model of the composite-material direction pipe should meet the temperature boundary condition. The environment temperature of the outer surface of the direction pipe was $T\left(x,y,z,t\right){|}_{S_{\text{out}}}=300\hspace{0.17em}\hspace{0.17em}\mathrm{K}$. The gas jet temperature field of the inner surface of the direction pipe was $T\left(x,y,z,t\right){|}_{S_{\text{inner}}}=T\left(z\right)$, where the unsteady temperature field of the inner surface of the direction pipe is the function of the axis. coordinate $z$ along the direction pipe.

The experiment indicated that during the launching process, the pressure curve measured at a certain location was the same as that measured in each location along the inner surface of the direction pipe at a certain moment [23]. Based on this conclusion, by transforming the gas jet pressure curve and the temperature distribution curve from the spatial domain to the time domain, the pressure or temperature that acts on the elements or nodes of the direction pipe can be obtained. As pressure is a type of surface load, the pressure should be imposed on the elements of the direction pipe, while the temperature should be imposed on the nodes.

5. Thermoelasticity Coupling Analysis of Composite-Material Direction Pipe

5.1. Response Characteristics under Dynamic Pressure

The von Mises stress nephograms at different times of the inner surface of the composite-material direction pipe (the left side is the end of the pipe and the right side is the nozzle) under dynamic pressure are shown in Figures 8(a)–8(g). As indicated from Figure 8(a), in the initial period of the launching process, although the engine ignited, the gas jet did not impinge with the direction pipe. The maximum stress occurred near the bourrelets by the action of gravity. From Figures 8(b)–8(g), the maximum stress occurred in the locations where the gas jet impinged with the inner surface of the direction pipe to form an impinging shock. When the rocket moved in the direction pipe, the maximum stress area moved quickly from the end of the direction pipe to the nozzle, and the maximum stress area nearly departed from the direction pipe after 80 ms as shown in Figure 8(g).

[figures omitted; refer to PDF]

Four representative points in the inner surface of the direction pipe were selected along the axial direction, and the stress responses were compared. The distribution of characteristic points A, B, C, and D are shown in Figure 7. The circumferential stress and von Mises stress temporal change curves of the characteristic points are shown in Figures 9(a) and 9(b). As shown in Figure 9(a), the circumferential stress is the tensile stress, along with the movement of the rocket. When the impinging shock reached the point, a phase step occurred for the circumferential stress, and the curve assumed a pulse shape similar to that of the pressure curve. In view of Figure 9(b), a similar distribution law was obtained. As the locations of the characteristic points were different, the peak values of the curves differed. The peak value of the stress of the composite-material direction pipe under dynamic pressure was not large enough (the maximum was about 35 MPa), the stress of the characteristic points increased quickly when the impinging shock reached a certain location, and a pulse shape formed and then decreased quickly.

[figures omitted; refer to PDF]

5.2. Response Characteristics under Unsteady Temperature Field

Under the unsteady temperature field, the von Mises stress nephograms of the inner surface of the composite-material direction pipe at different times are shown in Figures 10(a)–10(g). From Figure 10, when the rocket passed through the direction pipe, a large stress area occurred because of the temperature field of the gas jet, and the stress area moved from the end of the direction pipe to the nozzle. The distribution laws of stress nephograms differed from those of the stress nephograms under dynamic pressure (Figure 8). In Figure 8, the stress area moved along with the rocket, while in Figure 10, the stress area moved ahead with the temperature field and was located mainly in the middle of the first and second bourrelets and in the middle of the second and third bourrelets. A high-temperature gradient existed between inner and outer surfaces of the direction pipe because the temperature field of more than 1,000°C acted on the inner surface of the direction pipe. The direction pipe expanded under the temperature load, but there were restraints from the clapboards in the bourrelets, and the areas near the bourrelets could not expand freely, thus resulting in a high thermal stress.

[figures omitted; refer to PDF]

The circumferential stress and von Mises stress temporal change curves of the characteristic points are shown in Figures 11(a) and 11(b). As shown in Figure 11(a), the circumferential stress of the characteristic points was tensile stress, and the peak values in characteristic points A and D were larger than those of characteristic points B and C. From the von Mises stress curve in Figure 11(b), when the combined influence of the stress components is taken into consideration, the peak values of characteristic points B and C were slightly larger than those of A and D. The stress distribution along the axial direction was inhomogeneous under unsteady temperature with the pulse-shaped stress response curve. The stress caused by the temperature field was larger than that of the dynamic pressure, indicating that the thermal stress had a greater influence on the strength of the direction pipe.

[figures omitted; refer to PDF]

5.3. Response Characteristics under Coupling State

The von Mises stress nephograms of the inner surface of the composite-material direction pipe under dynamic pressure and an unsteady temperature field at different time periods are shown in Figures 12(a)–12(g). By comparison of the stress nephograms of Figures 8, 10, and 12, the stress nephograms under the coupling state matched well with those under the temperature field but differed significantly from those under dynamic pressure, indicating that the thermal stress significantly affects the coupling stress and the temperature field significantly affects the coupling stress.

[figures omitted; refer to PDF]

The circumferential stress and von Mises stress comparison curves of the characteristic points A and B under dynamic pressure, the temperature field, and the coupling state are shown in Figures 13 and 14. As shown in Figure 13(a), the dynamic pressure stress, the thermal stress, and the coupling stress all increased and decreased quickly. After that, the stress curves did not fluctuate wildly, and the peak values of the stress caused by dynamic pressure and the temperature field were similar, as was the shape of the ascending and the descending segment of the coupling stress curve and the thermal stress curve. As shown in Figure 13(b), the consistent degree was clearer, which indicates that the thermal stress significantly affected the coupling stress. As shown in Figure 14(a), the dynamic stress, thermal stress, and coupling stress quickly increased to the peak value, and the gradient change was large. Slight fluctuations occurred and then tended to be steady after the curves decreased. The fluctuations were more regular after the first peak value than those of the corresponding curve of characteristic point A and the amplitude was larger. Similarly, the coupling stress curve and thermal stress curve matched well in the fluctuation segment where the first peak value increased and decreased. The circumferential stress under the temperature field and the coupling load was the compressive stress in the initial time when the first wave crest rose. In a very short time, the stress changed to tensile stress, and the circumferential stress under the dynamic pressure was always the tensile stress. As shown in Figure 14(b), the von Mises stress caused by the temperature field and the coupling load was almost the same except for the wave crest, and the leading role of thermal stress in the coupling stress is clearly shown.

[figures omitted; refer to PDF]

[figures omitted; refer to PDF]

6. Conclusion

(1)

The numerical results obtained by the developed positive scheme method matched well with the experimental results. The correctness of the developed positive scheme method was verified.

(2)

The jet flow of the direction pipe numerically calculated by the developed positive scheme method shows that the positive scheme method has a high ability for capturing complex gas jet structures. The flow characteristics and pressure curves of the pipe were consistent with the results of the reference, which further shows that the development of the positive scheme method was successful.

(3)

For the stress response of the pipe under dynamic pressure, the stress curve shows a steep peak as the impinging shock wave arrives, then the stress decreases quickly, and the stress of the composite-material direction pipe under dynamic pressure was small.

(4)

For the stress response of the pipe under the temperature field, the thermal stress along the axial direction was uniform, the stress curve presented a pulse shape; the peak value of the thermal stress was large.

(5)

For the stress response of the pipe under dynamic pressure and a temperature field, at the initial time when the engine ignites, the stress caused by the dynamic pressure and temperature field had the same influence on the coupling stress of the direction pipe, and then the contribution to the coupling stress from the temperature field increased.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

The authors would like to gratefully thank the financial supports from the Innovation Partnership Fund for Universities of China Aerospace Science and Technology Corporation (CASC), Natural Science Funds for Young Scholar of Jiangsu Province (No. BK20170837), National Natural Science Foundation of China (No. 51303081), Natural Science Foundation of Jiangsu Province, China (No. BK20130761), and Fundamental Research Funds for the Central Universities (No. 309181B8807).

References

[1] M. D. Watson, "System exergy: system integrating physics of launch vehicles and spacecraft," Journal of Spacecraft and Rockets, vol. 55 no. 2, pp. 451-461, DOI: 10.2514/1.a33929, 2018.

[2] W. Barsukow, P. V. F. Edelmann, C. Klingenberg, "A numerical scheme for the compressible low-mach number regime of ideal fluid dynamics," Journal of Scientific Computing, vol. 72 no. 2, pp. 623-646, DOI: 10.1007/s10915-017-0372-4, 2017.

[3] L. Demkowicz, J. Gopalakrishnan, "A class of discontinuous Petrov–Galerkin methods. II. Optimal test functions," Numerical Methods for Partial Differential Equations, vol. 27 no. 1, pp. 70-105, 2015.

[4] I. Chang, "Adaptive, unstructured, finite element, multimaterial, thermal analysis," Journal of Spacecraft and Rockets, vol. 30 no. 1, pp. 43-50, DOI: 10.2514/3.25469, 2015.

[5] D. Bianchi, A. Turchi, F. Nasuti, "Coupled CFD analysis of thermochemical erosion and unsteady heat conduction in solid rocket nozzles," Proceedings of 48th AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit, .

[6] R. Vinze, S. Chandel, M. D. Limaye, "Local heat transfer distribution between smooth flat surface and impinging incompressible air jet from a chevron nozzle," Experimental Thermal and Fluid Science, vol. 78, pp. 124-136, DOI: 10.1016/j.expthermflusci.2016.05.017, 2016.

[7] N. Yasaswy, V. Katti, S. Prabhu, "Experimental study on local heat transfer distribution between smooth semi-cylindrical concave surface and impinging air jet from a circular straight pipe nozzle," Proceedings of AIAA Thermophysics Conference, .

[8] B. Qiu, J. Yan, J. Liu, "Experimental investigation on pressure oscillation frequency for submerged sonic/supersonic steam jet," Annals of Nuclear Energy, vol. 75, pp. 388-394, DOI: 10.1016/j.anucene.2014.08.055, 2015.

[9] X. D. Liu, P. D. Lax, "Positive schemes for solving multi-dimensional hyperbolic systems of conservation laws II," Journal of Computation Physics, vol. 187 no. 2, pp. 428-440, DOI: 10.1016/s0021-9991(03)00100-1, 2003.

[10] R. Fazio, A. Jannelli, "Second order positive schemes by means of flux limiters for the advection equation," International Journal of Applied Mathematics, vol. 39 no. 1, pp. 25-35, 2009.

[11] S. Do, Y. Ha, M. Kang, "Application of a multi-dimensional limiting process to central-upwind schemes for solving hyperbolic systems of conservation laws," Journal of Scientific Computing, vol. 69 no. 1, pp. 274-291, DOI: 10.1007/s10915-016-0193-x, 2016.

[12] A. Kurganov, M. Prugger, W. Tong, "Second-order fully discrete central-upwind scheme for two-dimensional hyperbolic systems of conservation laws," Siam Journal on Scientific Computing, vol. 39 no. 3, pp. A947-A965, DOI: 10.1137/15m1038670, 2017.

[13] H. Yang, F. Li, T. W. Chan, "Effect of nanofiller shape on viscoelasticity of rubber nanocomposite investigated by FEA," Composites Part B: Engineering, vol. 92, pp. 160-166, DOI: 10.1016/j.compositesb.2015.10.020, 2016.

[14] A. Lal, N. M. Kulkarni, B. N. Singh, "Stochastic thermal post buckling response of elastically supported laminated piezoelectric composite plate using micromechanical approach," Curved and Layered Structures, vol. 2 no. 1,DOI: 10.1515/cls-2015-0019, 2015.

[15] D. Palumbo, R. De Finis, G. P. Demelio, U. Galietti, "Study of damage evolution in composite materials based on the thermoelastic phase analysis (TPA) method," Composite Part B, vol. 117, pp. 49-60, DOI: 10.1016/j.compositesb.2017.02.040, 2017.

[16] A. Kylili, P. A. Fokaides, P. Christou, "Infrared thermography (IRT) applications for building diagnostics: a review," Applied Energy, vol. 134, pp. 531-549, DOI: 10.1016/j.apenergy.2014.08.005, 2014.

[17] T. Xiong, J. M. Qiu, Z. Xu, "A parametrized maximum principle preserving flux limiter for finite difference RK-WENO schemes with applications in incompressible flows," Journal of Computational Physics, vol. 252 no. 11, pp. 310-331, DOI: 10.1016/j.jcp.2013.06.026, 2013.

[18] A. J. Christlieb, Y. Liu, Q. Tang, "Positivity-preserving finite difference WENO schemes with constrained transport for ideal magnetohydrodynamic equations," SIAM Journal on Scientific Computing, vol. 37 no. 4, pp. A1825-A1845, DOI: 10.1137/140971208, 2015.

[19] Y. Guo, T. Xiong, Y. Shi, "A maximum-principle-satisfying high-order finite volume compact WENO scheme for scalar conservation laws with applications in incompressible flows," Journal of Scientific Computing, vol. 65 no. 1, pp. 83-109, DOI: 10.1007/s10915-014-9954-6, 2014.

[20] S. C. W. Stanley, "Recent advances in computational fluid dynamics," US/ROC (Taiwan) Joint Workshop on Recent Advances in Computational Fluid Dynamics, vol. 43, 2013.

[21] J. Agrell, R. A. White, An Experimental Investigation of Supersonic Axisymmetric Flow Over Boattails Containing a Centered Propulsive Jet, 1974.

[22] F. J. Mao, Numerical Simulation of Gas Dynamic Effect of Rocket Launching System and its Engineering Application, 1995.

[23] Q. Xu, The Experimental Research and Testing Technology for the Gas Jet Dynamic Effect of the Aero-Defense Rocket, 1994.

[24] C. J. Mitchell, L. M. Dangora, J. A. Sherwood, "Investigation into a robust finite element model for composite materials," Finite Elements in Analysis and Design, vol. 115,DOI: 10.1016/j.finel.2016.02.003, 2016.