Effect of Polymer Drag Reducer on Rheological Properties of Rocket Kerosene Solutions

Abstract

Adding drag reduction agent (DRA) to rocket kerosene is an effective way to reduce the pipeline resistance of rocket kerosene transportation systems. However, so far, there have been few research reports on the effect of DRA on the rheological properties of rocket kerosene solution, especially from a microscopic perspective. In this study, coarse-grained molecular dynamics simulations were conducted to investigate the rheological properties of rocket kerosene solutions with DRAs of different chain lengths and concentrations. The results showed that the viscosity of DRA—kerosene solution is generally higher than that of pure kerosene at a low shear rate, while with an increase in shear rate, the viscosity of DRA—kerosene solution decreases rapidly and finally tends to become similar to that of pure kerosene. The shear viscosity of DRA—kerosene solution increases with an increase in chain length and concentration of polymers. Through observing the morphologic change of DRA molecules and analyzing the radius of gyration and the mean-squared end-to-end distance of polymers, it was confirmed that the rheological properties of DRA—kerosene solutions are strongly related to the degree of entanglement of polymer chains. The simulation results provide microscopic insights into the rheological behavior of DRA—kerosene solutions and clarify the intrinsic relation between the morphologic change of polymer molecules and the rheological properties of DRA—kerosene solutions.

1. Introduction

The development of manned spaceflight, manned lunar landings, Mars exploration, and other major space activities requires the support of advanced rocket engines. According to the development trend for liquid propellant rocket engines [1,2,3], large-thrust liquid oxygen/kerosene engines and heavy carrier rockets will be the major design goals for the future. However, due to the increase in propellant flow rate and the influence of extreme heat flow, it is inevitable that pipeline resistance will increase greatly and that pump load will be large. Turbo pump power and generator temperature will greatly increase, which will affect the reliability of engines, limiting engine performance improvement [4]. Therefore, it is necessary to reduce the frictional resistance of rocket kerosene in the pipeline. Adding polymer additives to turbulent flow is a well-known method of hydrodynamic drag reduction.

The drag reduction effect was first observed by Tomes [5] in 1948, when he dissolved a small amount of polymer (polymethyl methacrylate) in monochlorobenzene and found that the flow pressure in turbulent pipes was significantly reduced compared with pure solvents. Similar observations were subsequently made in other scholars’ studies. For example, Virk’s research has shown that in turbulent boundary layers, dissolving one part per million of a long-chain flexible polymer in a solution can reduce turbulent friction losses by up to 80% compared to using the solvent alone [6]. This attracted wide attention and led to extensive research on drag reduction.

At present, the application of drag reduction agents in crude oil and refined oil pipeline transportation is very mature, but the drag reduction of rocket kerosene transportation systems is still in the stage of experimental demonstration. Lee et al. [7] investigated the turbulent drag reduction effect of the oil-soluble polymer polyisobutene in kerosene and explored the influence of molecular weight, Reynolds number (Re), and temperature on the time-dependence of drag reduction and explained the drag reduction phenomenon theoretically. Du et al. [8] conducted an experimental study on the flow resistance characteristics of rocket kerosene in pipelines and investigated the influence of various factors, such as the type of drag reduction agent (DRA), pipeline diameter, fluid flow rate, and concentration of DRA, on flow resistance. The results showed that flow resistance in the pipeline could be greatly reduced after adding DRA into the rocket kerosene and that the drag reduction efficiency could reach up to 75.05%. Zhang et al. [4] subsequently explored the drag reduction effect of DRA on high-temperature rocket kerosene under supercritical pressure and found that the drag reduction rate of drag reducers on rocket kerosene is affected by fluid temperature and Re and that the drag reduction rate of DRA could reach up to 60%. Li et al. [9] studied the drag reduction effect and heat transfer performance of supercritical kerosene with the addition of polymer additives by conducting CFD numerical simulations. The results showed that the pressure drop and heat transfer coefficient of supercritical kerosene flow was reduced by 46.8% and 37.5%, respectively.

