Modeling and Dynamic Radar Cross-Section Estimation of Chaff Clouds for Real-Time Simulation

Abstract

Chaff is a passive jammer widely used to disrupt radar or radio-frequency sensors. A mass of chaff fibers dispersed in the air is commonly referred to as a chaff cloud. It is nearly impossible to numerically simulate in real-time the enormous amount of chaff fibers composing the chaff cloud. In this paper, we model the behavior of numerically estimated chaff clouds as probability density functions (PDFs) and apply approximation techniques to estimate the radar cross-section (RCS) of the chaff cloud in real time. To model the aerodynamics of the chaff cloud, we represented the combination of PDFs as functions of time and wind speed. The applied approximation techniques—vector radiative transfer and generalized equivalent conductor method—showed a computation time that cannot be achieved by low-frequency methods such as the method of moments or finite-difference time-domain. Moreover, the dynamic RCS results of the approximation techniques showed a similar trend to those of other studies simulating similar situations. The proposed scheme is effective for real-time chaff cloud simulation, and the modeled dynamics and estimated dynamic RCSs can be a standard baseline for developing new analysis methods for chaff clouds. In the future, the proposed scheme will extend to more chaff fibers and more diverse environmental parameters.

1. Introduction

The principle of chaff as a radar countermeasure is to generate false targets on the radar by reflecting radar signals from a large amount of chaff. Chaff is typically depicted as a glass fiber coated with a metallic material such as aluminum and has a length of half a wavelength of radar waves to reflect them effectively. Sometimes, different lengths of chaff clusters are used at once to cover a wider frequency band. Chaff uses a large amount of chaff fibers for single use, and the chaff fibers spread in the air is called a chaff cloud. Due to being cost-effective and efficient compared to other radar countermeasures, chaff is regarded as an essential part of military vehicles.

The methods to disperse chaff fibers have evolved over the years, from basic ejection out of airplane windows to launching with spring-loaded or pneumatic machines. Currently, the methods can be categorized into (1) pyrotechnic charges, (2) rockets, (3) mortars, (4) air flows, and (5) motors. Chaff is ejected either mechanically or pyrotechnically. Mechanical ejection uses small, foil laminated cardboard boxes (2.8 by 4.8 by 0.8 inches) that are torn open during ejection. Debris from the cardboard boxes consists of the opened box, two high-impact polystyrene plastic support pieces (2.75 by 4.75 by 0.06 inches), and paper wrapping for each dipole cut. Cardboard specifications have changed from virgin kraft paper to recycled kraft paper, due to its biodegradable efficiency. The sealing adhesive for these boxes is aqueous-type polyvinyl acetate and a felt spacer. Pyrotechnic ejection uses hot gases generated by an explosive impulse cartridge. The gases push a small plastic piston down a chaff-filled tube with an 8-inch long and 1-inch square cross-section. This ejects a small plastic end cap, followed by the chaff fibers, and the tube remains in the aircraft. The debris that is ejected consists of two 1-inch square pieces of plastic l/8-inch thick (the piston and the end cap) [1].

In previous studies, various methods have been introduced to model the spatiotemporal distribution of chaff clouds. The most common method is simply to use a single probability distribution. For example, Bendayan et al. [2], Pandey [3], Yanchun et al. [4], and Zuo et al. [5] assumed the chaff clouds have a spherical shape and the chaff fibers are uniformly distributed, while Pinchot et al. [6] and Chae et al. [7] assumed a spherical shape and normal distribution. However, using only a single probability distribution does not reflect the dynamic characteristics of chaff, making it difficult to represent the distribution of chaff clouds in real situations. To realistically model the distribution of chaff clouds, numerical methods have been introduced. Zhu et al. [8] calculated the equation of motion for chaff by considering the initial velocity due to an explosion, while Huang et al. [9] and Wang et al. [10] calculated the equation of motion for chaff by considering the aerodynamic interference among chaff fibers using computational fluid dynamics (CFDs). These numerical methods represent the distribution of chaff clouds more realistically than methods using probability distributions but still do not reflect the impact of collisions among chaff fibers. By considering not only the impact of explosions and wind but also collisions among chaff fibers, in this paper, we model the spatiotemporal distribution of chaff clouds using the computational fluid dynamics and discrete element method (CFD-DEM) one-way coupling model that integrates CFDs and DEMs. In addition, we finally model the distribution of chaff clouds in various situations using probability density functions (PDFs) since the CFD-DEM analysis requires a large number of computer resources and computation time.

