# Fictitious Domain Technique for the Calculation of Time-Periodic Solutions of Scattering Problem

Authors: Ling Rao; Hongquan Chen
Journal: Mathematical Problems in Engineering
 (2011)

---

## Abstract

The fictitious domain technique is coupled to the improved time-explicit asymptotic method for calculating time-periodic solution of wave equation. Conventionally, the practical implementation of fictitious domain method relies on finite difference time discretizations schemes and finite element approximation. Our new method applies finite difference approximations in space instead of conventional finite element approximation. We use the Dirac delta function to transport the variational forms of the wave equations to the differential form and then solve it by finite difference schemes. Our method is relatively easier to code and requires fewer computational operations than conventional finite element method. The numerical experiments show that the new method performs as well as the method using conventional finite element approximation.


---

## Body

## 1. Introduction

Recently, aircraft design for military application has focused more and more attention on using stealth technologies. It is important to realize Rader stealth through reducing the intensity of scattering signals of Rader in stealth design. Theoretically, the stealth characteristics such as Radar Cross-Section (RCS) for a given aerodynamic body can be obtained by solving the fundamental electromagnetic Maxwell equations. The control method based on exact controllability has been successfully used in computing the time-periodic solutions of scattered fields by multibody reflectors (see [1–5]). An improved time-explicit asymptotic method is afforded through introducing an auxiliary parameter for solving the exact controllability problem of scattering waves [4].Fictitious domain methods are efficient methods for the solutions of viscous flow problems with moving boundaries [6]. In [7–9], fictitious domain method is combined with controllability method to compute time-periodic solution of wave equation, which is proved to be equivalent to the Maxwell equation in two dimensions for the TM mode. A motivation for using fictitious domain method is that it allows the propagation to be simulated on an obstacle free computational region with uniform meshes. In our paper, the fictitious domain technique is coupled to the improved time-explicit asymptotic method for calculating time-periodic solutions of wave equation. Conventionally, the practical implementation of fictitious domain method relies on finite difference time discretizations schemes and finite element approximation. Our new method applies finite difference approximations in space instead of conventional finite element approximation (see [7–9]). We use the Dirac delta function to transport the variational form of the wave equation to the differential form and then solve it by finite difference schemes. Our method is relatively easier to code and requires fewer computational operations than conventional finite element method does.In Section2, the formulation of the Scattering problem is presented. In Section 3, we introduce exact controllability problem of the Scattering problem and the corresponding improved time-explicit algorithm. In Section 4, we use fictitious domain method to solve the equivalent variational problem of the relevant time discretization of wave equations. In Section 5, we use the Dirac delta function to improve the computation procedure of the space discretization equations. Finally, the results of numerical experiments and conclusion are presented in Sections 6 and 7.


## 2. Formulation of the Scattering Problem

We will discuss the scattering of monochromatic incident waves by perfectly conducting obstacle inR2 [1]. Let us consider a scattering body ω with boundary γ=∂ω, illuminated by an incident monochromatic wave of period T and incidence β. We bound Rn∖ωby an artificial boundary Γ. We denote by Ω the region of Rn between γ and Γ (see Figure 1).  The scattered field u satisfies the following wave equation and boundary conditions: (2.1)utt-Δu=0,inQ(=Ω×(0,T)),u=g,onσ(=γ×(0,T)),∂u∂n+∂u∂t=0,onΣ(=Γ×(0,T)),
where  g=-Re[e-ikteik(xcosβ+ysinβ)], withi=-1,k=2π/T.Figure 1
Computational Domain.Due to the periodic requirement,u also should satisfy
(2.2)u(0)=u(T),ut(0)=ut(T).Equation (2.1) represent the electric field u satisfying the two-dimensional Maxwell equation written in transverse magnetic (TM) form.


## 3. Exact Controllability and Least-Squares Formulations