Although the drag reduction effect of DRA has been studied experimentally and numerically, the mechanism of reducing the flow resistance of rocket kerosene by DRA is still unclear. The main reason is that it is still difficult to understand the interaction between the DRA molecules and the flow field. Therefore, in this paper, we study the effect of DRA in rocket kerosene from two aspects: rheological properties and molecular structure. From the macroscopic perspective, the rheological properties of rocket kerosene play an important role in hydrodynamic performance and drag reduction, while from the microscopic perspective, the morphologic change of the molecular structure of DRA is closely related to the flow field and energy transfer. In recent decades, a molecular dynamics (MD) method has been widely used to study the rheological properties of dilute polymer solutions [10,11,12,13,14,15], so we use an MD method to conduct simulations of rocket kerosene in this study. By studying the change in rheological properties of rocket kerosene before and after adding DRA and its relationship with the microstructure of DRA molecules, we can better understand the drag reduction mechanism of rocket kerosene by adding DRA, evaluate the performance of DRA in rocket kerosene, and then guide the practical application in engineering.

The paper is organized as follows. The model construction and simulation methods are described in Section 2. The rheological properties of the system are analyzed in the first part of Section 3 and the morphologic change in the molecular structure of DRA and related explanations are given in the remaining parts of the section. Finally, a summary of the research results is given in Section 4.

2. Materials and Methods

2.1. System Model

Rocket kerosene is a mixture of hydrocarbons with a boiling point between 150–300 °C which includes hundreds to thousands of components, such as alkanes, naphthenes, and aromatic hydrocarbons [16]. Due to the complexity of rocket kerosene components, it is difficult to directly construct a MD model containing all the components of rocket kerosene and the effective method is to select several representative components from the components to build a simplified model [17]. In this study, a rocket kerosene with a known distribution of carbon mass fraction was analyzed and it was found that alkane, cycloalkane, and bicycloalkane accounted for 95.4% of the total mass fraction of rocket kerosene, so only these three components were considered in the simulation. The average carbon numbers of the three components were calculated to be 13, 13, and 12, respectively. Based on the average carbon number of these three substances and the most common rocket kerosene composition in the literature, n-tridecane [18], n-heptylcyclohexane [19], and naphthalene, decahydro-2,6-dimethyl [20,21], were selected as the representative components of rocket kerosene in this paper. Each component was added according to its mass percentage. The number of molecules and atoms of each component is shown in

Table 1. The number of molecules and atoms included in the rocket kerosene model.

Species	Components	Molecular Formula	Molecular Number	Atomic Number
Paraffins	n-tridecane	C13H28	1210	15,730
Monocycloalkanes	n-heptylcyclohexane	C13H26	2244	29,172
Dicycloalkanes	Naphthalene, decahydro-2,6-dimethyl	C12H22	2827	33,924
	sum	-	6281	78,826

.

In order to simplify the simulation model and save computational resources, the Transferable Potentials for Phase Equilibria-United Atom (TraPPE-UA) force field [22] was used. Compared with other force fields, TraPPE-UA can well reproduce liquid properties, especially liquid density [23,24,25]. The idea of the Trappe-UA force field is to unite each carbon and its bonded hydrogen into a single interaction site. According to this idea, the molecular structure of each component is shown in Figure 1. Experimental studies have shown that polyisobutene used in combustion chamber cooling channels can significantly reduce pressure loss [26], so polyisobutene was selected as the drag reduction agent (DRA). The polymer chains were composed of NS monomers. At the initial stage, Np polymer chains were randomly distributed in the simulation box. Two polymer chains with different lengths were considered: NS = 36 and NS = 60. For the short-chain polymer, i.e., NS = 36, two different mole fractions were considered: ϕ = 0.0287 and ϕ = 0.086. The corresponding polymer chain numbers were Np = 5 and Np = 15, respectively. For convenience, these two solutions are denoted as ns36np5 and ns36np15. For the polymer of NS = 60, the polymer concentration was set to be ϕ = 0.086 and the corresponding polymer chain number was Np = 9, denoted as ns60np9. The initial size of the simulation box was determined to be 140 Å × 140 Å × 140 Å. The studies of Packmol [27] and Moltemplate [28] were used to construct the kerosene and the DRA–kerosene solutions. A snapshot of the simulation system is shown in Figure 2.

2.2. TraPPE-UA Force Field

A TraPPE-UA force field was used to conduct the MD simulations and its potential function is expressed as follows [29]:

$$E_{total}=E_{nonbond}+E_{stretch}+E_{angle}+E_{torsion}$$

(1)