There are three main methods for estimating the RCSs of chaff clouds: empirical, statistical, and numerical methods. The empirical method is an important method that can consider all the complex factors that theoretical calculations cannot take into account by directly measuring the chaff cloud. Statistical methods include a method of statistically modeling numerical analysis or experimental results and a method to obtain the total scattered wave of a chaff cloud by multiplying the average RCS of the single chaff by the number of chaffs [11]. Methods to numerically estimate the radar cross-section (RCS) of a chaff cloud can be categorized into low-frequency methods and approximation methods. The method of moments (MoM) is one of the low-frequency numerical methods that provides high accuracy in estimating the RCS of a chaff cloud and has been used for this purpose in many works of literature [12,13]. However, estimating the average RCS of a chaff cloud using the MoM requires a large number of realizations, and moreover, the number of calculable chaff fibers does not exceed 2000. This is why the MoM is commonly only used to validate new methods for a small number of chaff fibers [14,15]. Due to these limitations, various approximation methods that have the advantages of small computer resources and a fast computation time have been attempted to estimate the average RCS of a large number of chaff fibers. Marcus [16] proposed the equivalent conductor (EC) method to calculate the average RCS of a conducting continuum by replacing the chaff cloud with it. However, the EC method is not applicable when the orientation of the chaff fibers is not uniformly distributed, so Seo et al. [17] proposed the generalized EC (GEC) method that can be applied to chaff clouds with arbitrary orientation distributions. Zuo et al. [18] estimated the average RCS of a chaff cloud composed of a large number of chaff fibers using the vector radiative transfer (VRT) method, which was aimed at chaff clouds with arbitrary orientations and densities.

In this paper, the aerodynamic behavior of a chaff cloud is estimated using the CFD-DEM to obtain a chaff cloud model that can correspond to various practical variables, and the aerodynamic behaviors are fitted into PDFs. The dynamic RCS of the PDF-modeled chaff cloud is estimated through the GEC and VRT methods.

2. Materials and Methods

2.1. Modeling of the Chaff Cloud

In order to model the chaff cloud, we use the CFD-DEM one-way coupling technique that calculates the aerodynamic behavior of each chaff. It is computationally expensive due to the consideration of factors such as drag, gravity, and collisions, and is only suitable for simulating specific scenarios, making it inappropriate for real-time chaff cloud simulation. Therefore, in this paper, we model the position and orientation distributions using the Gaussian distribution sum (GDS) and Gaussian mixture model (GMM) based on the results of the CFD-DEM. This method can generate distributions by adjusting parameters with different wind speeds, rather than the wind speed used for PDF model fitting. Using the CFD-DEM, we simulate the dynamic behavior of the chaff cloud for 5 s under seven different wind speeds and fitted the resulting distribution with statistical models.

2.1.1. CFD-DEM One-Way Coupling Technique

In the CFD-DEM technique for a chaff cloud, the CFD considers the effects of explosions and winds that affect the diffusion and distribution of chaff fibers, while the DEM treats chaff fibers as particles and calculates the aerodynamic behavior of each chaff particle. The CFD is the most representative computational technique for solving fluid mechanics problems, and it discretizes the governing equations of fluid flow, the Navier–Stokes equations, using finite volume methods for numerical calculations. In addition, the DEM is a computational technique suitable for systems composed of many particles, and it represents each particle as a node, assigning six degrees of freedom to each node to calculate the particle’s behavior. Additionally, collision models can be introduced to calculate collisions among particles. We calculate the actual behavior of chaff fibers by considering both the factors affecting the fluid, such as explosions and winds, and the collisions between chaff particles, by reflecting the fluid flow calculated by the CFD in the DEM analysis, and obtained the behavior of chaff clouds, which is modeled as a PDF in Section 2.1.2. The program used for the CFD-DEM technique is Star-CCM+ version 2020.3.

2.1.2. Statistical Modeling of Chaff Clouds

The GDS is a useful method when the distribution of data is clearly distinguished, allowing the probability density of the overall distribution to be obtained by simply summing the distinct data components. A single component of the GDS is given by:

$$N(x;{\mu }_k,{\sigma }_k^2)=\frac{1}{\sqrt{2\pi {\sigma }_k^2}}\exp \left(-\frac{{(x-{\mu }_k)}^2}{2{\sigma }_k^2}\right),$$

(1)

and the PDFs of a GDS with K components can be written by:

$$f_{GDS}(x)=\frac{1}{K}{\displaystyle \sum _{k=1}^{K}N(x;{\mu }_k,{\sigma }_k^2)},$$

(2)

where K is the number of the components, subscript k denotes the order of the components, μ is the mean, and σ is the standard deviation.

On the other hand, the GMM is a useful method when the distribution of data is not clearly separated, and additional consideration of covariance and mixing coefficients between Gaussian distributions is required. A single component of the GMM is:

$$N(x;{\overline{\mu }}_k,C_k)=\frac{1}{{(2\pi )}^{D/2}{\left|C_k\right|}^{1/2}}\exp \left(-\frac{1}{2}{(x-{\overline{\mu }}_k)}^T{C}_k^{-1}(x-{\overline{\mu }}_k)\right),$$

(3)

then, the GMM with K components has a PDF form [19] as follows:

$$f_{GMM}(x)={\displaystyle \sum _{k=1}^K{\pi }_{k}N(x;{\overline{\mu }}_k,C_k)},$$

(4)

where D is the dimension of the data, π is the mixing coefficient, $\overline{\mu }$ is the mean vector, and C is the covariance matrix.

