Abstract
Although the uniform theory of diffraction (UTD) could be theoretically applied to arbitrarilyshaped convex objects modeled by nonuniform rational B-splines (NURBS), one of the great challenges in calculation of the UTD surface diffracted fields is the difficulty in determining the geodesic paths along which the creeping waves propagate on arbitrarilyshaped NURBS surfaces. In differential geometry, geodesic paths satisfy geodesic differential equation (GDE). Hence, in this paper, a general and efficient adaptive variable step Euler method is introduced for solving the GDE on arbitrarilyshaped NURBS surfaces. In contrast with conventional Euler method, the proposed method employs a shape factor (SF) to efficiently enhance the accuracy of tracing and extends the application of UTD for practical engineering. The validity and usefulness of the algorithm can be verified by the numerical results.
1. Introduction
In the high-frequency electromagnetic problems, the UTD analysis of the diffraction of electromagnetic (EM) waves [1] by smooth convex surfaces is of great importance, such as the EM compatibility, the analysis of coupling between two antennas [2], and the estimation of scattering properties. And, as known to all, ray tracing of creeping waves plays an important role in the process of obtaining the solution of UTD surface diffracted flied. Therefore it is essential to investigate the creeping ray tracing (geodesic paths) on surfaces firstly. In fact, however, except some typically geometric objects which have analytical solution of the geodesic path, it is a great challenge in determining the geodesic paths along which the creeping waves propagate on arbitrarily shaped smooth convex surfaces.
Jha et al. [3–5] put forward a novel technique known as the geodesic constant method, an analytical method, and it is applicable to the creeping ray tracing on general paraboloids of revolution and hybrid surfaces of revolution. Usually, for some slightly more complex models such as simple aircraft, satellite, they are usually approximated by some typically geometric objects, such as boards, cylinders, cones, and spheres, which already have analytical solution in tracing the creeping rays, for engineering application. But it is difficult to use such typically geometric objects to approximate arbitrarily shaped objects, which apparently limits the application of UTD method.
Hence it is very essential to introduce the numerical ray tracing methods of creeping waves to handle arbitrarily shaped convex objects. In [6], the models are characterized by many discrete triangular meshes. And creeping ray tracing numerically on discrete triangular meshes is presented by Surazhsky et al. However, such methods cannot be used directly for the UTD solution. Except modeling arbitrarily shaped convex object by discrete triangular meshes, it can also be described in terms of free-form parametric surfaces such as NURBS surfaces [7]. Moreover, due to the advantages of high precision in modeling with low number of patches and great flexibility in simulation algorithms, NURBS has been introduced into the area of high-frequency electromagnetic analysis [8, 9]. When object is described by NURBS surfaces, geodesic path can be obtained by solving the geodesic differential equation [10]. Hence, some numerical schemes were employed to handle the GDE. In [11], to obtain the tracks of creeping rays on NURBS surfaces, the geodesic differential equations are solved by Euler method which is the first-order scheme for differential equations. And in [12], the Runge-Kutta method of order 4 was applied to solving the GDE for determining the tracing on NURBS surfaces. Euler method is highly efficient but with low precision. On the contrary, the Runge-Kutta method is of high precision but time-consuming in the process of creeping ray tracing.
In this paper, aiming at improving the accuracy and efficiency of creeping ray tracing on arbitrarily shaped NURBS surface, we present a creeping ray tracing algorithm based on a novel adaptive variable step Euler method. Since the adaptive variable step Euler method is based on the conventional Euler method, the efficiency is ensured. And in the process of numerically iteratively solving the GDE, due to the introduction of shape factor , the discrete step sizes are adaptively corrected timely. Hence, compared with the conventional Euler method, the proposed method can easily ensure the accuracy of creeping ray tracing on arbitrarily shaped NURBS surface. Namely, it is more suitable for engineering applications. The efficiency and accuracy of the proposed method will be analyzed and verified by comparing the results of different methods.
2. NURBS Modeling for Arbitrarily Shaped Convex Objects
NURBS are generations of nonrational B-splines and rational and nonrational Bezier curves and surfaces [10]. And NURBS surface is the rational generalization of the nonrational B-spline surface. It is defined as follows:where are the weights and are the control points; are the normalized B-spline basis functions of degree defined recursively aswhere are the so called knots forming a knot vector . The degree, number of knots, and number of control points are related by the formula .
In many applications of computational geometry, NURBS surfaces are often expanded in terms of rational Bezier patches by applying the Cox-De Boor transform algorithm [13] since Bezier patches are numerically more stable than NURBS surfaces, while NURBS are more efficient for the storage and representation of a model. The surface points of a rational Bezier surface are given bywhere are the control points and are the weights. and are Bernstein polynomials of degrees and , respectively, which are defined for all integers , by
A NURBS or Bezier surface can be defined as a vector-valued mapping from two-dimensional parameter domain to a set of three-dimensional coordinates. And they can be generally described asThe mapping of -parameter domain to three-space can be seen in Figure 1.
In Figure 1(b), when we fix and let vary, then depends on one parameter and is, therefore, a curve called a -parameter curve. Similarly, if we fix , then the curve is a -parameter curve. Note that both curves pass through in . Tangent vectors for the -parameter and -parameter curves are given by differentiating the components of with respect to and , respectively. We write
In this paper, if and are both constant on parametric surfaces, we call the surface uniform grid surface, otherwise nonuniform grid surface. Here some typical NURBS models are given in Figure 2. The surfaces of cylinder are uniform grid surfaces, and others are nonuniform grid surfaces.
(a) Cylinder
(b) Cone
(c) Sphere
The surfaces of models above are digitized by networks of -isoparametric curves. And the cylinder, cone, and sphere are, respectively, constructed by 4 Bezier patches highlighted in different colors. And in the following, for convenience, different Bezier patches will be shown in the same color.
3. Algorithm of Creeping Ray Tracing on Arbitrary NURBS Surfaces
According to the positions of the source and the observation points, the diffraction problems of convex surfaces can be classified as three types. First, the source and the observation points are both out of the surfaces and the observation points are far away from the surfaces. This case belongs to the scattering problem of the smooth convex surfaces. Second, the source is directly on the surfaces, and the observation points are far away from the surfaces. This can be called the radiation problem of antennas. Third, the source and the observation points are both on the surfaces. This is the mutual coupling problem among antennas on the surface. Hence, the ray tracing should also be divided into three cases, as in Figures 3–5. In all of these cases, the ray trajectories on the surfaces are called creeping rays, which are restrained to propagate along geodesic paths. And this paper mainly focuses on the tracing of creeping rays since it is the most difficult part in the process of tracing.
(a) Side view
(b) Top view
(a) Side view
(b) Top view
(a) Side view
(b) Top view
3.1. The Start Points and Exit Points of the Creeping Rays
As shown in Figure 6, at the starting point , the grazing incident condition satisfieswhere is the unknown on the surface and is the corresponding unit normal vector of at .
At the exiting point ,
According to (7) and (8), we can obtain a number of the effective starting and exiting points of creeping rays on NURBS surfaces.
3.2. Creeping Ray Tracing by Solving GDE with Adaptive Variable Step Euler Method
As we know, the paths of creeping rays on arbitrarily shaped surfaces satisfy GDE. Then the problem of obtaining creeping rays is changed into solving GDE. Generally, GDE are solved by Euler method which is a simple and fast method. However, according to study, in most of the cases conventional Euler method does not work because of its low accuracy and stability.
For the nonlinear problems, the step controlled correction procedure is essentially needed. Hence, in this paper, in order to improve the accuracy of ray tracing and make sure of the efficiency, an adaptive variable step Euler method is presented to solve the GDE.
Firstly, GDE are given as follows:where are the Christofel symbols of the second kind, is arc length along the geodesic path, and are parametric coordinates of creeping rays.
Then, the first and second derivatives of in (9) can be approximately expressed bywhere is the step size between two adjacent discrete points, the determination of which is of great importance. is shape factor (SF), which is introduced to adaptively control every discrete step size. And the value of is subject to the shape of the object; more details about will be presented in Section 3.3.
In numerical calculation, a number of discrete points can be computed to express the creeping ray, where . So, according to discretization, (11) and (12) can be rewritten aswhere , which is a test point. Sincethen
Next, substituting (16) into (14) we obtain
Similarly, the first and second derivatives of can be
Finally, substituting (13), (17), and (18) into the GDE (9) and (10), we havewhere
According to (19), to calculate on the creeping path, the previous two points , and are needed. Therefore the first two points , and should be calculated at the beginning of creeping ray tracing.
Here, the starting point can be obtained by (7). And according to the differential geometry [14], the second point is approximated by the sum of the two tangent vectors on tangent plane at the starting point shown in Figure 7.
Hence, the second point can be expressed bywhere is initial step size and and are the unit tangent vectors for the -parameter and -parameter curves at . And tangent vectors and for the -parameter and -parameter curves can be obtained by (6). is the unit vector along the incident direction.
The general expression of shape factor is derived in Section 3.3. According to the expression, can be obtained.
After the calculation of , and , (19) can iterate step by step with the increase of .
3.3. The Derivation of Shape Factor
As we know, the smaller the step size is, the more accurate that derivative approximate expression is. However, in the process of iteratively solving the discrete points on creeping ray path, the efficiency of algorithm will decrease with the increase of the number of discrete points. More importantly, the more the discrete points, the larger the accumulative error which may lead to wrong results.
Clearly, for the approximation of or , we could achieve greater efficiency if we could provide more discrete points where the variables are changing rapidly to achieve accuracy, while less discrete points should be provided in regions where variables are changing slowly. Therefore it is necessary to allocate discrete points during the process of Euler approximation.
However the fact is that the parametric expressions of the creeping ray path are unknown. We cannot directly obtain the variable change speed degree which the policy of allocating discrete points is based on. So, the shape factor is introduced to determine the discrete step size. The discrete value of shape factor is defined to reflect the relative change rate of parameters at two adjacent points on the creeping path. Since the creeping rays include two parametric expressions and , we should consider the change of , simultaneously for well determining . Hence,where , , , and are unknown, which need to be solved out. And the relationship between the SF and step size is
Let be an arc length parametric curve (creeping ray) on surface which pass through point as shown in Figure 8 and denote thatLet be a unit tangent vector of at , and let the number of the discrete points be .
Tangent vector for the creeping ray at the point can be obtained by differentiating equation (24) with respect to the arc lengthwhere and can be calculated by (6). Since ,However, and cannot be obtained from (26).
Now, we approximate and at point on by at point on -parameter curve () and at point on -parameter curve, respectively. For on -parameter curve, , (26) can beFor on -parameter curve, , (26) can beSo, according to (27) and (28),Similarly, at the point on ,Substituting (29) and (30) into (22), finally the shape factor expression isSo, based on the , we can obtain discrete variable step adaptively.
4. Results and Discussions
4.1. Creeping Ray Tracing Algorithm Validation
The creeping rays can be calculated theoretically on some typical objects, such as cylinder, cone, and sphere. Therefore, the proposed method can be verified from the analytical results of these canonical objects.
(1) Creeping Ray Tracing on Cylinder (Uniform Grid Surfaces). The cylinder which has a radius of 1 m and height of 3 m is considered for creeping ray tracing. We consider three different incident wave directions, all of which incident start from the same point . In order to validate the accuracy of the proposed creeping ray tracing method, it is compared against the analytical method and the conventional Euler method. Figure 9 shows the creeping ray tracing results obtained from the above mentioned methods, for three different incident directions, respectively. It is shown that, for each incident direction, the creeping ray obtained from the proposed method overlays with those from the other two methods. To further validate the proposed method, the creeping ray ending point and the length of the creeping rays are compared, as shown in Table 1.
From Table 1, the numerical results by the proposed method are in good agreement with the analytical method, which can confirm its rationality. And we can find that the numerical results by the proposed method and Euler method, respectively, are identical, and this is because cylinder is so special that the and of creeping ray on its surfaces are first-order linear. Namely, the discrete step sizes can be chosen equally. And in this special case (uniform grid surfaces), , the proposed method is reduced to conventional Euler method. (Comparing to the Euler method, the advantages of the proposed method can be demonstrated in the following examples.)
(2) Creeping Ray Tracing on Cone (Nonuniform Grid Surfaces). The cone which has a radius of 1 m and height of 2 m is considered for creeping ray tracing. We consider three different incident wave directions, all of which incident start from the same point . In Figure 10, the creeping ray tracing results on a cone are shown. The numbers 1, 2, and 3 in Figure 10(b) are three creeping rays of different incident directions.
(a) Side view
(b) Top view
Firstly, from Table 2, we can find the numerical tracing results by proposed method are in good agreement with the analytical method. As a comparison, the results obtained by the Euler method are in lower accuracy. Moreover, the creeping ray results obtained by the Euler method deviate from that of the analytical method when the creeping ray trajectory is closer to cone-tip. This is clearly shown in Figure 10 and Table 2, the Euler method results deviate from the analytical results for the ray number 2, and it even falls into false results when the ray runs too close to the cone-tip as for the ray number 3. This indicates that the conventional Euler method cannot be applied to nonuniform grids surfaces while the proposed method is capable of handling the nonuniform grids surfaces. Secondly, the proposed method is as efficient as conventional Euler method.
(3) Creeping Ray Tracing on Sphere (Nonuniform Grid Surfaces). The sphere which has a radius of 1 m is considered for creeping ray tracing. According to differential geometry, the geodesic paths on sphere are orthodromes. Hence, the numerical results of sphere can be validated by the theoretical results.
From Figure 11(a), the result by proposed method is in good agreement with the theoretical results. And from Table 3 and Figures 11(b) and 11(c), although, with the decrease of discrete step size, the error by Euler method is reduced, time consumption is increased a lot. More importantly, if the discrete step size continues to decrease, from Figure 11(d), the result is totally wrong, which means Euler method is unstable.
(a) Proposed method
(b) Euler method
(c) Euler method
(d) Euler method
4.2. Electromagnetic Calculation with UTD on the Base of Creeping Ray Tracing Algorithm
To assess the accuracy of UTD diffraction based on the proposed creeping ray tracing algorithm, the analytical results of canonical targets are given. Figures 13 and 15 show the bistatic scattering electric fields of PEC sphere and cylinder in shadow region, respectively, and the UTD solutions are compared with the analytical results.
Moreover, in reality, lots of targets are low detectable. When analyzing the scattering properties of those targets especially for some constructed by smooth convex curved surfaces, the contributions of creeping waves cannot be neglected. Here, taking an antiradar stealth screen component of a radar stealth satellite as example, since the antiradar stealth screen component is constructed by some smooth convex curved surfaces (see Figure 16), the contributions of creeping waves will play important roles in this case. In the following example, the monostatic RCS of the antiradar stealth screen are studied by comparing the PO + UTD solution with the numerical solution to demonstrate the contribution of creeping waves.
(1) The bistatic diffracted electric fields of a PEC cylinder in shadow region are calculated by UTD. The radius of the cylinder is 1 m, and the distance between the observation point and the origin is 2 m. The frequency of incident plane wave is 3 GHz. The sketch map of the diffracted contribution is given in Figure 12.
From Figure 13, the UTD results of PEC cylinder in the shadow region based on the proposed creeping ray tracing algorithm are in good agreement with the analytical results, and the good agreement confirms the validity of the creeping ray tracing algorithm.
(2) The bistatic diffracted electric fields of a PEC sphere in shadow region are calculated by UTD. The radius of the sphere is 1 m, and the distance between the observation point and the origin is 2 m. The frequency of incident plane wave is 3 GHz. The sketch map of the diffracted contribution is given in Figure 14.
As shown in Figure 15, the UTD results of PEC sphere in the shadow region based on the proposed creeping ray tracing algorithm are in good agreement with the analytical results, which denotes the validity of the creeping ray tracing algorithm.
(3) Monostatic RCS of the antiradar stealth screen component of a radar stealth satellite are calculated by the PO method combined with the UTD method based on the proposed creeping ray tracing algorithm, and these results are then compared with that of the numerical method to validate the accuracy of the proposed method. The side and 3D view of the antiradar stealth screen component are shown in Figure 16, where m, m, and m.
Figure 16 shows the trajectories of creeping rays on the antiradar stealth screen component. In space, the conical tip of the antiradar stealth screen is directed to the earth. Hence we focus on the analysis of scattering property about the antiradar stealth screen nearby region of the conical tip. Here, the RCS will be calculated in the range of . The frequency of incident plane wave is 3 GHz.
Next, we compare the monostatic RCS of the target in the range of by three different methods, namely, the numerical method Multilevel Fast Multipole Algorithm (MLFMA), the PO method, and the PO + UTD method. In the PO + UTD scheme, the UTD solution based on the proposed algorithm is used to obtain the diffracted contribution of creeping rays. In the following, the monostatic RCS of the antiradar stealth screen are given in Figure 17.
From Figure 17, it is observed that the contribution of creeping waves is very obvious when the incident direction . In this range, thanks to the contributions of the creeping waves, the PO + UTD results are about 5 dBsm higher than the PO results. So, the contributions of creeping waves cannot be neglected. And compared with PO, the PO + UTD solution has a better agreement with the MLFMA solution.
Hence, the rationality and validity of the creeping ray tracing algorithm presented in this paper are proved.
5. Conclusion
An accurate and efficient creeping ray tracing algorithm based on an adaptive variable step Euler method is presented in this paper. The proposed algorithm is applicable to arbitrary NURBS surface. This adaptive variable step Euler method, based on the conventional Euler method, employs a shape factor (SF) to solve the GDE. The SF can timely reflect the relative change rate of the creeping ray parameters at two discretely adjacent points on surfaces. And according to , every adaptive step size, which is very important for the accuracy of creeping ray tracing, can be obtained in turn. Finally based on the fast and accurate creeping ray tracing, the contribution of UTD diffraction can be calculated. Bearing in mind the accuracy and efficiency, the algorithm of NURBS-UTD will be very useful in practical engineering.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This work was supported by the National Natural Science Foundation of China (Grant no. 60671040, Grant no. 61001059, and Grant no. 61301061), the China Postdoctoral Science Foundation, and the Fundamental Research Funds for the Central Universities.