Generative adversarial networks with physical evaluators for spray simulation of pintle injector

Due to their wider throttling range and greater combustion stability, pintle injectors are especially suitable for liquid rocket engines that require deep, fast, and safe throttling,^1–3 such as the descent propulsion system in the Apollo program⁴ and the reusable Merlin engine of SpaceX.⁵ 

In practical throttleable engine applications, the pintle is movable to alter the injection area so that the mass flow rate of the injected propellants can be varied continuously according to the economical and safe thrust curve in a given situation.⁶ However, in the previous spray simulations of pintle injectors, the changes were only considered under discrete condition combinations over a limited number of select operating points.^7–9 For the traditional discrete methods they used, simulations have to be conducted repeatedly to vary the operating conditions and the computational cost becomes prohibitively expensive.¹⁰ Innovations for the spray simulation of the pintle injector are needed to address this issue.

Contrarily, the machine learning approach, especially the Neural Network (NN), has demonstrated its efficiency to predict the flow fields under different conditions with a single surrogate model.^11,12 Previous research studies on flow field prediction using the NN are mainly focused on the data-driven method. In addition to the indirect way using the closure model,^13,14 the field solution can also be directly obtained from the network model, which is trained with a large number of samples.^15–18 However, some predictive results obtained by data-driven methods may still exhibit considerable errors against physics laws or operating conditions.^19–21 In addition, in some sparse data regimes, some machine learning techniques lack robustness and fail to provide guarantees of convergence.²² For the purpose of remedying the above-mentioned shortcomings of data-driven methods, the physics-driven/informed methods are proposed recently.^23,24 By providing physics information, NNs are able to directly obtain field solutions that obey physical laws and operating conditions.²⁵ In these works, Partial Differential Equations (PDEs) were employed in the loss function to explicitly constrain the network training.^26,27

In the state-of-the-art neural network methods, Generative Adversarial Networks (GANs) proposed by Goodfellow et al.²⁸ are efficient to generate the instantaneous flow fields.^29,30 Despite the impressive performance for unsupervised learning tasks, the quality of generated solutions by GANs is still limited for some realistic tasks.³¹ In addition, as shown in the training results later, the transient nature of the spray injection and liquid sheet break results in the extreme difficulty of usual networks to qualify the place and intensity of dominating characterizations.

In this paper, based on one specific type of GAN and the idea of physics constraint, a novel Generative Adversarial Networks with Physical Evaluators (GAN-PE) framework is proposed. By introducing mass conversation and spray angle models as the two evaluators, this framework has a better training convergence and predictive accuracy. The trained model is able to simulate the macroscopic morphology and characterization of the instantaneous flow fields under different conditions. This paper is organized as follows: We first introduce the experimental settings and dataset acquisition. Second, the architecture of GAN-PE and the detailed parts are described. Then, the learning results of numerical experiments are presented for validation. Finally, conclusions are drawn.

Our training data are extracted from the spray experimental results of the pintle injectors.

The non-reactive cold experiments were conducted at atmospheric pressure. The dry air is used for axial flows and filtered water for radial flows. The schematic of experimental facilities is shown in Fig. 1. A back-lighting photography technique is used for instantaneous spray image visualization. The image acquisition system consists of a light-emitting diode (LED) light source, a high-speed camera, and a computer. The exposure time is 10 µs, and the frame rate is 50k fps.

The detailed gas–liquid pintle injector is shown in Fig. 2. In order to study the influence of the momentum ratio on the spray angle, the experimental device is designed to use the replaceable parts. In the experiment, the height of the radial liquid jet outlet L_open and the thickness of the axial gas sheet T_gs are adjusted by changing the height of the sleeve and the axial gap distance, respectively. When the liquid propellant is injected radially from the two sides of the pintle end through the manifold, the liquid columns are formed. These columns are broken by the axial gas propellant injection from the gap cling to the pintle. Finally, due to impingement and collision, the liquid columns break and form a plane conical spray like a hollow-cone atomizer. This design induces vigorous mixing of the gaseous and liquid propellants, which yields a high combustion efficiency.³² 

