$30
EE6506-Assignment 3 Solved
The aim of this exercise is to compute the radar cross-section of an object of arbitrary shape. This would help us to understand the assumptions involved and design choices made in real world scenarios where the radar cross-section is computed for different aircraft designs.
In this exercise, we try to find the radar cross-section (or RCS) of a regular pentagon for an incident wave at a particular frequency and angle of incidence. Since the object can be represented in 2D as the properties of the object are invariant in the Z-direction, we will be using Surface Integral Methods to calculate the scattered fields from the object.
The incident field is taken to be a S-band radar wave at 3GHz, incident at an angle of 45◦ from the x-axis. The analysis is done in 2 different scenarios:
• Assuming the object to be made of a Perfect Electrical Conductor (PEC)
• Assuming the object to be made of Carbon Fiber
We calculate the 2D RCS of the object taking the incident wave to be of Ez polarisation. For this, we will be using Pulse basis, Delta testing approach to solve the SI equations. We define the discretisation length or the width of the basis function to be λ/15 where λ is the wavelength of the incident wave.
2 Modelling the object
We first model the pentagonal object in the 2D cartesian plane and then discretise it to get a set of points along the edges of the object. The size of the object is taken to 3λ on each side.
The 2D model is generated assuming that one vertex of the object is always along the + Y-axis. The following code generates the starting points of the pulse basis functions and the evaluation points of the delta testing function. Here ‘theta’ is the direction of the ‘first’ vertex.
%generates a regular polygon in 2D, oriented along theta, centred at origin ang = 360/n_sides; c = a/(2*sind(ang/2)); %length of line segment connecting centre to vertex vrtx_A = [c*cosd(theta); c*sind(theta)]; %vertex pointing along theta dir_AB = [-cosd(90 - (ang/2) - theta); sind(90 - (ang/2) - theta)]; rotat_ang = [cosd(ang) -sind(ang); sind(ang) cosd(ang)]; n_e = round(a/da); %total no. of discretisation points
test_pt = zeros(2, n_sides*n_e); %points of delta testing strt_pt = zeros(2, n_sides*n_e); %starting points of pulse basis test_pt(:,1:n_e) = vrtx_A + [da/2:da:a].*dir_AB; strt_pt(:,1:n_e) = vrtx_A + [0:da:a-da].*dir_AB; for i = 1+n_e : n_e:n_sides*n_e test_pt(:,i:(i+n_e-1)) = rotat_ang*test_pt(:,i-n_e:i-1); strt_pt(:,i:(i+n_e-1)) = rotat_ang*strt_pt(:,i-n_e:i-1);
end
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
Figure 1: Modelling the Pentagonal Object
The function is written such that it will generate the discretisation points of any regular polygon. Since we analyse a pentagonal object, nsides=5 here.
3 Framing the Surface Integral Equations
In the 2D problem setup, we are analysing the scattered fields for TM Polarisation (Ez,Hx,Hy are the only non-zero components of the incident field). As the object is Z-invariant, we define a new variable φ = Ez. As the incident and scattered fields are assumed to be plane waves, if we know the transverse Electric field component, all the other components can be found out. We define φ1 as Ez outside the object and φ2 as Ez inside the object. So, φ1 and φ2 will obey the Helmholtz Equations:
= 0 (1)
where Q represents the contribution from the induced currents inside the object. We define the 2D Green’s functions as the solutions of the equation:
) (2) (3a) (3b)
On algebraically manipulating the above equations, we get:
4 Formulating the MOM Matrix
We first use the Extinction theorem to solve for the fields on the boundary. We try to solve the following equations:
I
φi(r0) =
[g1(r,r0)∇φ1(r) − φ1(r)∇g1(r,r0)] · ndlˆ
(4a)
I 0 =
[g2(r,r0)∇φ2(r) − φ2(r)∇g2(r,r0)] · ndlˆ
(4b)
On applying Boundary Conditions of tangential field continuity at the interface, we get the conditions:
φ1(r) = φ2(r); ∇φ1 · nˆ = ∇φ2 · nˆ (5) ) (6)
Using the pulse basis function, we discretise the two unknown variables, φ(r) and ∇φ · nˆ as
where pn(r) is the n’th pulse basis function that is non-zero only in the n’th interval. So, we have 2N unknown variables: and . We evaluate the equations 4a and 4b at the mid points of the N segments on the boundary (delta testing) and use the boundary integral method to frame 2N equations (N each in 4a and 4b).
So, we need to frame a 2N×2N matrix to solve for the 2N unknown variables. The corresponding n’th equations of 4a and 4b are:
) (7a) = 0 (7b)
where In is the n’th interval between the n’th and (n+1)’th discretisation points. Therefore, when we represent this system of equations in Ax = B form, A has dimensions 2N×2N, x has dimensions 2N×1, B has dimensions 2N×1.
Finding elements of A, B and x:
Elements of A
We can divide the matrix A into 4 quadrants and have indices i,j each ranging from 1 to N such that in equations 7a and 7b, • The 1st term in 7a will have its coefficients in 1st quadrant. So, A[i,j] = RIj g1(r,ri0) • The 2nd term in 7a will have its coefficients in 2nd quadrant. So, A[i,j + N] = −RIj ∇g1(r,ri0) • The 1st term in 7b will have its coefficients in 3rd quadrant. So, A[i + N,j] = RIj g2(r,ri0) • The 2nd term in 7a will have its coefficients in 4th quadrant.
So, A[i + N,j + N] = −RIj ∇g2(r,ri0)
Note that for the segment corresponding to i = j, there will be a singularity in g and ∇g (8)
where ~r = r~i0. We have ), and for x << 1
where γ is the euler constant. The singularity due to ln(x) is an integrable singularity and does not pose a problem. So, A[i,i] and A[i + n,i] are calculated in the same way as per the above conditions. Whereas, in the second term, we have ), and for x << 1
(9)
and the singularity due to 1/x gives a residue. To solve this problem we deform the boundaries into semicircles near the testing points. The sides on which the semicircles lie are governed by condition of the extinction theorem. Integrating on these semicircles gives,
Z
− dl∇g1 = 0.5
S1
Z
and, − dl∇g2 = −0.5 (10)
S2
Also note that the semicircles S1 and S2 are assumed to be infinitesimally small in radius. So, we need to numerically integrate ∇g1 and ∇g2 in the interval excluding the singularity. We split the interval length into 3 parts: 0 to (0.5 − 10−5)Dr, semicircle and (0.5 + 10−5)Dr to Dr. So,
Z Z
A[i,i + n] = 0.5 − ∇g1(r,ri0)dr −
Ileft Iright
∇g1(r,ri0)dr
(11a)
Z Z
A[i + n,i + n] =
−0.5 − ∇g1(r,ri0)dr −
∇g1(r,ri0)dr
(11b)
Ileft Iright
Elements of B
The matrix B will have the constant terms. Since the upper half of matrix A corresponds to equation 7a and lower half to equation 7b, the top N terms in B will have the incident field (RHS of 7a) and bottom N terms will be 0 (RHS of 7b). So, for i ranging from 1 to N, B[i] = φinc(i)
Elements of x
Aligning with the elements in A, the left half of A has the variable an and right half of A has the variable bn. So, top N elements in x will be ∇φ1 · nˆ and bottom N elements will be φ
So, for i ranging from 1 to N, x[i] = ai and x[i + N] = bi
The following code frames the matrices and finds the unknown matrix x.
A = zeros(2*n, 2*n); b = zeros(2*n,1); for i = 1:n for j = 1:n
Dr = strt_pt(:,j+1) - strt_pt(:,j); t_hat = Dr/norm(Dr); n_hat = [-t_hat(2), t_hat(1)];
if i==j %considering cases of singularity separately
A(i, j) = integral(@(d)green2d(k0, test_pt(:,i), strt_pt(:,j), Dr , d),0.0,1,’AbsTol’,tolabs,’RelTol’,tolrel); %The g1 term ( integrable singularity)
1
2
3
4
5
6
7
8
9
A(i+n, j)= integral(@(d)green2d(sqrt(eps_r)*k0, test_pt(:,i),
strt_pt(:,j), Dr, d),0.0,1,’AbsTol’,tolabs,’RelTol’,tolrel); %The g2 term (integrable singularity)
A(i,j+n) = 0.5 -1*integral(@(d)gradgreen2d_dot_n(k0, test_pt(:,i),
strt_pt(:,j), Dr, d, n_hat),0.0,0.5-a_sc,’AbsTol’,tolabs,’RelTol’,tolrel) ...
-1*integral(@(d)gradgreen2d_dot_n(k0, test_pt(:,i), strt_pt(:,
j), Dr, d, n_hat),0.5+a_sc,1,’AbsTol’,tolabs,’RelTol’,tolrel); %The grad(g1).n term
A(i+n,j+n) = -0.5 -1*integral(@(d)gradgreen2d_dot_n(sqrt(eps_r)*k0 , test_pt(:,i), strt_pt(:,j), Dr, d, n_hat),0.0,0.5-a_sc,’AbsTol’,tolabs,’ RelTol’,tolrel)...
-1*integral(@(d)gradgreen2d_dot_n(sqrt(eps_r)*k0, test_pt(:,i) , strt_pt(:,j), Dr, d, n_hat),0.5+a_sc,1,’AbsTol’,tolabs,’RelTol’,tolrel);
%The grad(g2).n term else
A(i, j) = integral(@(d)green2d(k0, test_pt(:,i), strt_pt(:,j),
Dr, d),0.0,1,’AbsTol’,tolabs,’RelTol’,tolrel); %The g1 term
A(i+n, j) = integral(@(d)green2d(sqrt(eps_r)*k0, test_pt(:,i),
strt_pt(:,j), Dr, d),0.0,1,’AbsTol’,tolabs,’RelTol’,tolrel); %The g2 term
A(i, j+n) = -1*integral(@(d)gradgreen2d_dot_n(k0, test_pt(:,i),
strt_pt(:,j), Dr, d, n_hat),0.0,1,’AbsTol’,tolabs,’RelTol’,tolrel);
%The grad(g1).n term
A(i+n, j+n) = -1*integral(@(d)gradgreen2d_dot_n(sqrt(eps_r)*k0,
test_pt(:,i), strt_pt(:,j), Dr, d, n_hat),0.0,1,’AbsTol’,tolabs,’RelTol’, tolrel); %The grad(g2).n term
end
end b(i) = inc_field(test_pt(:,i), theta_i, k0);
end
A = da*A; %scaling by the length of the segment since the integral iterated over d from 0 to 1 x = A\b; %solve by gaussian elimination
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
where, test pt are the testing centres and strt pt are the starting points of the pulse basis functions.
5 Computing RCS
As mentioned in the introduction, we analyse the RCS of the object for 2 cases, assuming it is made of a PEC and assuming it is made of Carbon fiber.
When we consider a PEC object, the epsilon = inf and the boundary conditions get simplified as the tangential electric field in a PEC has to be 0. So, we separately solve for this case. When we consider it to be a PEC, we need to find the relative permittivity of the material. Assuming it is made of double-carbon microfoils, we get its electric permittivity as 12 - 5.5j at a frequency of 3GHz. This value is obtained from the following plot by Tengfei Chen et. al.
From the plot, the real part of epsilon is taken as 12 and imaginary part as 5.5 at 3GHz. Since they follow a +jωt convention for forward propagating waves, the complex permittivity is obtained as 12-5.5j.
The 2D Radar Cross Section (RCS) is given by,
(12)
So, we need to find the total fields due to scattering in the far field region. Once we know the fields at the boundary, which are calculated using the extinction theorem, we can calculate the
Figure 2: Finding eps of CF Source : Reference 1
far fields using Huygen’s principle according to Equations 3a and 3b. So,
I
[g1(r,r0)∇φ1(r) − φ1(r)∇g1(r,r0)] · ndlˆ = φscat(r0) (13)
Once we get φ2(r0), we can use that in the RCS expression. Also note that we assume the incident field magnitude to be 1.
The following code uses the information about the field at the boundary to calculate the RCS.
for i = 1:n_ff for j = 1:n
Dr = strt_pt(:,j+1) - strt_pt(:,j); t_hat = Dr/norm(Dr); n_hat = [-t_hat(2), t_hat(1)]; integral_g1 = da*integral(@(d)green2d(k0, ff_pt(:,i), strt_pt(:,j), Dr,
d),0.0,1,’AbsTol’,tolabs,’RelTol’,tolrel); integral_grad_g1_dot_n = da*integral(@(d)gradgreen2d_dot_n(k0, ff_pt(:,
i), strt_pt(:,j), Dr, d, n_hat),0.0,1,’AbsTol’,tolabs,’RelTol’,tolrel); scat_field_ff(i) = scat_field_ff(i) - (fields_bndry(j)*integral_g1 -
fields_bndry(j+n)*integral_grad_g1_dot_n);
end
end
RCS = 2*pi*a_ff*abs(scat_field_ff).^2; % (the magnitude of the incident field is just 1)
1
2
3
4
5
6
7
8
9
10
11
12
We assume that the object is heading towards the radar, hence it is oriented as in Fig. 3.
Figure 3: Object pointed towards the source
The RCS plots are as shown in Fig 4 and 5. The field is incident at an angle of 45◦, i.e., it is coming from the bottom left corner. We can see that the plots are symmetric about the y = x line, as was to be expected since the object is symmetric about this line and the incident field is along this line. We can see that PEC does not scatter much in the direction of the source, which was to be expected from the geometry of a pentagon. To our surprise, carbon fiber scatters a lot in the direction of the radar.
Figure 4: RCS for PEC with object oriented as in Fig. 3
Figure 5: RCS for carbon fiber with object oriented as in Fig. 3
Figure 6: Object pointed away from the source
We also calculated for the case in which the object is oriented as in Fig. 6. For this orientation, the RCS plots are in Fig. 7 and 8.Here again we can see that PEC scatters strongly in the direction of the source, since the incident field encounters a normal incidence on the edge.
Figure 7: RCS for PEC with object oriented as in Fig. 6
Figure 8: RCS for carbon fiber with object oriented as in Fig. 6
