The diffraction of waves in free space is governed by the sourceless Helmholtz equation $$\nabla^2u +k^2u=0$$ There is a well stablished solution called Kirchhoff´s diffraction integral (KDI) and it´s deduction is straightforward from greens theorem, the solution is $$ u(r) = \int_{S} \left [ u \space \partial_n \left(\frac{e^{iks}}{s} \right) - \left(\frac{e^{iks}}{s} \right) \space \partial_n u\space \right] dS$$

Where $S$ is a close surface containing the aperture and an imaginary boundary at infinity that fulfills Sommerfeld radiation condition and $s$ is the distance between any point in the surface and the point $r$ inside the enclosed volume.

Now I am trying to solve the Helmholtz equation numerically using the Finite Element Method (FEM). For this, I am assuming that a scalar plane wave propagating in the z-direction is illuminating a 1D screen with an aperture at the center. To set up the model, I am assuming the following Neumann conditions for the dependent variable u(x,z):

Absorbing boundary condition (ABC) to model an open system (boundaries at infinity) and ABC for the diffraction obstacle too, $n\cdot \nabla u=-ik_0 u$

radiation condition for the aperture to model illumination by a plane wave, $n \cdot \nabla u= ik_0 ((1+n\cdot e)-u)$

where $n$ is the boundary normal and $e$ is the incoming wave propagation direction. The solution obtained with this model is similar to the one using the KDI, but there is a mismatch in the scale and phase of the solution

numerical results of the model vs the analytical Kirchhoff diffraction formula

I have tried different boundary conditions, for example, Sommerfeld's and Kirchhoff's, the ones used to deduce the KDI, but the results are similar. The solutions differ in a factor of $\sqrt{s}$ in the 2D case.

I found that the numerical and theoretical solutions match in the 3D case, i.e. u(x,y,z) and the aperture being a surface instead of a line as in the 2D case. The problem remains in how to reduce the KDI to 2D instead of 3D, normally when one consider a 2D diffraction problem, one assumes that, for example, the aperture extend to infinity in the $y$ direction, therefore it is not necessary to take the $y$ coordinate into account. I think that maybe the $\sqrt{s}$ factor comes from the reduction to 2D of the KDI. I tried doing the integration from $[-\infty,\infty]$ numerically, but it didn't work.

I will appreciate if someone knows what is wrong with the model or how to reduce de KDI to 2D rigorously.

Note: The boundary conditions were taken from Jian Ming's "The Finite Element Method in Electromagnetism".

  • 1
    $\begingroup$ A guess: in 3D, the electric potential goes as -1/s; but in 2D, it should go as Log(s). Should come from Gauss divergence theorem/Helmoltz eq , hence the Log(s) should appear in the Green's function. Put the Log(s) in the KDI instead of 1/s? Needs some trickery as elecromagnetism in 2D is divergent; your attempt to do the integral in y probably crashed against this divergence. $\endgroup$ – patta Apr 12 at 6:49
  • $\begingroup$ You may try to integrate the KDI / (2y) from -y to +y; then let y go to infinity...? $\endgroup$ – patta Apr 12 at 6:58

Short instruction. In 2D case one needs instead of $\frac{e^{ikr}}{r}$ to put a Bessel function of the second kind. The 2D version of KDI is the following: $$ u(r)=\frac{1}{4} \left[ \int_S u \partial_n Y - Y \partial_n u \right]\quad(1) $$ where $$ Y=-Y_0(k r), $$ where $Y_0(x)$ is a Bessel function of the second kind $$ \pi Y_0(x)=\int_0^\pi \sin(x \sin{\theta})d \theta - \int_0^\infty e^{-x \sinh(t)}dt. $$ For large $x$ one can use an approximation $$ Y_0(x)\approx-i\frac{1}{\sqrt{2\pi x}}e^{i(x-\frac{\pi}{4})}. $$ You can notice that your phase shift is about $\pi/4$ everywhere, so this formula is quite likely applicable to your case.

Some background. For any two functions, $u$ and $u_1$, satisfying the Helmholtz equation, the following integral can be shown to be zero $$ \int_S u \partial_n u_1 - u_1 \partial_n u = 0.\quad (2) $$ Function $u_1=Y_0(k r)$ satisfies Helmholtz equation everywhere except at $r=0$. Picking an area that includes $r=0$, and cutting out a tiny disk out around $r=0$ ensures that the integral is zero over such area. The integral over the perimeter of the tiny disk gives us $u(r=0)$ (with some coefficient), while the integral over the big outside perimeter gives us the right side of the Eq. (1).


Another good method is the Fourier approach. First decompose your aperture function in terms of a set of plane waves propagating in different directions, then propagate each plane wave to the observation plane, then add them up again. Thus two Fourier transforms with a simple phase factor introduced in between them. Sorry I don't have time now to write out the details but I hope this might be helpful anyway.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.