The spray experiments are carried out with the throttling level L_t of 40%–80%. L_t is varied by the linear adjustment of the height of the radial liquid jet outlet and the thickness of the axial gas sheet. The radial liquid jet outlet heights at throttling levels of 40%, 60%, and 80% are 2, 3, and 4, respectively. When L_t is fixed, the height of the radial liquid jet outlet is fixed and equal to the thickness of the axial gas sheet. Table I shows the operating conditions of the experimental campaign and the corresponding key specifications of the pintle injector. The temperature of the liquid and gas is 298.15 K. The mass flow rate m_l is determined by the variation of liquid pressure upstream. Eventually, 35k raw images were captured and 29k of them were used for training and the others for validation.

As shown in Fig. 3, to measure the spray angle, the spray images obtained in the experiment are post-processed to clarify the spray boundary. The average of ten images with the same time interval is used to measure the spray angle manually. Every operation condition has 1k raw images, so the 100th, 200th, …, 900th, and 1000th images are averaged. Then, the spray angles of the time-averaged spray images, defined as $θ=12(θ1+θ2)$⁠, are manually measured. Since the raw images are all captured in the steady injection stage and no temporal fluctuations, the measured angle value is unique per operating condition. Note that the average images and the corresponding manually measured angles are only used to train the spray angle estimator, i.e., the down-sampling CNN in the spray angle evaluator. The raw images are used in the training of GAN-PE.

The resolutions of the instantaneous spray image are 640 × 480. In order to reduce the training cost, the images are interpolated to the images with a resolution of 128 × 128. While the measured angle values, which represented the nature of the spray phenomenon, are fixed despite the image scaling.

Here, a Generative Adversarial Networks with Physical Evaluators (GAN-PE) framework is proposed. As shown in Fig. 4, the GAN-PE is composed of four parts: generator (G), discriminator (D), and two physical evaluators. The field solutions are generated by G, and the other three parts are employed to guarantee that the outputs catch the spray morphology and obey the operating conditions. The GAN is the base of the proposed network framework; the G captures the real spray data distribution, which corresponds to the operation conditions; and the D estimates the probability that a condition–sample pair came from the training data rather than G. There are also two evaluators designed to improve the performance of GANs. The first evaluator, Mass Conservation Evaluator (E_MC), is used to improve the generation robustness by calculating the ring error between output and the corresponding average target. The second evaluator, spray angle evaluator (E_SA), is used to improve the predictive accuracy in the specific operating conditions by comparing the output angle with the theoretical one. Fed with the outputs from G, the losses of D, E_MC, and E_SA are calculated, respectively. After that, backpropagation is applied to adjust the U-net Convolutional Neural Network (CNN) of G to generate a new spray field that more satisfies the conditions and prior physics knowledge. After enough iterations, the network will be able to generate a “correct” spray field.

From inputs toward outputs, the network of G consists of two parts: encoding and decoding.³³ In the encoding process, the operating conditions L_open, T_gs, m_g, and m_l are resized as four feature channels of the input tensor for convolutional down-sampling with corresponding kernels.^34,35 After that, the matrices with a size of 128 × 128 are progressively reduced to 512 single-value vectors. Each layer of the network consists of a convolution operation, batch normalization, and a non-linear activation function. By the convolutional calculation, along with the increasing number of feature channels, the matrix size is down-sampled by a factor of 2. In this way, the information of operating conditions is translated into the extracted features in the next layer. In addition, skip-concatenations from input to output feature channels are introduced to ensure that operating condition information is available in the following up-sampling process for inferring the solution. Then, the decoding part works in an opposite way, which can be regarded as an inverse convolutional process mirroring the behavior of the encoding part. Along with the increase in spatial resolution, the spray fields are reconstructed based on the single-value vectors by up-sampling operations. For more details of the U-net architecture and convolutional block, including active function, pooling, and dropout, see Ref. 27.

The weighted loss function considering the following discriminator and evaluators is written as

where $LD$⁠, $LEMC$⁠, and $LESA$ are the loss terms that are calculated by D, E_MC, and E_SA, respectively. In addition, α and β are the constant hyperparameters that are manually tuned before training to adapt the scales of these loss terms. Here, the orders of magnitude of α and β are 2 and 1, respectively. After proper training, both generator and discriminator losses remain stable and the generator is able to map a spray sample from a random distribution to the desired one that obeys the physical knowledge and conditions.