Solving problem (2.1)-(2.2) is equivalent to finding a pair {v0,v1} such that
(3.1)u(0)=v0,ut(0)=v1,u(T)=v0,ut(T)=v1,
where u is the solution of (2.1). Problem (2.1), and (3.1) is an exact controllability problem which can be solved by the following controllability methodology given by [1].LetE is the space containing {v0,v1}(3.2)E=Vg×L2(Ω),
with Vg={φ∣φ∈H1(Ω),φ|γ=g(0)}.  Least-squares formulations of (2.1), and (3.1) are given by
(3.3)min⁡v∈EJ(v),
with
(3.4)J(v)=12∫Ω[|∇(y(T)-v0)|2+|yt(T)-v1|2]dx,∀v={v0,v1},
where y is the solution of
(3.5)ytt-Δy=0,inQ=(Ω×(0,T)),(3.6)y=g,onσ(=γ×(0,T)),(3.7)∂y∂n+∂y∂t=0,onΣ(=Γ×(0,T)),(3.8)y(0)=v0,yt(0)=v1.
The problem (3.3)–(3.8) may be solved by the conjugate algorithm [1]. Because this method looks some complicated, we use an alternative improved time-explicit asymptotic algorithm [4] to solve it. This method introduces an auxiliary parameter to control the time-explicit asymptotic iteration, and the auxiliary parameter is updated during the iteration based on the existing or current iterated solution of the wave equation. The algorithm is presented as follows.Algorithm 3.1.
We have the following steps.Step 1 : (initialization).
(1) Givenv={v0,v1}∈E as an initial guess.
(2)  compute the first periodic solutionyT: solving wave equation problem (3.5)–(3.8) to have solution yT={y(T),yt(T)}.
(3) compute the second periodic solutiony2T: solving wave equation problem (3.5)–(3.8) to get solution y2T={y(2T),yt(2T)} with initial condition y(0)=y(T),yt(0)=yt(T).
Step 2 : (computeβ* and update v,yT).
(1) Computeβ* by：(3.9)β*=12-J(yT)-J(v)∫Ω[|∇(δTTy(T))|2+|δTTyt(T)|2]dx,
where δTTy(T)=y(2T)-2y(T)+y(0).
(2) Updatev and yT by：(3.10)v={[v0+β*(y(T)-v0)],[v1+β*(yt(T)-v1)]},yT={[y(T)+β*(y(2T)-y(T))],[yT(T)+β*(yt(2T)-yt(T))]}.
Step 3 : (solve wave equation to obtainy2T).
Solve (3.5)–(3.8) for the second periodic solution y2T={y(2T),yt(2T)} with initial condition y(0)=y(T),yt(0)=yt(T).
Step 4 (test of the convergence).
Compute control functionJ(yT). If the value of J(yT) satisfies a given accuracy, then v=yT is taken as final solution, otherwise return to Step 2.


## 4. Fictitious Domain Method for Solving Wave Equation

