Skip to article content
Report

Project Elara 2024 research and development report

An end-of-year overview of research conducted by Project Elara

Appendix

Appendix A

Here, we show the derivation of the weak form of the vector-valued Helmholtz equation for our boundary conditions. As a starting point, recall that the general variational form of the vector Helmholtz equation is given by:

Ωk2ΦEdVΩJE:JΦdV+ΦJEdA=0\int_{\Omega} k^2 \Phi \cdot \mathbf{E} \hspace{0.17em} dV - \int_{\Omega} \nabla_J \mathbf{E} : \nabla_J \Phi \hspace{0.17em} dV + \int_{\partial } \Phi \cdot \nabla_J \mathbf{E} \hspace{0.17em} d \mathbf{A} = 0

Where : denotes the double dot-product. The boundary conditions we applied within our model were given by:

EnΩP=EnΩW=0,ΩR3limrE=ϵ\begin{align*} \frac{\partial \mathbf{E}}{\partial n} |_{\partial \Omega_P} = \frac{\partial \mathbf{E}}{\partial n} |_{\partial \Omega_W} = 0, & \Omega \in \mathbb{R}^3 & \\ \lim_{r \rightarrow \infty} \mathbf{E} = \epsilon & \end{align*}

To simplify modelling of far-field radiation, we used a coordinate transformation as shown in Frei, 2015, where we applied a coordinate transform xαtanhx,yβtanhyx \to \alpha \tanh x, y \to \beta \tanh y to mathematically extend the domain to infinity without needing to construct a larger grid, where, α and β are respectively the half-length and half-width of the domain bounding box. Physically, such a coordinate transform approximately replicates an open boundary that allows outgoing electromagnetic waves to radiate outwards in all directions, allowing the radiative boundary conditions to be imposed without needing to discretize an infinite domain.

However, the direct application of the coordinate transformation results in significant mathematical complexities due to the dependence of the representation of the variational form on the choice of coordinates. Thus tensors will be used instead to give expressions of physical laws that are invariant under coordinate transformations. From this point on, all expressions of the weak form will be given in tensors unless otherwise specified. All tensors are assumed to be within Euclidean space where upper and lower indices are equivalent, that is, Ai=AiA^i = A_i. The Einstein summation convention is assumed, in which repeated indices are implicitly summed over, and all indices take the numerical values of i=1,2,3i = 1, 2, 3 unless otherwise specified.

To begin, the weak form may be expressed in tensor notation as:

Ωk2ΦiEidVΩjEijΦidV+ΩΦijEidAj=0\int_{\Omega} k^2 \Phi_i E^i \hspace{0.17em} dV - \int_{\Omega} \partial_j E^i \partial^j \Phi_i \hspace{0.17em} dV + \int_{\partial \Omega} \Phi_i \partial_j E^i \hspace{0.17em} dA^j = 0

One may reduce by one dimension to just consider the 2D case, as was done to simplify the problem for this initial stage of research:

Ωk2ΦiEidAΩjEijΦidA+ΩΦijEidxj=0\int_{\Omega} k^2 \Phi_i E^i \hspace{0.17em} dA - \int_{\Omega} \partial_j E^i \partial^j \Phi_i \hspace{0.17em} dA + \int_{\partial \Omega} \Phi_i \partial_j E^i \hspace{0.17em} dx^j = 0

The old coordinates are denoted xi=x=(x,y)x^i = \mathbf{x} = (x, y), and new coordinates denoted xk=x=(x,y)x^k = \mathbf{x}' = (x', y') where x=x(x)x' = x' (x) and y=y(y)y' = y' (y). A change of variables was applied on the weak form xixkx^i \to x^k. On the first integral term, the Kronecker delta was used to relabel indices from iki \to k, by the relationships Φi=δikΦk,Ei=δkiEk\Phi_i = \delta_i^k \Phi_k, \quad E^i = \delta^i_k E^k, from which one can substitute into the first integral term to obtain:

Ωk2ΦiEidA=k2ΩδikΦkδkiEkdA\int_{\Omega} k^2 \Phi_i E^i \hspace{0.17em} dA = k^2 \int_{\Omega} \delta_i^k \Phi_k \hspace{0.17em} \delta^i_k E^k \hspace{0.17em} dA

The two Kronecker deltas are contracted over both of their indices, so δikδki=δii=δ11+δ22+δ33=3\delta_i^k \hspace{0.17em} \delta^i_k = \delta^i_i = \delta^1_1 + \delta^2_2 + \delta^3_3 = 3. Therefore the first term further simplifies to 3k2ΩΦkEkdA\displaystyle 3 k^2 \int_{\Omega} \Phi_k \hspace{0.17em} E^k \hspace{0.17em} dA. For the second term, the chain rule j=xkxjxk=jxkk\partial_j = \dfrac{\partial x^k}{\partial x^j} \dfrac{\partial}{\partial x^k} = \partial_j x^k \hspace{0.17em} \partial_k may be used to rewrite derivatives in terms of the new coordinates xkx^k. After substitution one finds:

ΩjEijΦidA=Ω(jxk)kEi(jxk)kΦidA\int_{\Omega} \partial_j E^i \partial^j \Phi_i \hspace{0.17em} dA = \int_{\Omega} (\partial_j x^k) \partial_k E^i (\partial^j x_k) \partial^k \Phi_i \hspace{0.17em} dA

Where, after rearrangement and using the Kronecker delta for relabeling the EE and Φ indices (as shown previously with the first term), becomes:

Ω(jxk)(jxk)kδkiEkkδikΦkdA\int_{\Omega} (\partial_j x^k) (\partial^j x_k) \partial_k \delta^i_k E^k \partial^k \delta_i^k \Phi_k \hspace{0.17em} dA

Taking out the Kronecker deltas allows for further simplification, where δkiδik=3\delta^i_k \delta_i^k = 3 due to the double contraction previously explained:

Ωδkiδik(jxk)(jxk)kEkkΦkdA=3Ω(jxk)(jxk)kEkkΦkdA\int_{\Omega} \delta^i_k \delta_i^k (\partial_j x^k) (\partial^j x_k) \partial_k E^k \partial^k \Phi_k \hspace{0.17em} dA = 3 \int_{\Omega} (\partial_j x^k) (\partial^j x_k) \partial_k E^k \partial^k \Phi_k \hspace{0.17em} dA

The third term may be simplified via the tensor transformation law dxj=xjxkdxk=kxjdxkdx^j = \dfrac{\partial x^j}{\partial x^k} dx^k = \partial_k x^j dx^k, from which one may obtain:

ΩΦijEidxj=ΩΦijEikxjdxk\int_{\partial \Omega} \Phi_i \partial_j E^i \hspace{0.17em} dx^j = \int_{\partial \Omega} \Phi_i \partial_j E^i \partial_k x^j dx^k

Now substituting in j=jxkk\partial_j = \partial_j x^k \hspace{0.17em} \partial_k and the two Kronecker delta identities Φi=δikΦk\Phi_i = \delta_i^k \Phi_k and Ei=δkiEkE^i = \delta^i_k E^k the result is:

ΩδikΦkjxkkδkiEkkxjdxk\int_{\partial \Omega} \delta_i^k \Phi_k \partial_j x^k \hspace{0.17em} \partial_k \delta^i_k E^k \hspace{0.17em} \partial_k x^j dx^k

Where, after another double contraction over the Kronecker deltas and rearranging, noting that jxkkxj=1\partial_j x^k \hspace{0.17em} \partial_k x^j = 1, gives:

3ΩΦkjxkkxjkEkdxk=3ΩΦkkEkdxk3 \int_{\partial \Omega} \Phi_k \partial_j x^k \hspace{0.17em} \partial_k x^j \hspace{0.17em} \partial_k E^k \hspace{0.17em} dx^k = 3 \int_{\partial \Omega} \Phi_k \partial_k E^k \hspace{0.17em} dx^k

Altogether, after substitution of all simplified terms and ordering the terms such that the quasi-linear term is second due to software requirements, one finally obtains the weak form in the transformed coordinates:

3ΩjxkjxkkEkkΦkdA+3k2ΩΦkEkdA+3ΩΦkkEkdxk=0- 3 \int_{\Omega} \partial_j x^k \hspace{0.17em} \partial^j x_k \hspace{0.17em} \partial_k E^k \partial^k \Phi_k \hspace{0.17em} dA + 3 k^2 \int_{\Omega} \Phi_k \hspace{0.17em} E^k \hspace{0.17em} dA + 3 \int_{\partial \Omega} \Phi_k \partial_k E^k \hspace{0.17em} dx^k = 0

After which some algebraic simplifications gives the most simplified general form:

ΩjxkjxkkEkkΦkdA+k2ΩΦkEkdA+ΩΦkkEkdxk=0- \int_{\Omega} \partial_j x^k \hspace{0.17em} \partial^j x_k \hspace{0.17em} \partial_k E^k \partial^k \Phi_k \hspace{0.17em} dA + k^2 \int_{\Omega} \Phi_k \hspace{0.17em} E^k \hspace{0.17em} dA + \int_{\partial \Omega} \Phi_k \partial_k E^k \hspace{0.17em} dx^k = 0

In typical vector calculus notation rather than tensor notation, this may alternatively be written as:

Ω(Jx:Jx)(E~)(Φ)dA+k2ΩΦE~dA+ΩΦ(E~)dx=0- \int_{\Omega} (\nabla_J \mathbf{x}' : \nabla_J \mathbf{x}') (\nabla \cdot \tilde{\mathbf{E}}) (\nabla \cdot \Phi) \hspace{0.17em} dA + k^2 \int_{\Omega} \Phi \cdot \tilde{\mathbf{E}}\, d \mathbf{A} + \int_{\partial \Omega} \Phi \cdot (\nabla \cdot \tilde{\mathbf{E}}) \hspace{0.17em} d \mathbf{x} = 0