The discriminator in a GAN is simply a classifier. It tries to distinguish real samples from the data created by the generator, i.e., fake data. The discriminator’s training data come from two sources. One consists of the real data instances, here are the real experimental images. D uses these instances as positive examples during training. The other are the fake data instances created by the generator. The discriminator uses these instances as negative examples during training. Then, D is used to feed the possibility that samples come from the targets rather than generation distribution back to G. We use Least Squares Generative Adversarial Networks (LSGANs) settings to train the D and the G simultaneously.³⁶ This special type of GAN helps to remedy the gradient vanishing by using the least squares loss function instead of the sigmoid cross entropy loss function.³⁷ 

Here, D is modified by the encoder of the G, which means that the generating solutions from G are down-sampled by the re-convolutional calculation so that the spray field information is concluded into the linear 1D tensor. Then, this 1D tensor will be used to be trained to maximize the probability of assigning the correct label, real or fake, to both training targets and generating solutions. Similar to the work in Ref. 38, we use the input–output and input–target pairs to feed D instead of only output and target in the random image generation tasks. The operating conditions and the outputs/targets are concatenated as the different feature channels in a 4D data tensor. In this way, D not only discriminates between the real and fake but also helps to judge whether the outputs accord with the corresponding conditions.

The loss functions for LSGANs are defined as

where x is the training data and z is the input variables. In addition, a and b are the labels for fake data and real data, respectively; and c denotes the value that G wants D to believe for the generated solutions. Here, we apply a = 0 and b = c = 1. Therefore, $LD$ is equal to the second part of Eq. (2).

As shown in Fig. 5, we assume that there are a few spherical volumes with different diameters that are tangent at the middle point of the upper boundary in both generating images and the average images. The idea originates from that the spray phenomenon obeys the mass conservation law. Because the propellants are injected from the two flanks of the central pintle and the width of the image’s normal direction is fixed, the 3D effect is ignored and the spheres are simplified to the 2D rings. The mass fluxes of droplets through one specific ring in every instantaneous frame are equivalent. The spray field generation is indeed a 2D image reconstruction task. Therefore, the “mass” here is a broader concept, which also involves the background shadow, i.e., the quasi-mass in one position is represented by the gray value of the corresponding pixel. Following the definition of “L1loss” that is widely used in the machine learning community, we define a mass conservation loss here. The difference is that the former measures the sum of absolute errors between each element in the generation and target,³⁹ but ours first calculates a gray value summation of every element in one concerning ring and then compares absolute error between the corresponding rings in the generation and target and, at last, the summation of the errors.

The mass conservation error, i.e., the loss term from E_MC, is defined as

where x and y are the gray values in generated images and average targets, respectively. In addition, m is the number of concerning rings and n is the number of data points in one concerning ring in the matrix. $E$ is the error between the output and average matrix and defined as

where N is the number of all the data points in the matrix. To improve the generation randomness and in view of the error caused by light transmission and reflection, the loss threshold $Ethr$ is introduced herein. Once the loss is less than the threshold, this loss term will be ignored.

The present evaluator is composed of two parts: one is the theoretical model of the spray angle and the other is a CNN encoder to estimate the spray angles from the generated field solutions.

The schematic diagram of the theoretical model of the spray angle is shown in Fig. 6. Several basic hypotheses must be declared before carrying out the theoretical analysis: (a) An element of fluid emerging from the jet exit is assumed to have the constant length and width equal to the jet exit length L_open and width W, as it moves along the trajectory; (b) liquid jet deformation, evaporation, and droplet dispersion are ignored; (c) the fluid element has a constant y-component velocity and an initial angle equal to the central propellant deflection angle φ; (d) the spray angle is assumed to be equal to the slope of the liquid jet at the position where it passes through the gas film; and (e) surface tension, gravity, friction, heat transfer, and phase change are ignored.

Following the definition in Ref. 40, the aerodynamic drag F_d of the liquid element in the x-direction is

where C_d is the drag coefficient. v_g and v_l are the gas velocity and initial liquid element velocity, respectively. According to Newton’s second law in the x-direction, there is