Note that the above algorithm needs solve wave equations (3.5)–(3.8). The equivalent variational formulation of (3.5)–(3.7) is
(4.1)∫Ωyttzdx+∫Ω∇y⋅∇zdx+∫Γ∂y∂tdΓ=0,∀z∈V0,y=g,onσ,
where V0={φ∣φ∈H1(Ω),φ|γ=0}.The implementation used in [1] is based on an explicit finite difference scheme in time combined to piecewise linear finite element approximations for the space variables. Time discretization is carried out by a centered second-order difference scheme with time step Δt=T/N. After time discretization, (4.1) with (3.8) becomes(4.2)1Δt2∫Ω(yn+1-2yn+yn-1)zdx+∫Ω∇yn+1⋅∇zdx+12Δt∫Γ(yn+1-yn-1)zdΓ=0,∀z∈V0,yn+1=g(tn+1),onγ,y0=v0,y1-y-12Δt=v1.The fully discrete system can be obtained by the corresponding space discretization. BecauseΩ is irregular, if we directly use fitted meshes of Ω as in [1], we will meet great trouble of constructing meshes and difficulty of computation especially to those shape optimization problems with several scatters. So, we consider the problem (3.5)–(3.8) in the extended rectangular domain B=ω̅∪Ω with boundary Γ by the following boundary Lagrangian fictitious domain method. It allows the propagation to be simulated on B with uniform meshes. By introducing Lagrangian multipliers to enforce the Dirichlet boundary condition on γ, (3.5)–(3.8) is equivalent to the following variational problem.Find{y,λ}∈H1(B)×L2(γ),  such that
(4.3)∫Byttzdx+∫B∇y⋅∇zdx+∫Γ∂y∂tdΓ+∫γλzdγ=0,∀z∈H1(B),∫γμ(y-g)dγ=0,∀μ∈L2(γ),y(0)=v0,yt(0)=v1.LetΔt=T/N, discretize (4.3)  with respect to time with
(4.4)y0=v0,y0-y-1Δt=v1,
for n=0,1,…,N, we compute yn+1, λn+1 via the solution of(4.5)1Δt2∫B(yn+1-2yn+yn-1)zdx+∫B∇yn⋅∇zdx+12Δt∫Γ(yn+1-yn-1)zdΓ+∫γλn+1zdγ=0,∀z∈H1(B),(4.6)∫γμ(yn+1-gn+1)dγ=0,∀μ∈L2(γ).
Below, we consider conjugate gradient method for solving (4.5) and (4.6).For givenyn, yn-1, define linear functional f on H1(B)(4.7)f(z)=1Δt2∫B(-2yn+yn-1)zdx+∫B∇yn⋅∇zdx-12Δt∫Γyn-1zdΓ,∀z∈H1(B).
Let(4.8)a(w,z)=1Δt2∫Bwzdx+12Δt∫ΓwzdΓ,∀w,z∈H1(B).
Suppose z0 satisfies
(4.9)a(z0,z)+f(z)=0,∀z∈H1(B).
Then, (4.5) is
(4.10)a(yn+1-z0,z)+∫γλn+1zdγ=0,∀z∈H1(B).
Define A:L-1/2(γ)→L1/2(γ),Aμ=yμ|γ, for all μ∈L2(γ), where yμ satisfies
(4.11)a(yμ,z)+∫γμzdγ=0,∀z∈H1(B).
Let 〈·,·〉 denote scalar product in L2(γ), then
(4.12)a(yμ′,yμ)+〈μ′,Aμ〉=0∀μ′,μ∈L2(γ).-A is symmetric and positive definite. Then, in L2(γ) (4.5) (or (4.10)) becomes       (4.13)Aλn+1=(yn+1-z0)|γ.
By (4.6), (4.14)yn+1|γ=gn+1.Then,(4.15)Aλn+1=gn+1-z0|γ.
Its variational form is
(4.16)〈-Aλn+1,μ〉=〈z0|γ-gn+1,μ〉,∀μ∈L2(γ).A conjugate gradient algorithm for the solutionλn+1 of (4.16) is given by the following.Step 1 : (initialization).
(1) Give initial valueλ0∈L2(γ) and a real number ε>0 small enough.
(2) Findu0∈H1(B) such that
(4.17)a(u0,z)+f(z)+∫γλ0zdγ=0,∀z∈H1(B),
that is,
(4.18)1Δt2∫B(u0-2yn+yn-1)zdx+∫B∇yn⋅∇zdx+12Δt∫Γ(u0-yn-1)zdΓ+∫γλ0zdγ=0,∀z∈H1(B).
(3) Calculated0∈L2(γ) by
(4.19)∫γd0μdγ=∫γ(gn+1-u0)μdγ,∀μ∈L2(γ).
(4) Setw0=d0.Step 2.
For allk>0, calculate λk+1, dk+1, wk+1 from λk, dk, wk.
(1) Findu̅k∈H1(B) such that
(4.20)a(u̅k,z)+∫γwkzdγ=0,∀z∈H1(B),
that is,
(4.21)1Δt2∫Bu̅kzdx+12Δt∫Γu̅kzdΓ+∫γwkzdγ=0,∀z∈H1(B).
(2) Calculateρk:ρk=∫γ|dk|2dγ/-∫γu̅kwkdγ.
(3) Calculateλk+1: λk+1=λk-ρkwk.
(4) Calculate the new gradientdk+1∈L2(γ) by
(4.22)∫γdk+1μdγ=∫γdkμdγ+ρk∫γu̅kμdγ,∀μ∈L2(γ).Step 3 (test of the convergence).
If∥dk+1∥L2(γ)/∥d0∥L2(γ)≤ε, then takeλn+1=λk+1 and solve (4.10) for the corresponding solution yn+1, take yn+1 as the final solution; else, compute γk by
(4.23)γk=‖dk+1‖L2(γ)‖dk‖L2(γ),
and update wk by
(4.24)wk+1=dk+γkwk.
Set k=k+1, return to Step 2.


## 5. Improving the Computation Procedure of the Space Discretizations