We analyze the spatiotemporal distribution of 2124 chaff fibers assumed to be used in a warship during 5 s using the CFD-DEM and model them as PDFs in Cartesian coordinates, zenith, and azimuth angles. For a warship chaff, it is assumed that the blast pressure inside the chaff cartridge is 1 GPa, and the wind is blowing at the wind speed of 0, 2.5, 5, 7.5, 10, 20, and 30 m/s in the X-direction. If we simulate the aerodynamic behavior of the chaff cloud by adjusting the blast pressure, initial chaff cartridge angle, and other parameters, we can obtain a PDF model that includes those parameters. The PDFs of the X-axis and Z-axis distributions are modeled using the GDS with three and two components, respectively, while the GMM is used to model the PDFs of the Y-axis. The zenith and azimuth angles of the chaff are modeled as a Gaussian distribution and a uniform distribution in degrees, respectively. The fitting parameter values are hand-tuned using a data analysis program, and the parameter values of the fitted functions are summarized in

Table 1. The parameter values of the PDF model.

Category		Average	Variance (GDS)/Covariance(GMM)	Mixing Rate
X	1	$vt-0.09v-0.195$	$(97.2v^{-0.2}-72v^{-0.4}·{0.3}^t)·{10}^{-4}$	1/3
	2	$vt-0.0899v-0.028$	$(75.6v^{-0.07}-54v^{-0.2}·{0.3}^t)·{10}^{-4}$	1/3
	3	$vt-0.0899v+0.146$	$(57.6v^{0.04}-36v^{-0.08}·{0.3}^t)·{10}^{-4}$	1/3
Y	1	$-0.252$	$\left(216+1.08v\right)·t^{0.01+0.03·{0.9}^v}·{10}^{-4}$	1/2
	2	$0.252$		1/2
Z	1	$-0.1104-0.243525t$	$324+\frac{216{·10}^{-4}-324}{1+{\left(\frac{t}{239}\right)}^{(2-0.013v)}}$	1/2
	2	$0.3704-0.213475t$	$2.52+\frac{144{·10}^{-4}-2.52}{1+{\left(\frac{t}{13}\right)}^{(2.3+0.01v)}}$	1/2
$\theta$		$90$	$1414·t^{0.1}$	1
$\varphi$		Uniform distribution

, in which v represents wind speed, and t represents the elapsed time after the chaff cartridge explosion.

The similarities of the PDF model can be evaluated with the Bhattacharyya coefficient, in which the Bhattacharyya distance D is commonly used to evaluate the accuracy of statistical distribution to represent experimental data [20] and is given by:

$$D=-\ln \left(\rho \right),$$

(5)

$$\rho ={\displaystyle \sum _{x\in \chi }\sqrt{P\left(x\right)Q\left(x\right)}},$$

(6)

where ρ is the Bhattacharyya coefficient, P and Q are the probability distributions, x is the random variable, and $\chi$ is the domain of x. The Bhattacharyya distance D ranges from zero to infinity, and it is used to evaluate which PDF model is more suitable for evaluating the order of similarity to the measured values among multiple PDF models. Therefore, in this paper, we compare the distribution of chaff cloud data interpreted using the CFD-DEM to the PDF model using the Bhattacharyya coefficient to evaluate the suitability of the PDF model we created. The Bhattacharyya coefficient has a range of zero to one, and the closer it is to one, the more similar the distribution of the measured values and the PDF model is.

2.2. RCS Estimation of the Chaff Cloud

In this paper, the aerodynamic behavior of the chaff cloud is represented using the PDF model for real-time simulation. However, in practical situations, low-frequency numerical methods such as the MoM and finite-difference time-domain (FDTD) are unable to calculate these massive chaff fibers, as a chaff cloud consists of millions of chaff fibers. The approximation methods such as the GEC and VRT do not provide an exact solution, but they are sufficient for obtaining an approximate average RCS of a chaff cloud consisting of large quantities of chaff fibers.

2.2.1. Generalized Equivalent Conductor (GEC) Method

The EC and GEC are methods for calculating the electromagnetic characteristics of the entire scattering body by replacing the chaff cloud formed by the dispersion of the chaff fibers with an effective medium [21]. Especially, the GEC can be applied to chaff clouds with arbitrary orientation distributions. Figure 1 shows the brief methodology of a GEC. The blue points are the chaff fibers and the red dash lines are the location of the effective medium. By replacing a chaff cloud consisting of numerous chaff fibers with a medium with an effective permittivity, parameters such as effective permittivity, volume, and area can be obtained, and the RCS of the replaced chaff cloud can be estimated. It is not possible for the GEC method to obtain exactly the same results as those obtained from the MoM. Nevertheless, the GEC can provide good accuracy with a very low computation burden [17].

The effective permittivity of the effective medium is [17]:

$${\overline{\overline{\epsilon }}}_{eff}=\left[\begin{array}{cc}{\epsilon }_{\theta \theta }& {\epsilon }_{\theta \varphi }\\ {\epsilon }_{\varphi \theta }& {\epsilon }_{\varphi \varphi }\end{array}\right]=\left[\begin{array}{cc}{\epsilon }_0\left(1-j\frac{{\eta }_0}{2\pi }\cdot \rho {\lambda }^3\cdot \frac{l}{\lambda }\cdot {\overline{I}}_{\theta \theta }\right)& {\epsilon }_0\left(1-j\frac{{\eta }_0}{2\pi }\cdot \rho {\lambda }^3\cdot \frac{l}{\lambda }\cdot {\overline{I}}_{\theta \varphi }\right)\\ {\epsilon }_0\left(1-j\frac{{\eta }_0}{2\pi }\cdot \rho {\lambda }^3\cdot \frac{l}{\lambda }\cdot {\overline{I}}_{\varphi \theta }\right)& {\epsilon }_0\left(1-j\frac{{\eta }_0}{2\pi }\cdot \rho {\lambda }^3\cdot \frac{l}{\lambda }\cdot {\overline{I}}_{\varphi \varphi }\right)\end{array}\right],$$

(7)

where:

$${\overline{I}}_{pq}({\theta }_i,{\varphi }_i)=\frac{{\displaystyle {\int }_{\mathsf{\Omega }}{\tilde{I}}_{pq}\left(\theta ,\varphi ;{\theta }_i,{\varphi }_i\right)W\left(\theta ,\varphi \right)d\theta d\varphi }}{{\displaystyle {\int }_{\mathsf{\Omega }}W(\theta ,\varphi )d\theta d\varphi }}$$

(8)

is the effective average current, the subscript p and q are the vertical or horizontal polarization, and ${\tilde{\mathbf{I}}}_{pq}(\theta ,\varphi ;{\theta }_i,{\varphi }_i)$ is the actual current component of a single chaff with orientation (θ, ϕ) to an incident wave with direction (θ_i, ϕ_i). The actual current component can be obtained with the MoM. The parameters ${\epsilon }_0$, ${\eta }_0$, ρ, l, and λ denote the electric constant, the intrinsic impedance of free space, the average density of chaff fibers, the length of chaff fibers, and the wavelength in free space, respectively [22]. The effective average current is calculated based on the weight function W(θ, ϕ), which represents the orientation distribution of chaff.

The average RCS of the effective medium is obtained as the sum of incoherent RCS [23] and coherent RCS [22], as follows:

$${\sigma }_{total}={\sigma }_{incoh}+{\sigma }_{coh}.$$

(9)

2.2.2. Vector Radiative Transfer (VRT) Method

The theory of radiative transfer deals with the transfer of energy as it passes through the spaces where particles are distributed. It was initially developed in the field of astrophysics and is now also used for the remote sensing of terrains or scattering media. A chaff cloud is a space distributed with a large number of chaff fibers, and it can be applied to calculate changes in energy that are covered in radiative transfer theory. The VRT method, which is an extension of the theory in the vector case to consider polarization, can well reflect the changes in polarization components due to the shape of the chaff. When specific intensity is incident on a certain space, the change in specific intensity passing through that space is represented by:

$$\frac{d\overline{I}(\overline{r},\widehat{s})}{ds}=-{\overline{\overline{\kappa }}}_{\epsilon }\overline{I}(\overline{r},\widehat{s})-{\kappa }_{ag}\overline{I}(\overline{r},\widehat{s})+{\displaystyle {\int }_{4\pi }^{}d{\mathsf{\Omega }}^{\prime }\overline{\overline{P}}(\widehat{s},{\widehat{s}}^{\prime })\overline{I}(\overline{r},{\widehat{s}}^{\prime })}+{\overline{J}}_e,$$

(10)

where $\overline{r}$ is a position vector, $\widehat{s}$ denotes the direction of propagation, $\overline{I}$ is the Stokes parameters of specific intensity, ${\stackrel{=}{\kappa }}_{\epsilon }$ is the extinction matrix, ${\kappa }_{ag}$ is the power absorption coefficient of the background material, $\stackrel{=}{P}$ is the phase matrix, and ${\overline{J}}_e$ is the emission vector.

The radiative transfer equation is formulated in terms of three constitutive functions: (1) the extinction matrix, which describes the attenuation of the specific intensity due to absorption and scattering; (2) the phase matrix, which characterizes the coupling between the incident and scattered intensities at every point in the medium; and (3) the emission function, which accounts for the black body emission by the medium [24]. However, since the background material of a chaff cloud is air and the emission term is also negligible in radar remote sensing, we can neglect the loss term and the emission term in the VRT equation, then the VRT equation for a chaff cloud can be rewritten as:

$$\frac{d\overline{I}(\overline{r},\widehat{s})}{ds}=-{\overline{\overline{\kappa }}}_{\epsilon }\overline{I}(\overline{r},\widehat{s})+{\displaystyle {\int }_{4\pi }^{}d{\mathsf{\Omega }}^{\prime }\overline{\overline{P}}(\widehat{s},{\widehat{s}}^{\prime })\overline{I}(\overline{r},{\widehat{s}}^{\prime })}.$$