where m_l and a are the mass and acceleration of the liquid element, respectively. v_x is the x-component of the liquid element velocity. It can be rewritten as

Then, the integration of Eq. (7) with respect to time is

and then,

Since v_x(0) = v_lcosφ, there is

As v_x = dx/dt, the integration with respect to time is

where x is the x-coordinate of the element’s trajectory. Since the y-component velocity remains constant and v_y = v_lsinφ = y/t, where y is the y-coordinate, the mathematical expression of the element trajectory is derived as

For the collision between the gas sheet and the rectangular liquid jet, the momentum ratio is

Therefore, Eq. (12) is expressed in terms of the momentum ratio C_TMR as follows:

At the position of penetrating through the gas sheet where y = H, the trajectory expression is modified to

and the slope of the liquid jet θ is obtained as

The theoretical model assumes that the liquid jet does not deform, but in reality, it will deform under aerodynamic forces, which results in a reduction of the effective momentum of the liquid jet. Consequently, the liquid jet deformation factor γ, which is obtained through the experimental results, is introduced to modify the spray angle theoretical model. In this way, Eq. (16) is modified to

Note that the spray angle model derived above has been demonstrated by the experimental results. In fact, when there is no central propellant deflection and φ = 90°, Eq. (17) could be simplified as $θ=γ90◦−arctan(Cd/4CTMR)$⁠. It formally corresponds to the experimental fitting model θ = C₁ arctan(C₂C_TMR) of the liquid–liquid pintle injector in Ref. 41, where C₁ and C₂ are the fitting parameters.

In the field of medical image analysis, the machine learning approach, especially the deep neural networks, has been employed for automated scoliosis assessment.^42–44 In these publications, the x-ray images are fed into the neural network estimator and the spinal Cobb angles are obtained. Similarly, inside the E_SA, there is a well-trained spray angle estimator to output the angle values from the predictive images. The architecture of this down-sampling CNN is like D, except the addition of one fully connected layer in the end to output the estimated spray angle θ′.

The loss term from E_SA is calculated as

Figure 7 shows the generated results in one typical operating condition of the CNN generator and original GAN, which consist of only G and D. The L1 loss used by the CNN generator compares the difference between the generations and targets. This absolute error loss performs very well in some steady or mean state field prediction tasks, such as the work in Refs. 21 and 27. However, when the training cases have a multi-modal distribution, this loss will fail. In our spray field prediction task, although the morphology under one specific operating condition is similar, the detailed droplets are distinguishable. Therefore, the instantaneous spray field solution has many possibilities, which are like the various frames at different times. However, the L1 loss averages all the possibilities and produces a very blurry average image instead. However, the discriminator in GANs, which can be regarded as the loss of generator, is not an explicit loss function. Instead of the pixel-wise loss, D is an approximation loss, which discriminates between the real and fake data distributions and guides the spray generation with detailed morphology.

However, in the training process of GANs, the generator and discriminator have to be balanced trained and the convergence is often an unstable state. The discriminator and generator are always in a seesaw battle to undercut each other. Similar to the work in Ref. 45, we maintain a dynamic ratio between the number of gradient descent iterations on the discriminator and the generator using Exponential Moving Average (EMA).

For the spray simulation task, the discriminator has difficulties capturing the detailed feature of the small droplets. The $LD$ has a possibility of becoming less meaningful through the training process. G will update itself based on the random feedback and the quality of generation may collapse. The G outputs low-quality images through many epochs, and some of them show faint spray patterns in the background but are easily identified as fake. It will be very easy for the discriminator to distinguish the targets and generation so the values of the loss from D drop to zero rapidly.

The comparison of generated spray images from different frameworks demonstrates the superior performance of the mass conservation evaluator as shown in Fig. 8. In some generated images, the background does not agree with the real target; the introduction of $LEMC$ helps G identify the position and intensity of the droplet as well as the background shadow.