where E_total, E_nonbond, E_stretch, E_angle, and E_torsion represent total potential energy, non-bonded potential energy, bond-stretching potential energy, bond-angle bending potential energy, and torsional potential energy, respectively. Non-bonded interactions, which represent interactions between pseudoatoms separated by more than three bonds or belonging to different molecules, can be described by pairwise-additive Lennard-Jones 12-6 potentials:

$$E_{nonbond}={\displaystyle \sum _{i,j}4}{\epsilon }_{ij}\left({\left(\frac{{\sigma }_{ij}}{r_{ij}}\right)}^{12}-{\left(\frac{{\sigma }_{ij}}{r_{ij}}\right)}^6\right)$$

(2)

where r_ij is the distance between atoms i and j and ε_ij and σ_ij represent the potential well depth and equilibrium distance between atoms, respectively. The LJ interaction parameters are given in

Table 2. Lennard-Jones potential parameters for the TraPPE-UA force field.

United Atom Type	Molecule	σ (Å)	ε/kB (K)
CH3 [22]	ALL	3.75	98
CH2 [22]	n-alkane/DRA	3.95	46
CH [30]	DRA	4.68	10
C [30]	DRA	6.4	0.5
CH2(cyc) [31]	Monocylic and dicycloalkanes	3.91	52.5
CH(cyc) [29]	Monocylic and dicycloalkanes	4.68	12

. The standard Lorentz–Berthelot combination rules are adopted for the cross parameters of the Lennard-Jones interaction between different pseudoatoms.

$${\sigma }_{ij}=\left({\sigma }_{ii}+{\sigma }_{jj}\right)/2$$

(3)

$${\epsilon }_{ij}={\left({\epsilon }_{ii}{\epsilon }_{jj}\right)}^{1/2}$$

(4)

In the TraPPE-UA force field, pseud-atoms are connected by fixed bond lengths (all bond lengths are 1.54 Å). In order to preserve the flexibility of the molecule, bond stretching was considered. The bond connecting two atoms is described by the harmonic bond-stretching potential:

$$E_{stretch}=\frac{1}{2}K{\left(R-R_0\right)}^2$$

(5)

where K is the elastic constant of bond stretching, R is the distance between adjacent particles, and R₀ is the equilibrium bond length. The parameters of bond-stretching potential energy are shown in

Table 3. Bond stretching parameters for the TraPPE-UA force field.

Bond Type	Molecule	R0 (Å)	K (kcal/mol/Å2)
CHx-CHy [32]	n-alkane/DRA	1.54	899.52
CHx(cyc)-CHy(cyc) [29]	Monocylic and dicycloalkanes	1.54	536

.

Two adjacent bonds form a bond angle θ whose motion is governed by a harmonic potential:

$$E_{angle}=\frac{1}{2}{K}_{\theta }{\left(\theta -{\theta }_0\right)}^2$$

(6)

where θ, θ₀, and K_θ are the measured bending angle, the equilibrium bending angle, and the force constant, respectively. The values of θ₀ and K_θ are listed in

Table 4. Bond angle parameters for the TraPPE-UA force field.

Bend Type	Molecule	θ (°)	Kθ/kB (K·rad−2)
CHx-CH2-CHy [22]	n-alkane/DRA	114	62,500
CHx(cyc)-CH2(cyc)-CHy(cyc) [31]	Monocylic and dicycloalkanes	114	62,500
CHx(cyc)-CH(cyc)-CH3 (or CH2) [30]	Monocylic and dicycloalkanes	112	62,500
CHx-C-CHy [30]	DRA	109.47	62,500

.

Four interconnected particles in the molecule form a dihedral angle, and the dihedral angle constantly twists in the simulation process. The potential energy that describes this torsion is the torsional potential energy, and its potential energy function is expressed as follows:

$$E_{torsion}=C_0+C_1\left(1+\cos \left(\phi \right)\right)+C_2\left(1-\cos \left(2\phi \right)\right)+C_3\left(1+\cos \left(3\phi \right)\right)$$

(7)

where $\phi$ is the dihedral angle and $C_{i=0:3}$ are the Fourier coefficients. Values of these coefficients are given in

Table 5. Dihedral torsion parameters for the TraPPE-UA force field.

