Document Type : Research Paper
Authors
Electrical and Electronic Engineering Department , Shahed University , Tehran , Iran
Abstract
Keywords
Controlling electromagnetic waves has broad range of applications in various area of science and engineering, for instance the design of perfect flat lenses, waveguides, invisibility cloaks or in general metamaterials. The idea of controlling electromagnetic waves by introducing materials inhomogeneity was introduced in [1] by Pendry et. al. It stems on this fact that the Maxwell equations are invariant under the coordinate transformation. Therefore, by finding an analytic mapping from the physical (undeformed) space to the clocked (deformed) space in which the transformation maps a point onto a ball with the finite reduce, the mentioned ball will be hidden from the electromagnetic waves. Later, this method has been extensively used by various researchers for the design of metamaterials [2]-[4], photonic band gap devices [5], active nanoplasmonic metamaterials [6], extreme-angle broadband metamaterial lens [7], all-dielectric, polarization-independent optical cloaks [8].
There are several efforts in the literature for the design of invisibility cloaks, in particular within the frequencies spectrum of the visible light, for instance the optical cloaking with metamaterials [9], full-wave invisibility of active devices at all frequencies [10], macroscopic invisibility cloak for visible light [11]-[13]. However, as it is already warned by Pendry et. al. [1], the need for highly anisotropic inhomogeneous media is the main difficulty of their method. This makes it difficult to realize these materials for real world applications. To cope this limitation, some authors use conformal mapping to generate nearly isotropic cloaks, for instance see [14]-[16]. The use of constrained optimization and quasi-conformal mapping [17]-[19] is another alternative approach to overcome this limitation. Computing the mapping based on the laplace equation [20]-[21] was another idea to reduce the material anisotropy.
The use of topology optimization method is an alternative way to control the propagation of electromagnetic waves in a heterogeneous media. In [22]-[25], the topology optimization approach was used for the design of directional cloaks considering a plane wave with a specific frequency. In these works the optimal design problem leads to the solution of a distributed parameter identification problem under the Maxwell equation as a constraint which was computationally very expensive. In [26] a simple method for the design of directional cloaks has been suggested based on the laplace equation. In this method, the light trajectories are computed by the solution of the laplace equation and then the material properties, the index of refraction, was computed directly using the eikonal equation. This method produces the isotropic material properties and does the perfect cloaking under the geometric optics approximation. However, there is no bound on the material properties.
In the present study, the topology optimization method is adapted to formulate the original problem as a distributed parameter identification problem. However, unlike formerly mentioned works, it is assumed that the light trajectories follow the eikonal equation. Therefore, we avoid to solve the Maxwell equation which is computationally very expensive. Instead, we solve the eikonal equation that can be solved very efficiently by the fast marching or fast sweeping methods.
The remainder of this paper is organized as follows. The problem formulation is presented in the next section. The sensitivity analysis and necessary optimality conditions are presented in section 3. The solution algorithm is presented in section 4. Section 5 is devoted to the space discretization and numerical algorithm. The numerical results for a model problem is presented in section 6. Finally, section 7 summarizes this work and presents the outline of further studies.
Consider () as an open bounded domain with sufficiently regular boundaries occupied with two isotropic materials. Assume the characteristic length scale of is denoted by. Consider arbitrary electromagnetic wave with wavelengths (considering non-monochromatic waves includes no technical difficulty). It is well known that under assumption of the wave (light) trajectories obeys the following eikonal equation:
(1)
Where denotes the first arrival time of wave at , denotes the index of refraction and denotes the set of point source in. Note that isosurfaces of are equivalent to the wave fronts, i.e. the wave trajectory at is orthogonal to -isosurface.
The goal of optimal design problem in the present study is to determine the optimal distribution of in the design domain such that the waves propagate according to a-priori desired trajectories in the control domain. For this purpose, the objective functional, to be minimized, is defined as the distance of and target wave fronts over,
(2)
It is worth mentioning that the objective function can defined in more general form to result other desires. For instance, to make hidden from the light, i.e. to make invisibility cloaks for visible light spectrum, we can maximize the arrival time of light into , and on the other hand minimize the arrival time of light in the remained areas, i.e. in . For this purpose, the objective function can be defined alternatively as follows:
(3)
To make the optimal design manufacturable, is constrained to lies in the admissible design space which is defined as follows:
(4)
Where denotes the admissible lower and upper bound for pointwise refractive index respectively. Therefor the optimal design problem in the present study can be summarized as follows:
(5)
To minimize the objective functional, the constrained gradient flow of (2) or (3) is used in the present study. To do this, it is required to compute the fist variation of objective functional (2) or (3) with respect to the design variable. Because after the space discretization, the number of design parameters will be possibly very large, the adjoint approach is adapted here to compute the first variation of (2) or (3) with respect to .
Consider as the Lagrange multiplier corresponding to PDE constraint (1). Ignoring the bound contestants (they are treated by the projected gradient approach later), the local solutions of (5) are subsets of stationary points of the following augmented lagrangian:
(6)
Taking the first variation of with respect to and equating it to zero will be identical to the direct PDE (1). The first variation of with respect to will be computed as follows (consider objective functional (2)):
(7)
Consider as a scalar field and as a vector field, applying identity to the second term in right hand side of (6) results in (take and):
(8)
Applying the divergence theorem on the first term on the right hand side of (8) results in:
(9)
Where denotes the boundaries of, i.e. and denotes the outward unit normal vector on. Considering (7) to gather with (9) and equating to zero the first variation of with respect to results in the following adjoint PDE:
(10)
Where denotes the characteristic function corresponding to, i.e. in and in. Note that the boundary conditions of adjoint PDE is computed by equating to zero the first term in the right hand side of (9). If alternatively the objective functional (3) is considered, by similar derivation we reach to the following adjoint PDE:
(11)
Finally, taking the first variation of with respect to leads to:
(12)
Therefore, equating it to zero the first variation of with respect to results in:
(13)
Treating the bound constraints by the projected gradient approach (cf. [27], [28]), the first order KKT necessary optimality conditions (cf. [2]) corresponding to the optimization problem (5) can be expressed as follows:
where operator denotes the orthogonal projection of trial fields onto design space . By straightforward computation, the explicit form of operator can be computed as follows (cf. [2]):
For the case of objective functional (3) the necessary first order optimality conditions will has the following form:
Considering results of section III, the solution for (5) based on the constrained gradient flow could be expressed as follows:
|
Algorithm 1. Constrained gradient flow of (5):
|
where denotes the characteristic function of the design space, i.e. in and in . Because the gradient is computed based on the norm, the regularity of the new design may not be better than (is a square integrable field). In fact, it can possesses discontinuity in which is not desirable from the practical and computational point of views. It can decrease the convergence rate of algorithm 1, or even leads to oscillatory behavior in the course of computations. To overcome this problem, we perform a regularization step in algorithm 1 and lift the new design variable from the function space to the more regular space in which both of function and its first order derivative are square integrable. For this purpose, the following Helmholtz equation should be solved at every step:
(14)
whereis a regularization parameter that controls the smoothness of solution and denotes the regularized design variable. According to the above discussion, algorithm 1 could be modified as follows:
|
Algorithm 1. Regularized constrained gradient flow of (5):
|
For the purpose of numerical solution, the spatial domain is assumed to be a rectangular domain with dimension unites and is discretized into a uniform squared grid with edge length. The spatial variables are defined at center of finite difference cells, and is used to denote the discretized counterpart of at . Because, direct equation (1) is a nonlinear first order partial differential equation, using conventional methods for its solution leads to a nonlinear system of equations which is computational expensive to solve in practice. However, exploiting the hyperbolic and causality nature of this equation leads to a very efficient solution algorithm called the fast marching method [28] with computational complexity; where. Later, the fast sweeping method has been introduced in [31] by Zhao
Fig. 1. Contour plot of first arrival time for iteration 10 (legend in time unit).
and it has been shown that it exhibits computational complexity to achieve a moderate (reasonable) accuracy. Due to the high computational performance and ease of implementation, the fast sweeping method will be used in the present study to solve (1). The numerical method in the
present study is closely related to method introduced in [32]. This method will be briefly discussed in this paper in the appendix section.
In this section, the performance and success of the presented algorithm is studied by a numerical experiment. For this purpose consider as a rectangular domain with edge lengths and units along and directions respectively. Moreover, assume that the spatial domain is decomposed into uniform square grids with grid size. The control domain, , is defined as a circle with radius located at the center of the spatial domain. The design domain, cloak in this example, is defined as a ring around the control domain with external radius. Our goal is to hide the control domain from light rays propagated from a point source which is located at the center of left edge of the spatial domain. For this purpose, the objective function is defined such that its minimizer leads to very large arrival time in and minimal arrival time in other regions, i.e. objective functional (3) is considered here. Figs. 1-5 show the contour plot of first arrival times after 10, 20, 30, 50 and 100 iterations.
Fig. 2. Contour plot of first arrival time for iteration 20 (legend in time unit).
Fig. 3. Contour plot of first arrival time for iteration 30 (legend in time unit).
Fig. 4. Contour plot of first arrival time for iteration 50 (legend in time unit).
Fig. 5. Contour plot of first arrival time for iteration 100 (legend in time unit).
Fig. 6. The variation of the objective function in the course of optimization.
According to plots the first arrival time is increased inside the control domain as the optimization proceeds. This observation confirms the success of the optimization algorithm to move toward the optimal solution. Fig. 6 shows the variation of objective functional as a function of iterations during the first 100 iterations. As the figure shows, the objective functional decreased monotonically by iterations. This observation confirms the convergence of the presented numerical algorithm.
A mathematical model and its corresponding computational method is presented in this study to control the light propagation in a media by means of optimal determination of material inhomogeneity. It is based on solving a distributed parameter identification problem constrained to an eikonal equation. In contrast to available alternative methods, the main benefit of the presented approach is its ability to produce designs with bounded isotropic material properties. The success of the presented method is shown by solving a model problem corresponding to the invisibility cloaking. According to our numerical results the performance of optimized design, the objective functional, is improved by an order of magnitude by the presented numerical algorithm.
In this section the numerical solution algorithm for the eikonal equation and its corresponding adjoint equation used in the present study are briefly discussed.
Applying the Godunov numerical Hamiltonian to the eikonal equation (1), for , we have:
(A.1)
where
and denotes the positive part of, i.e. . The one sided finite difference method is used at the boundary points. The fast sweeping algorithm can be summarized as follows:
Following [30], the fast sweeping method is used in the present study to solve (11). Considering two dimensional case, equation (11) can be written in the following form:
(A.2)
where , and are function of and , i.e. , , . The space discretized form of (A.2) can be expressed as follow:
(A.3)
where and are determined according to the direction of light propagation. For instance, when, . Moreover,
Following [30], we introduce the following notation:
Using this notation results in:
(A.4)
(A.4) can be written as follows:
(A.5)
(A.5) is in fact a system of linear equation that can be solved efficiently to a reasonable accuracy by the fast sweeping method.
)