In the machine learning field, the parallel training usually is conducted with a batch form. Inputting a certain number of samples to the neural networks in a batch, the solutions are parallel generated and compared to the targets and an error is calculated. In our framework, when obtaining the predicted spray field, all the outputs, ten generated images, in one training batch will be averaged to one image. Then, this average image will be used to feed the down-sampling CNN to obtain the spray angle of this batch. To validate the spray angle estimator, we use the average images from experiments as test samples to output the angle values. Table II shows the comparison of spray angles obtained by Manual Measurements (MM) and CNN estimator. With the increase of C_TMR, the deviations between the two tend to be smaller. The error of all the test cases is less than 5.8%.

After this CNN estimator is trained, it only takes some milliseconds to output an angle value, which is according to the operating conditions of this batch. Then, the estimated value is employed in $LESA$ to update G. Due to the quick decrease at the beginning of the training, this error no longer affects G, only except for the abnormal generated solutions with angles diverged awfully from the theoretical model.

As shown in Fig. 9, the results predicted by the model without E_SA have large value ranges. The difference between maximum and minimum angles in a small momentum ratio level, C_TMR ∈ [0.78, 0.80], even approach 12°. Meantime, the predicted angles have an “even-bias.” Compared with the experimental measurements, all values are closer to the mean value of the whole domain. For small momentum ratios, as shown in Fig. 9(a), most of the predicted angle values are larger than the experimental results. For large momentum ratios, as shown in Fig. 9(b), most of the values are smaller than the experimental results. With the help of E_SA, this phenomenon is reduced and all three groups with different momentum ratios have the predicted angles closer to the experimental results. In addition, the predicted values of spray angles have narrower ranges.

The literature showed that the macroscopic morphology study is important to characterize a spray.^46,47 Here, the simulated spray morphology is analyzed and compared with the experimental results.

Figure 10 compares the simulated and experimental spray morphology under different operating conditions. The typical spray morphology, i.e., liquid column formation, breakup of the column, and lateral expansion of the spray, can be clearly noted. The generated spray field shows that the droplets experience a size reduction before approaching a uniform tiny size distribution because of the impingement and collision. As shown in Fig. 10, the generated spray is fairly comparable with the experimental results; meanwhile, the hollow-cone-shaped profiles with a specific spray angle are well reproduced. For those cases that are out of the training domain, the last column in the figure, the generation also presents a similar good quality compared with those insides. Imperfectly, the background still represents grainy and the values of adjacent data points are not as continuous as those in the real images. However, all in all, the simulation succeeds to report the macroscopical morphology of the spray.

Figure 11 shows the curves of the spray angle vs momentum ratio at different throttling levels; it can be observed that the GAN-PE results coincide with the experimental results well in a wide momentum ratio range. According to theoretical analysis that has been explained in Sec. III C, the spray angle is mainly determined by the momentum ratio, and the simulated solutions also represent this. Due to the difficulties of E_SA to estimate the small angles, the predictions of the cases with small momentum ratios have relatively large errors. The values of different throttling levels are close to each other. For the test cases that are not in the learning domain, the results have a small deviation and all the predicted angle values are less than the manually measured ones. It is because these test cases are not constrained by the targets so the prediction has a trend to approach the mean value of the adjacent operating points, which is happened to be larger.

In this paper, we proposed a novel deep learning framework constrained by physical evaluators to directly predict spray solutions based on generative adversarial networks. The normal discriminator and the mass conservation and spray angle evaluators are used to constrain the CNN to generate the spray solution, including macroscopical morphology and spray angle. The former evaluator is able to improve the training convergence and the latter one helps to obtain more accurate spray angles that are consistent with the operating conditions. It is noteworthy that the related network architecture and spray problem are generic and the proposed framework is potentially suitable for other fluid field simulations that have proper prior physics knowledge. Further research will be carried out for spray droplet size analysis and prediction with the present network framework.

This article was supported by the TUM Open Access Publishing Fund. Hao Ma was supported by the China Scholarship Council (Grant No. 201703170250). Chi Zhang was supported by Deutsche Forschungsgemeinschaft (Grant Nos. DFG HU1527/6-1 and DFG HU1527/10-1). The authors thank Xiangyu Hu, Nils Thuerey, Yujie Zhu, and Wei Wang for beneficial discussions and Anlong Yang for the experimental facilities.

The data that support the findings of this study are available from the corresponding author upon reasonable request.