Torsion	Molecule	C0/kB (K)	C1/kB (K)	C2/kB (K)	C3/kB (K)
CHx-(CH2)-(CH2)-CHy [22]	n-alkane	0	355.03	−68.19	791.32
CH2(cyc)-CH2(cyc)-CH2(cyc)-CH(cyc) [29]	Monocylic and dicycloalkanes	0	355.03	−68.19	791.32
CHn(cyc)-CH2(cyc)-CH(cyc)-CHn(cyc) (n = 1 or 2) [29]		−251.06	428.73	−111.85	441.27
CHx-(CH2)-(C)-CHy [30]	DRA	0	0	0	461.29

. It should be noted that the parameters of bond stretching, bending angle, and dihedral angle torsion of branches in monocycloalkanes and dicycloalkanes are consistent with those of linear chains with the same group. For the sake of distinction, pseudoatoms on normal alkanes are represented by CH_x, and pseudoatoms on cycloalkanes are represented by CH_x(cyc).

2.3. Simulation Details

In this paper, the RNEMD method (also known as Müller-Plathe method) proposed by Müller-Plathe [33] was used to calculate shear viscosity. The core idea of this method is to obtain shear viscosity by applying momentum flux, which is difficult to measure, and calculating the shear rate or velocity profile, which are easy to measure. The calculation formulas involved are shown in Equations (8) and (9):

$$j_y\left(p_x\right)=-\eta \frac{\partial v_x}{\partial y}$$

(8)

$$j_y\left(p_x\right)=\frac{p_{\mathrm{total}}}{2t{L}_x{L}_z}$$

(9)

where $\frac{\partial v_x}{\partial y}$ is the gradient of the velocity in the x direction along the y direction (also known as the shear rate). The momentum flux j_y(P_x) is the x component of momentum transported in the y direction per unit area at a given time. The proportional coefficient between shear rate and transverse linear momentum flux is shear viscosity η. P_total is the accumulation of momentum exchange in the x direction, t is the total simulation time, L_x and L_z represent the length of the simulation box perpendicular to the direction of momentum flux. Due to the periodic boundary conditions of the system, the coefficient 2 in Equation (9) is required.

According to the idea of the method, the simulation box is divided into M slabs along the y direction and then a momentum flux is artificially imposed along the y direction on the system. The process is as follows: firstly, find out the atom in the first slab with the largest momentum component in the −x direction and the atom in the (M/2) + 1 slab with the largest momentum component in the +x direction; secondly, exchange the x components of the momentum between the two atoms with the same mass and perform a momentum exchange every N time steps. As the simulation proceeds, this will produce a velocity profile in the simulation box. The shear rate is obtained by calculating the slope of the velocity profile, then the shear viscosity at a given shear rate is calculated according to Equation (8).

In this study, the simulation box was divided into 20 slabs along the y direction. During the simulation, the shear rate was controlled by adjusting the momentum exchange interval N. The simulation process of RNEMD was carried out at 293 K and 0.1 MPa. The whole simulation process was divided into two stages. In the first stage, the energy minimization process was performed and then a relaxation process of 3 ns was carried out in the isothermal and isobaric (NPT) ensemble to obtain a reasonable initial simulation configuration. The Nosé–Hoover method [34,35] was used for controlling the temperature and pressure of the system. The second stage involved data collection, including velocity gradient and momentum flux. The momentum exchange intervals were set to be 10, 20, 50, 100 and 200. The simulation time of this stage for the pure kerosene and the DRA–kerosene solution was 6 ns and 32 ns, respectively. For all simulation processes, the time step was 1fs. The trajectories of molecules were recorded every 20 ps. The snapshots of trajectories were visualized by OVITO [36]. All simulations were performed in LAMMPS (Large-Scale Atomic/Molecular Massively Parallel Simulator) [37].

3. Results and Discussion

3.1. Shear Viscosity of Different Solutions

In order to obtain the rheological characteristics and shear viscosity of DRA–kerosene solutions, the flow field of different shear rates was induced to pure kerosene and DRA–kerosene solutions by choosing different momentum exchange intervals. Figure 3 shows the velocity profile generated in the simulation box under different momentum exchange intervals. It can be seen that as momentum exchange interval decreases, velocity gradient and shear rate imposed on the simulation system increase. According to the velocity profile in Figure 3, the range of the magnitude of shear rate was found to lie between 1.0 × 10⁹ and 3.5 × 10¹⁰ s⁻¹.