When applying boundary conditions to the boundary integral, the only contribution is that of the radiative boundary condition, which reduces to the constant Dirichlet boundary condition E~ΩB=ϵ\widetilde{\mathbf{E}} |_{\partial \Omega_B} = \epsilon and thus:

ΩΦ(E~)dx=0\int_{\partial \Omega} \Phi \cdot (\nabla \cdot \tilde{\mathbf{E}}) \hspace{0.17em} d \mathbf{x} = 0

And therefore the boundary integral vanishes. We have therefore obtained our most simplified weak form:

Ω(Jx:Jx)(E~)(Φ)dA+k2ΩΦE~dA=0- \int_{\Omega} (\nabla_J \mathbf{x}' : \nabla_J \mathbf{x}') (\nabla \cdot \tilde{\mathbf{E}}) (\nabla \cdot \Phi) \hspace{0.17em} dA + k^2 \int_{\Omega} \Phi \cdot \tilde{\mathbf{E}}\, d \mathbf{A} = 0

Appendix B

In this section, we show the domain parametrization used for our model. To start, the Laplacian in a generalized Euclidean coordinate transformation xixjx^i \to x^j takes the form:

2=iixjj+ixjixjjj\nabla^2 = \partial^i \partial_i x^j \partial_j + \partial^i x_j \partial_i x^j \partial^j \partial_j

This result is obtained using the transformation rule for partial derivatives, which is simply the chain rule:

xj=xixjxi=jxii\dfrac{\partial}{\partial x^j} = \dfrac{\partial x^i}{\partial x^j} \dfrac{\partial}{\partial x^i} = \partial_j x^i \partial_i

Which is then substituted into the definition of the Laplacian 2=kk\nabla^2 = \partial^k \partial_k and distributed to obtain the result in (17). Using this expression, with the coordinate transforms x=αtanhxx' = \alpha \tanh x, y=βtanhyy' = \beta \tanh y, the transformed Laplacian becomes:

2=α2(1x2)[x((1x2)x)]+β2(1y2)[y((1y2)y)]\nabla^2 = \alpha^2 {(1 - x'}^2) \left[ \frac{\partial}{\partial x'} \left( {(1 - x'}^2) \frac{\partial}{\partial x'} \right) \right] + \beta^2 {(1 - y'}^2) \left[ \frac{\partial}{\partial y'} \left( {(1 - y'}^2) \frac{\partial}{\partial y'} \right) \right]

Thus the Helmholtz equation in the aforementioned transformed coordinates has the form:

α2(1x2)[x((1x2)x)]+β2(1y2)[y((1y2)y)]E~(x,y)+k2E~(x,y)=0\begin{align*} \alpha^2 {(1 - x'}^2) \left[ \frac{\partial}{\partial x'} \left( {(1 - x'}^2) \frac{\partial}{\partial x'} \right) \right] &+ \beta^2 {(1 - y'}^2) \left[ \frac{\partial}{\partial y'} \left( {(1 - y'}^2) \frac{\partial}{\partial y'} \right) \right] \tilde{\mathbf{E}} (x', y') \\ &+ k^2 \tilde{\mathbf{E}} (x', y') = 0 \end{align*}

The geometry of the simulation was parametrized by four constants, the simulation volume width DD and length LL, the primary reflector radius R1R_1, secondary reflector radius R2R_2, opening gap radius bb (where the opening is located at the rear of the primary reflector dish), and the primary reflector anchor point d1d_1. The primary collector and secondary mirror were respectively parametrized as follows:

x1(t)=d112R1t2,y1(t)=t,t[R1,b][b,R1]x2(t)=d1R12+12R1t2,y2(t)=t,t[R2,R2]\begin{align*} x_1 (t) = d_1 - \frac{1}{2 R_1} t^2, \quad y_1 (t) = t, \quad &t \in [- R_1, - b] \cup [b, R_1] \\ x_2 (t) = d_1 - \frac{R_1}{2} + \frac{1}{2 R_1} t^2, \quad y_2 (t) = t, \quad &t \in [-R_2, R_2] \end{align*}

The bounding box of the simulation was parametrized as follows, for which t[0,1]t \in [0, 1]:

xT(t)=LtL2,yT(t)=D2xB(t)=LtL2,yB(t)=D2xL(t)=L2,yL(t)=DtD2xR(t)=L2,yR(t)=DtD2\begin{align*} x_T (t) &= Lt - \frac{L}{2}, & y_T (t) = \frac{D}{2}\\ x_B (t) &= Lt - \frac{L}{2}, & y_B (t) = - \frac{D}{2}\\ x_L (t) &= - \frac{L}{2}, & y_L (t) = Dt - \frac{D}{2}\\ x_R (t) &= \frac{L}{2}, & y_R (t) = Dt -\frac{D}{2} \end{align*}
References
  1. Frei, W. (2015). Using Perfectly Matched Layers and Scattering Boundary Conditions for Wave Electromagnetics Problems.