(11)

We can obtain the solution of the VRT Equation (1) using an iterative method. The procedure of the iterative method begins with the computation of the zeroth-order solution. If the zeroth-order solution is obtained, then it can be used as the source function to calculate the first-order solution. In this work, we use solutions up to the first order, and these solutions are used to calculate the transformation matrix that relates the Stokes parameters of the incident and scattered waves. The relationship between the incident and scattered waves is:

$${\overline{I}}_s({\theta }_s,{\varphi }_s)={\overline{\overline{T}}}_t({\theta }_s,{\varphi }_s;{\theta }_i,{\varphi }_i){\overline{I}}_i({\theta }_i,{\varphi }_i),$$

(12)

where ${\stackrel{=}{I}}_s$ and ${\stackrel{=}{I}}_i$ are the scattered and incident waves, respectively, and ${\stackrel{=}{T}}_t$ is the total transformation matrix and is expressed as the sum of the zeroth- and first-order transformation matrices:

$${\overline{\overline{T}}}_t({\theta }_s,{\varphi }_s;{\theta }_i,{\varphi }_i)={\overline{\overline{T}}}_0({\theta }_s,{\varphi }_s;{\theta }_i,{\varphi }_i)+{\overline{\overline{T}}}_1({\theta }_s,{\varphi }_s;{\theta }_i,{\varphi }_i)$$

(13)

The zeroth- and first-order transformation matrices are given by [24]:

$${\overline{\overline{T}}}_0({\theta }_s,{\varphi }_s;{\theta }_i,{\varphi }_i)=e^{-{\overline{\overline{\kappa }}}_{\epsilon }^{s}d/\cos {\theta }_s}\overline{\overline{R}}({\theta }_i,{\varphi }_i)e^{-{\overline{\overline{\kappa }}}_{\epsilon }^{i}d/\cos {\theta }_i},$$

(14)

$$\begin{array}{cc}{\overline{\overline{T}}}_1\left({\theta }_s,{\varphi }_s;{\theta }_i,{\varphi }_i\right)& =\sec {\theta }_i\overline{\overline{E}}\left({\theta }_s,{\varphi }_s\right)\overline{\overline{D}}\left(-\beta \left({\theta }_s,{\varphi }_s\right)d\sec {\theta }_s\right){\overline{\overline{E}}}^{-1}\left({\theta }_s,{\varphi }_s\right)\overline{\overline{R}}\left({\theta }_i,{\varphi }_s\right)\overline{\overline{E}}\left(\pi -{\theta }_i,{\varphi }_s\right){\overline{\overline{A}}}_1\\ & ·{\overline{\overline{E}}}^{-1}\left({\theta }_s,{\varphi }_i\right)\overline{\overline{R}}\left({\theta }_i,{\varphi }_i\right)\overline{\overline{E}}\left(\pi -{\theta }_i,{\varphi }_i\right)\overline{\overline{D}}\left(-\beta \left(\pi -{\theta }_i,{\varphi }_i\right)d\sec {\theta }_i\right){\overline{\overline{E}}}^{-1}\left(\pi -{\theta }_i,{\varphi }_i\right)\\ & +\sec {\theta }_i\overline{\overline{E}}\left({\theta }_s,{\varphi }_s\right){\overline{\overline{A}}}_2{\overline{\overline{E}}}^{-1}\left({\theta }_i,{\varphi }_i\right)\overline{\overline{R}}\left({\theta }_i,{\varphi }_i\right)\overline{\overline{E}}\left(\pi -{\theta }_i,{\varphi }_i\right)\overline{\overline{D}}\left(-\beta \left(\pi -{\theta }_i,{\varphi }_i\right)d\sec {\theta }_i\right){\overline{\overline{E}}}^{-1}\left(\pi -{\theta }_i,{\varphi }_i\right)\\ & +\sec {\theta }_i\overline{\overline{E}}\left({\theta }_s,{\varphi }_s\right)\overline{\overline{D}}\left(-\beta \left({\theta }_s,{\varphi }_s\right)d\sec {\theta }_s\right){\overline{\overline{E}}}^{-1}\left({\theta }_s,{\varphi }_s\right)\overline{\overline{R}}\left({\theta }_i,{\varphi }_s\right)\overline{\overline{E}}\left(\pi -{\theta }_i,{\varphi }_s\right){\overline{\overline{A}}}_3{\overline{\overline{E}}}^{-1}\left(\pi -{\theta }_i,{\varphi }_i\right)\\ & +\sec {\theta }_i\overline{\overline{E}}\left({\theta }_s,{\varphi }_s\right){\overline{\overline{A}}}_4{\overline{\overline{E}}}^{-1}\left(\pi -{\theta }_i,{\varphi }_i\right),\end{array}$$

(15)

In Figure 2, the arrow marked ‘a’ represents the zeroth-order transformation, while the remaining arrows marked ‘b–e’ denote the first-order transformation matrix. The arrows in Figure 2 represent the transfer of energy, and the small dots are particles in the medium. Because the chaff cloud is suspended in the air, ground term $\stackrel{=}{R}$ should be excluded from (14) and (15). This allows (13) to be rewritten as:

$${\overline{\overline{T}}}_t\left({\theta }_s,{\varphi }_s;{\theta }_i,{\varphi }_i\right)=\sec {\theta }_i\overline{\overline{E}}\left({\theta }_s,{\varphi }_s\right){\overline{\overline{A}}}_4{\overline{\overline{E}}}^{-1}\left(\pi -{\theta }_i,{\varphi }_i\right),$$

(16)

where θ and ϕ are the zenith and azimuth angle of the specific intensity, respectively, the subscript s and I represent scattering and incidence, respectively, and:

$$\begin{array}{l}{\left[{\overline{\overline{A}}}_4\left({\theta }_s,{\varphi }_s;{\theta }_i,{\varphi }_i\right)\right]}_{jk}=\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\frac{1-e^{-{\beta }_j\left({\theta }_s,{\varphi }_s\right)d\sec {\theta }_s-{\beta }_k\left(\pi -{\theta }_i,{\varphi }_i\right)d\sec {\theta }_i}}{{\beta }_j\left({\theta }_s,{\varphi }_s\right)\sec {\theta }_s+{\beta }_k\left(\pi -{\theta }_i,{\varphi }_i\right)\sec {\theta }_i}{\left[{\overline{\overline{E}}}^{-1}\left({\theta }_s,{\varphi }_s\right)\overline{\overline{P}}\left({\theta }_s,{\varphi }_s;\pi -{\theta }_i,{\varphi }_i\right)\overline{\overline{E}}\left(\pi -{\theta }_i,{\varphi }_i\right)\right]}_{jk}.\end{array}$$

(17)

Here, $\stackrel{=}{E}(\theta ,\varphi )$ is the eigenmatrix and β_j is the jth eigenvalue of the extinction matrix.

2.3. Simulation Parameters

For practical chaff simulation, it is assumed that a typical Short-Range Air Defense (SHORAD) radar detects a chaff cloud launched from a warship and the SHORAD radar has a beam width of 2° elevation ang azimuth angles with a pulse length (or gate range) of 400 ns at a typical range of search targets of 6 km [1]. In addition, it is assumed that the operating frequency of the SHORAD radar is in a Ku band, the angular resolution including elevation is about 210 m, and the range resolution is 60 m. The radar resolution volume is shown in Figure 3.

For chaff fibers, it is assumed that a bunch of 2124 chaff fibers in a chaff cartridge are launched from a warship and exploded at a height of 35 m above sea level, forming a chaff cloud that falls to the sea level. Because the height of a ship depends on its type, speed, course, size, and naval environment, the typical height at which chaff cartridges detonate is not fixed, so we assumed a height of 35 m, which corresponds to being slightly higher than the ship’s height and the velocity of the chaff cartridge is ignored. On the other hand, the length of each chaff fiber is a half wavelength of the operating frequency in the Ku band, and the diameter of each chaff fiber is 25 μm. A total of 2124 chaff fibers are used to generate the chaff cloud, which is subject to movement by wind. The wind speeds acting on the chaff cloud are 0, 2.5, 5, 7.5, 10, 20, and 30 m/s in the +x direction. Therefore, the chaff cloud can move out of the radar resolution volume due to wind, and the RCS estimation of the chaff cloud participating in a single radar resolution volume is performed. The simulation was conducted using MATLAB R2021b on a Windows 10 PC with an Intel i9-11900 CPU and 64 GB of RAM. The simulation parameters are summarized in

Table 2. The simulation parameters.

Category	Simulation Parameter	Value
Radar specifications	Frequency	Ku band
	Pulse length	400 ns
	Beam width	2 degrees
	Range of search targets	6 km
Simulation setup	Wind speed	0, 2.5, 5, 7.5, 10, 20, 30 m/s
	Simulation interval	0.1 s
	Total simulation time	5 s
	Initial location of chaff cloud	(0, 6000, 35) m 1
Properties of chaff	Length	0.01 m
	Diameter	25 μm
	Quantity	2124
PC specifications	OS	Windows 10
	CPU	Intel i9-11900
	Program	MATLAB R2021b

¹ Coordinates in a Cartesian coordinate system: (X, Y, Z).

.

3. Results

3.1. Model Evaluation and Shape

We modeled the PDF model of the chaff cloud based on interpretation of CFD-DEM and summarized the PDF model in Table 1. In this section, the similarity between our PDF model and the original data is evaluated using the Bhattacharyya coefficient of (6), and the Bhattacharyya coefficients between the datasets are shown in Figure 4. The colors of the lines in Figure 4a–e correspond to the designated wind speeds, and the relation with wind speeds and the colors are shown in Figure 4f. As shown in Figure 4a–e, the Bhattacharyya coefficients are very close to 1 in all distributions of position and orientation, indicating that very similar distributions to the original CFD-DEM data can be obtained. In addition, no particular common trend can be found as time or wind speed vary, and almost constant values are shown. Therefore, our model allows us to obtain the aerodynamic behavior of a chaff cloud similar to those obtained from the CFD-DEM by adjusting the time or wind speed of the PDF model.