The shear viscosities at different shear rates for different systems are shown in Figure 4. As can be observed from the figure, with the increase in shear rate, the viscosities of both the pure kerosene and the DRA—kerosene solutions decrease significantly, resulting in the phenomenon of shear thinning. The viscosities of the DRA—kerosene solutions are higher than that of the pure kerosene at low shear rates. With the increase in shear rate, the viscosities of the DRA—kerosene solutions decrease rapidly and finally approach that of the pure kerosene. A possible reason is that after adding DRA to the pure kerosene, the polymer chains appear to be bent and curled at a low shear rate and are prone to entangle with each other, which increases the resistance of polymer chains in the solutions and presents an increase in fluid viscosity. As the shear rate increases, the polymer chains are stretched along the flow direction and the entangled polymer chains undergo chain disentanglement [38] and the relevant flow resistance decreases rapidly, resulting in a decrease in viscosity of the solutions with the increase in shear rate.

It can also be seen from Figure 4 that shear viscosities of the four systems at low shear rates are listed in descending order as follows: ns60np9 > ns36np15 > ns36np5 > pure kerosene, which indicates that the viscosity of the DRA–kerosene solution increases with increasing polymer concentration. A possible reason is that when the polymer concentration of the solution is low, the polymer molecules are scattered in the solution and the intermolecular forces of polymers are weak; with the increase in polymer concentration, the distance between polymer molecules decreases, which increases their probability of contacting each other, along with the degree of entanglement, so that the hydrodynamic volume of polymer molecules increases and the viscosity of the solution increases [39,40,41]. Of all the solutions, ns60np9, with a higher molecular weight of polymer, has the greatest viscosity at low shear rates. This can be understood from two aspects. Firstly, the larger molecular weight of a polymer indicates that more segments are contained in a polymer chain, and the polymer chain is more likely to slip between the chains when moving. Therefore, more energy needs to be consumed during the shearing process to overcome the synergy between the chains [42,43]. Secondly, the increase in molecular weight of the polymer also enhances the probability of entanglement and thus leads to greater flow resistance and eventually increases the shear viscosity [44].

3.2. Microstructure Evolution of Polymer Chains in DRA—Kerosene Solutions

The conformational morphology of polymer molecules in the shearing process for ns36np5, ns36np15, and ns60np9 solutions at a low shear rate (N = 200) is shown in Figure 5. Different degrees of entanglement can be found in different solutions. In the ns36np5 solution, the polymer molecules were sparsely dispersed due to the small number of polymer chains. However, under the effect of Brownian motion of the polymer molecules, a slight entanglement of polymer molecules occurred in the shearing process. In the ns36np15 and ns60np9 solutions, due to the increase in the concentration and molecular weight of polymer chains, their chance of contacting each other was greatly increased, so the polymer chains tended to tangle with each other and form a network-like structure.

In order to offer an intuitive understanding of the conformational change of polymer molecules under the shearing effect, the evolution process of a single molecule in the ns36np15 solution with different shear rates corresponding to N = 200 and N = 10 is shown in Figure 6 and Figure 7, respectively. It can be seen that under the weak shearing flow (N = 200), the polymer molecule curls up into roughly ellipsoid form in most cases, while under the strong shearing flow (N = 10), the form of the polymer molecule alternates between curling and stretching states. Therefore, it can be considered that due to the presence of interactions between polymer molecules, entanglement between molecules will occur under low shear action and the degree of entanglement is enhanced with the increase in concentration and the molecular weight of polymers. With an increase in shear rate, the polymer molecules tend to stretch into elongated structures and change their orientation to the shearing direction, which destroys the entanglement of polymers, reduces flow resistance, and accordingly results in a decrease in the shear viscosity of solutions. Moreover, polymers in different solutions exhibit different stretching degrees under different shear rates, which will be introduced quantitatively in Section 3.3.

3.3. Quantitative Description of the Morphologic Change of Polymer Chains

In order to quantitatively compare the stretching degree of polymers under different shear rates, the radius of gyration and the mean-squared end-to-end distance of polymer molecules were calculated as follows:

$$R_{\mathrm{g},i}^2=\frac{1}{M}{{\displaystyle \sum m_{ia}\left({\overrightarrow{r}}_{ia}-{\overrightarrow{r}}_{cm,i}\right)}}^2$$

(10)

$$〈R_{\mathrm{ete}}^2〉=\frac{1}{N_p}\frac{1}{N_t}{{\displaystyle \sum _{i=1}^{N_p}{\displaystyle \sum _{j=1}^{N_t}\left(r_{i,1}\left(j\right)-r_{i,n}\left(j\right)\right)}}}^2$$

(11)

where M is the molecular mass of a single polymer molecule and m_ia and ${\overrightarrow{r}}_{ia}$ represent the mass and position of the monomer a in the polymer molecule i, respectively. ${\overrightarrow{r}}_{cm,i}$ is defined as the center-of-mass position of the polymer molecule i. The radius of gyration R_g is the averaged calculation result for all polymer molecules in the system. N_p is the total number of polymer chains in solution. N_t is the sampling times in the total simulation time. r_i,1 and r_i,n represent the first atom and the last atom of the polymer molecule i, respectively.

Figure 8 and Figure 9 show the calculation results of radius of gyration 〈R_g〉 and mean-squared end-to-end distance 〈R_ete²〉 at different shear rates. It can be seen from the figure that with the increase in shear rate, the variation trend of 〈R_g〉 and 〈R_ete²〉 is similar, i.e., with the increase in shear rate, both of them increase sharply first and then reach a plateau, indicating that with the increase in shear rate, the curled and entangled polymer molecules in solutions are rapidly stretched. When the shear rate reaches a certain level, the entanglement between polymer chains is almost destroyed, and the extension degree of polymer chains reaches its limit. Due to a longer chain of polymers in the ns60np9 solution, the 〈R_g〉 and 〈R_ete²〉 of the ns60np9 solution are larger than those of the other two solutions.

The variation of the normalized root-mean-squared end-to-end distance 〈R_ete〉/L with shear rate for different solutions is plotted in Figure 10, where L represents the length of fully stretched polymers. It can be seen that for all three systems, with the increase in shear rate, 〈R_ete〉/L increases sharply first and finally reaches a plateau. The value of the plateau region for the ns36np5 solution is close to 0.5, which is consistent with the results obtained by Xu et al. [38]. However, the value of the plateau region for the ns36np15 and ns60np9 solutions is slightly smaller, which is due to the fact that when the concentration or the chain length of polymers increase, the repulsive force between different polymer chains is stronger than that between different segments of the same polymer chain, resulting in incomplete stretching of polymers.

4. Conclusions

In this paper, the coarse-grained molecular dynamics method was used to study the rheological behavior of DRA—kerosene solutions and the conformational change of polymer molecules under steady shear flow conditions. The results show that DRA molecules in the solutions are easy to entangle and form network-like structures at low shear rates and that the degree of entanglement increases with the increase in the concentration and molecular weight of the polymer. This entanglement phenomenon increases the resistance to the motion of kerosene molecules, resulting in an increase in solution viscosity at low shear rates. With an increase in shear rate, the polymer molecules tend to be oriented along the flow direction and their molecular structure changes from a roughly ellipsoid form to a form alternating between curling and stretching, which is manifested by a sharp increase in the mean-squared end-to-end distance and the radius of gyration of polymers. At high shear rates, the network structures between polymer molecules are destroyed and the flow resistance decreases, so the shear viscosities of the solutions decrease greatly compared to those at low shear rates. With further increase in shear rate, the orientation of polymer chains is basically close to the limit state and the entanglement between molecules is completely destroyed. Correspondingly, the radius of gyration and the mean-squared end-to-end distance of polymers finally reach a plateau region. At this stage, the viscosities of DRA–kerosene solutions are consistent with that of pure rocket kerosene. The simulation results provide insight into the relation between the morphologic change of DRA molecules and the rheological properties of rocket kerosene solutions and have guiding significance for the application of DRA in aerospace working fluids.

Future research should focus on the mechanism of drag reduction of rocket kerosene and the exploration of better molecular structures of DRA. Meanwhile, rocket kerosene also plays a role as a coolant in regenerative cooling engines, and the addition of DRA will affect the heat transfer capacity of rocket kerosene. Therefore, it is also important to study the relationship between drag reduction and heat transfer in DRA–kerosene solutions.