Conventionally, we solve (4.18) and (4.21) by the finite element method (see [7–9]). In the computation procedure of the finite element discretizations, the mesh of the extended domain is regular, but the boundary is irregular. We will meet the trouble of computing the boundary integrals which leads to complex set operations like intersection and subtraction between irregular boundary γ and regular mesh of B. In order to avoid these difficulties and solve (4.18) and (4.21) more efficiently, we use the Dirac delta function to improve the computation procedure of the discretizations. We discuss this method as follows.We construct a regular Eulerian mesh onB(5.1)Bk={xij∣xij=(x0+ih,y0+jh),0≤i,j≤I},
where h is the mesh width (for convenience, kept the same both in x- and in y-directions). Assume that the configuration of the simple closed curve γ is given in a parametric form (s),0≤s≤L.  The discretization of the boundary γ employs a Lagrangian mesh, represented as a finite collection of Lagrangian points {Xk,0≤k≤M} apart from each other by a distance Δs, usually taken as being h/2. Let δ(·) be a Dirac delta function. In the following calculation procedure, δ is approximated by the distribution function δh. The choice here is given by the product
(5.2)δh(x)=dh(x1)dh(x2),
where x=(x1,x2) and dh is defined by
(5.3)dh(z)={0.25h[1+cos⁡⁡(πz2h)],|z|≤2h,0,|z|>2h.
Using the above Dirac delta function, we can transport the variational form (4.18) to the difference form. We write ∫γλ0zdγ in (4.18) as the following form:
(5.4)∫γλ0zdγ=〈L0,z〉H-1(B)H1(B),
where
(5.5)L0(x)=∫0Lλ0(s)δ(x-X(s))ds,∀x∈B,
that is, λ0 calculated over the Lagrangian points are distributed over the Eulerian points. Thus, we can write (4.18) in the difference form as follows:
(5.6)u0-2yn+yn-1Δt2-Δyn+L0=0,inB,∂yn∂n+u0-yn-12Δt=0,onΓ.
Thus, the solution of (4.18) is
(5.7)u0=2yn-yn-1+Δt2(Δyn-L0),inB,u0=yn-1-2Δt∂yn∂n,onΓ.
The discrete form of (5.5) is
(5.8)L0(xij)=∑kλk0δh(xij-Xk)Δs,∀xij∈Bh.
So, we can obtain u0(xij) for all xij∈Bh.In the same way, let(5.9)Wk(x)=∫0Lwk(s)δ(x-X(s))ds,∀x∈B.
Then, (4.21) also can be written in the difference form as follows:
(5.10)u̅k=-Δt2Wk,inB,u̅k=0,onΓ.
Calculate
(5.11)Wk(xij)=∑mwmkδh(xij-Xm)Δs,∀xij∈Bh.
Then, we can get u̅k(xij), for all xij∈Bh.Thus(5.12)u̅k|γ=u̅k(X(s))=∫Bu̅k(x)δ(x-X(s))dx,∀0≤s≤L.
Its discrete form is
(5.13)u̅mk=∑iju̅k(xij)δh(xij-Xm)h2,∀1≤m≤M.
And by (4.22), we have  (5.14)dk+1=dk+ρku̅k|γ.It can be seen from the above discretization process that most of the calculations are done over the Lagrangian points and the neighboring Eulerian points of the boundaryγ. The solutions of (4.18) and (4.21) are given explicitly by (5.7) and (5.10). And we only need do the evaluation in (5.8), (5.11), and (5.13) to obtain the solutions of (4.18) and (4.21). So, our method is easier to code and requires fewer computational operations than conventional finite element method (see [7–9]).


## 6. Numerical Experiments

In order to validate the methods discussed in the above sections, we apply our algorithm to simulate the scattering of planar monochromatic incident waves by a perfectly conducting obstacle. The obstacle is a Semiopen rectangular cavity; the internal dimensions of the cavity are4λ×1.4λ, and the thickness of the wall is 0.2λ as shown in Figure 2. Wavelength λ=0.25 m and incidence of illuminating waves is 0°. The corresponding scattered fields and convergence histories of control function J are shown in Figures 3 and 4.  Figures 3 and 4 show that our method performs as well as the method discussed in [7–9] does where fictitious domain method and obstacle fitted meshes were used.Figure 2
Semiopen rectangular cavity.Figure 3
Contours of the scattered field.Figure 4
Convergence histories.


## 7. Conclusions

In this paper, the fictitious domain technique is coupled to the improved time-explicit asymptotic method for calculating the time-periodic solutions of wave equations. It allows the propagation to be simulated on an obstacle free computational region with uniform meshes. One of the main advantages of the fictitious domain approach is that it is well suited to those shape optimization problems with several scatters that minimize, for example, a Rader Cross Section. We use the Dirac delta function to improve the computation procedure of space discretizations.  Numerical experiments invalidate that our algorithms are efficient and easy to implement alternative to more classical wave equation solvers.



---

*Source: 503791-2011-01-26.xml*