The chaff clouds generated by the CFD-DEM and PDF model at several time steps are plotted in Figure 5 and Figure 6, respectively. The chaff cloud simulated using the CFD-DEM has a distinct circular ring shape on the YZ plane at 0.1 s, which gradually breaks down over time. On the other hand, the PDF model exhibited an indistinct rectangular ring shape initially and slowly spread out and eventually became similar in shape to the CFD-DEM chaff cloud. This is because the chaff fibers spread out circularly from the center of the chaff cartridge cross-section, while our PDF model was modeled in each axis direction in a Cartesian coordinate system.

Both datasets exhibit nearly constant diffusion along the X- and Y-axes, while diffusion along the Z-axis is observed due to the relationship between the chaff orientation and drag. This pattern indicates that chaff dispersion by explosion is completed before 0.1 s and thereafter, the chaff spread due to gravity, drag, and collisions.

3.2. RCS Estimation with GEC and VRT

We estimate the average RCS of the chaff cloud generated through our PDF model using the GEC and VRT. Figure 7 represents the estimated dynamic RCS of the chaff cloud for several wind speeds. The wind speeds and colors of the figure are the same as in Figure 4f. For the VV polarization, the dynamic RCSs estimated by the two methods start with similar values but the RCS levels gradually increase, and at 5 s, the difference between the RCSs of the two methods reaches about 1.2 dB. The HH-polarized dynamic RCSs show a slight decrease only at the beginning, but from then on, a constant difference of about 2 dB are observed until 5 s. The reason for these differences is that the higher-order terms in the VRT equation’s first-order iterative solution are ignored, resulting in smaller values compared to the exact numerical solution. The reason why the HH polarized dynamic RCSs show this tendency is that the orientation distribution of zenith angle θ changes relatively quickly from 0.1 to about 1 second, causing the HH polarized RCS to decrease. However, the change in the distribution of θ after that time is small. The increase in the VV polarized RCS is similar to the factor that causes the decrease in the HH polarized RCS, but the rate of change is greater than that of the HH polarized RCS. The orientation distribution of θ along the time is shown in Figure 8. Except for minor differences depending on the wind speed, generally, similar dynamic RCS values are obtained. The abrupt drop in RCSs at a wind speed of 30 m/s is because all or part of the chaff cloud is calculated to be outside the radar resolution volume.

As shown in Figure 9, the average computation time for the average RCS estimation by the chaff cloud at any one moment is approximately 0.04 s for the GEC and 0.11 s for the VRT. The wind speeds and colors of Figure 9 are the same as Figure 4f. Although the GEC and VRT achieve computation times that are impossible to achieve with low-frequency numerical techniques, they are not sufficient to be called real-time when considering the radar’s pulse length. Considering that the chaff cloud moves according to the wind and gravity while almost maintaining the shape of the chaff cloud formed immediately after the explosion, it is applicable to real-time simulation with a fixed step of about 0.1 s. Moreover, the current simulations were performed using MATLAB, but if ported and optimized using low-level languages such as C or C++, better performance can be achieved.

The red line and magenta line in Figure 9 show slightly different patterns, representing wind speeds of 0 m/s and 30 m/s, respectively. The reason why the initial computation time of the red line took longer is due to the simulation order. The simulation proceeds in the order of wind speed from 0 m/s to 30 m/s, causing an increase in the initial computation time of the red line due to the influence of variable declaration and initialization. After 3.7 s of the magenta line, the computation time drops to almost 0 s. This is because the chaff cloud has completely moved out of the radar resolution volume by the effect of the wind.

Figure 10 represents the estimation of a large number of chaff fibers. The number of chaff fibers increases from 40,000 to 2 million. The distribution of chaff fiber positions inside a sphere with a radius of 15 m follows a uniform distribution, and the incident and scattering angles of the electromagnetic wave on the sphere are both 90 degrees. The distribution of the chaff fiber orientations, θ and ϕ, ranges from 0 to 180 degrees and 0 to 360 degrees, respectively, following a uniform distribution. The RCS pattern in Figure 10a shows that the results obtained using VRT are smaller than those obtained using the GEC, as previously estimated. Of particular note is the computational time. Figure 10b represents the results from a subsequent run for accurate comparison. When the simulation time, including the entire process, was measured, the simulation time for chaff cloud analysis using GEC ranged from approximately 0.0501 s to 0.401 s, while for VRT, it ranged from approximately 0.132 s to 0.653 s for a single chaff cloud analysis. However, when measuring the estimation time of GEC and VRT themselves, excluding the process of handling the chaff cloud data, the results showed that the estimation time remained constant regardless of the increasing number of chaff fibers. Ultimately, when calculating for two million chaff fibers, the simulation time for GEC and VRT was 9.931 times and 6.127 times longer, respectively, compared to the self-computation time. These results indicate that both GEC and VRT are capable of handling calculations for a large number of chaff fibers, but they also highlight the importance of efficient processing of large-scale data for real-time simulations.

4. Discussion

We simulated the aerodynamic behavior of the chaff clouds using numerical methods and modeled the results in PDF format to make them suitable for real-time simulation. There are not many studies that have simulated the aerodynamic behavior of chaff clouds over time. If we limit the range to chaff clouds formed using explosives, the number of such studies is even smaller. Few studies clearly present information on the shape, size, and growth of a chaff cloud formed by an explosion, or until it grows.

In [25], a chaff cloud diffusion model was presented for the chaff cloud with a diameter of over 80 m when diffusing for 30 s. Their chaff cloud diffusion model is quite effective for high-density chaff cloud aerodynamics, but it is hard to represent the initial explosion and collision among chaff fibers. In [22], the authors simulated the aerodynamics of a chaff cloud using six degrees of freedom equations and estimated the dynamic RCS of the chaff cloud using the GEC. However, because the chaff fiber shape changes as well as collisions among chaff fibers are not considered, two types of chaff, straight and bent chaff, were analyzed. While the volume of the chaff cloud consisting of the straight chaff fibers kept increasing, the bent chaff fibers stopped diffusing sideways and their chaff cloud maintained a constant size once it reached a certain size. In this work, the shape of the chaff cloud maintained a constant size from the very beginning, as the diffusion caused by the explosion was stopped by air resistance before 0.1 s. Our analysis considered the same factors as previous studies, including air resistance, gravity, and wind, as well as the direct effects of the explosion on the chaff and collisions between chaff particles. Therefore, we cannot determine which chaff cloud is more similar to a real chaff cloud, but one thing is clear: our PDF model considers more variables than previous studies.

In our analysis, the radar looks at a chaff cloud almost from the side. Among the orientation components of the chaff, the zenith angle θ has a significant impact on the RCS due to the linear shape. Zhang et al. [25] did not consider the variation of the orientation distribution to investigate the shielding effect, so after the chaff cloud reached a certain size, there was no change in the RCSs. However, by selecting an orientation distribution more horizontal to the surface, the HH polarized RCS was higher than the VV polarized RCS. Seo et al. [22] showed that the dynamic RCS of the chaff cloud was initially overshooting, and in the case of bent chaff fibers, as time passed, the orientation distribution of the chaff changed to be close to horizontal, resulting in a larger RCS for the HH polarization than the VV polarization. Both papers showed a rapid increase in the RCS as the size of the chaff cloud increased, and represented a larger RCS level for the HH polarization than the VV polarization. Our dynamic RCS results did not show a rapid increase due to diffusion by an explosion ending within a short time of 0.1 s, but as with the previous two papers, a higher RCS was calculated for the HH polarization than the VV polarization, confirming a similar trend to the existing methods when comparing the estimated RCS results in similar situations.

To use the dynamic RCSs of chaff clouds in hardware-in-the-loop simulation (HILS) for RF sensors, two elements are necessary: a chaff cloud behavior model that can consider various environmental variables in real situations, and the ability to calculate the real-time scattering characteristics of the chaff cloud. While the aerodynamic behavior of chaff clouds can be precisely analyzed using numerical methods such as the CFD-DEM and six degrees of freedom equations, the computation time cannot meet the requirements for real-time simulations. Therefore, we modeled the aerodynamic behavior of the chaff cloud as PDFs that can be adjusted for various environmental variables. We estimated the dynamic RCSs of the PDF model using the GEC and VRT as approximate methods in a short amount of computation time. These methods will enable real-time simulations of chaff clouds in HILS. Currently, the PDF model can only consider wind, gravity, explosion, and collision, but we plan to expand it to consider various environmental variables in the future and will release it. Since the PDF model is presented in closed form in this paper, it allows for the easy generation of chaff clouds, and the RCS estimation results of the PDF model can be used as a comparison in future research.

5. Conclusions

In this paper, we proposed a framework for the real-time simulation of chaff clouds. The aerodynamic behavior of the chaff cloud was numerically estimated by considering drag, gravity, explosions, and collisions between chaff fibers, and was fitted to PDF models for application in various situations. The chaff cloud model generated through the PDFs showed a good fit to the original data, and the dynamic RCSs of the model was estimated by the GEC and VRT, which can quickly calculate a large amount of chaff fibers. Although a discrepancy existed due to the limitations of the approximation method, it was found to be able to calculate the dynamic RCSs of the chaff cloud that is unattainable with numerical methods such as the MoM and FDTD. If we consider our chaff cloud model as a single cluster, in reality, there are hundreds to thousands of clusters dispersing in the air. Therefore, the actual size of the chaff cloud is much larger. Through our chaff cloud model, it would be possible to extend it to represent the actual chaff cloud. The chaff cloud’s PDFs, which can be adjusted with a small number of parameters, were used, but obtaining a chaff cloud model that can be utilized in a wider range of situations is possible by adding more parameters and setting more realistic conditions.