JP2014026524A - Simulation method and simulation program - Google Patents

Abstract

PROBLEM TO BE SOLVED: To provide a simulation method and a simulation program, for more stably performing simulation by an FDTD (Finite Difference Time Domain) method.SOLUTION: In a simulation method in one embodiment, sizes and a time step unit of a cell 22 of an FDTD area 18 are set, and an electromagnetic field is calculated according to an FDTD method. Electromagnetic fields are calculated by using first and second equations for electromagnetic field analysis with respect to an object model 22A and a medium model 22B of an analysis model 20, and the first equation for electromagnetic field analysis is derived by incorporating a recursive convolution method in the Maxwell equations. When lengths of first to third sides of a cell at least in the object model are denoted by Δs, Δs, and Δs, and the time step unit is denoted by Δt, Δs, Δs, Δsand Δt are set based on a stabilization condition: Δt/Δs≤0.354/v, where j=1,2,3, and v is the light velocity in vacuum.

Description

  The present invention relates to a simulation method and a simulation program using the FDTD method.

An electromagnetic wave simulation method that uses a computer to calculate the behavior of an electromagnetic wave traveling through a system in which a predetermined object is surrounded by a medium is useful for the following various designs.
・ Design of display materials such as diffusers and color filters ・ Design of electronic components such as high-frequency circuits, radars and antennas ・ Design of materials in the color industry such as ink, paint, plastic and dyeing ・ Remote sensing and meteorological science Design of various measuring devices in the field of medicine and medicine ・ Calculation of LED light extraction efficiency required for designing LED devices, especially LED devices containing nanoparticles or nanostructures ・ Light-electric conversion efficiency in solar cells Design for improvement

  In the simulation of electromagnetic waves, a simulation in a three-dimensional space that is closer to the progress of actual electromagnetic waves is required than in a simulation in a two-dimensional plane. Therefore, an electromagnetic wave simulation method is proposed in which a lattice cell is set in a virtual space in which an electromagnetic wave travels in a system surrounded by a medium, and calculation is performed using the FDTD (Finite Difference Time Domain) method. (See Patent Document 1). In the FDTD method, the FDTD region, which is the virtual space, is divided into a plurality of cells, and the Maxwell equation is directly solved by approximating the space and time differentiation in the Maxwell equation by a finite difference.

  However, if the object in the system contains metal, the real part of the dielectric constant of the object will be negative, so as the field update proceeds in the FDTD calculation, ie as the time step update proceeds Tend to diverge.

  To suppress such divergence, various analytical models for permittivity, such as Debye, Drude, or Lorentz, use frequency domain convolution, for example, Can be incorporated into FDTD. Such a method using convolution is known as a recursive convolution method (hereinafter also simply referred to as RC method).

JP 2010-204859 A

  Even if such an inductive convolution method is used, the simulation tends to become unstable depending on the setting of conditions such as a spatial discretization value.

  Therefore, an object of the present invention is to provide a simulation method and a simulation program capable of more stably performing a simulation by the FDTD method.

The simulation method according to one aspect of the present invention is applied to an analysis object in which an object having a real part of a dielectric constant of 0 or less and a positive value of an imaginary part is arranged in a medium having a positive real part of the dielectric constant. This is a simulation method for three-dimensionally analyzing an electromagnetic field when an electromagnetic wave is irradiated according to the FDTD method. In this simulation method, an FDTD region in which an analysis model of an analysis object is arranged, the cell size of the FDTD region composed of a plurality of cells, and a time step unit for updating an electromagnetic field in the FDTD method, An analysis condition setting step set by a computer, and an analysis step of calculating an electromagnetic field when electromagnetic waves are incident on the analysis model by a computer according to the FDTD method. In the analysis step, for the object model corresponding to the object in the analysis model, the electromagnetic field is calculated by a computer using the first electromagnetic field analysis formula derived from the Maxwell equation, and the medium is included in the analysis model. For the medium model corresponding to, the electromagnetic field is calculated by a computer using the second electromagnetic field analysis formula derived from the Maxwell equation. The first electromagnetic field analysis equation is derived by incorporating an inductive convolution method into the Maxwell equation, where the light velocity in vacuum is v, and at least the first to first orthogonal cells included in the object model are orthogonal to each other. 3 sides of the length of each Delta] s _1, and Delta] s ₂ and Delta] s _3, when the step unit was Δt time, the FDTD region setting _{step, Δs 1, Δs 2, Δs} 3 and Δt the following stabilization conditions It is set based on the formula (1) indicating
However, j = 1, 2, 3

In this method, at least Δs ₁ , Δs ₂ , Δs ₃ , which are cell sizes of the object model, and a time step unit Δt are set based on the stabilization condition expressed by Expression (1). As a result, the simulation can be stably performed.

  In one embodiment, the electromagnetic wave may be a monochromatic electromagnetic wave.

  In one embodiment, an actual value or an interpolated value obtained by interpolating the actual value may be used as the real part and imaginary part values of the dielectric constant in the object and the medium.

  In this way, it is possible to carry out the simulation more accurately by using the actual measurement value or the interpolation value.

  In one embodiment, the material of the object may be at least one of gold, silver, aluminum and copper.

  Another aspect of the present invention relates to a simulation program for causing a computer to execute the simulation method according to one aspect of the present invention.

  According to the present invention, simulation can be stably performed according to the FDTD method.

It is a schematic diagram of the analysis target to which one embodiment of the simulation method according to the present invention is applied. It is a schematic diagram of the FDTD area | region which arrange | positions the analysis model corresponding to the analysis target object shown in FIG. FIG. 3 is an enlarged view of the cell shown in FIG. 2. It is a block diagram of the analysis device for applying one embodiment of the simulation method concerning the present invention. It is a flowchart of one Embodiment of the simulation method which concerns on this invention. It is a figure which shows the comparison with the value calculated by the dispersion | distribution formula of Drude using the optimized parameter and the interpolation value. It is drawing which shows the result of an example of simulation.

  Hereinafter, embodiments of a simulation method and a simulation program according to the present invention will be described with reference to the drawings. In the description of the drawings, the same reference numerals are given to the same elements, and duplicate descriptions are omitted. The dimensional ratios and the like in the drawings do not necessarily match those described.

  The simulation method and the simulation program according to the present invention provide an electromagnetic field in a case where an electromagnetic wave is incident on an object to be analyzed in an object having a real part of permittivity of 0 or less in a medium having a positive real part of permittivity. This is for three-dimensional analysis according to the FDTD method (Finite Difference Time Domain Method). Examples of the object may include metal fine particles and a member in which a plurality of metal fine particles are connected. Examples of the medium may include a dielectric such as a resin, air, and vacuum. In the following description, unless otherwise specified, the object in the medium is a metal fine particle, and the medium is a dielectric. The present invention can be applied to, for example, the design of an optical sheet in which metal fine particles are used as a diffusing agent.

  FIG. 1 is a schematic diagram illustrating an example of an analysis object according to an embodiment. The analysis object 10 can be composed of metal fine particles 12 and a propagation medium 14 for propagating light around the metal fine particles 12. An example of the metal fine particle 12 is a nano fine particle composed of a metal. The metal fine particles 12 may be made of at least one of silver, gold, aluminum, and copper, for example. An example of the propagation medium 14 as a dielectric is resin.

  The simulation method and the simulation program according to the present embodiment use the FDTD method (Finite Difference Time) to calculate the electromagnetic field when monochromatic light (hereinafter also referred to as incident light) 16 is incident on the analysis target 10 as shown in FIG. For analysis according to Domain Method). The wavelength λ of the incident light 16 can be selected from the range of the ultraviolet region, the visible light region, and the near infrared light region. The wavelength λ is preferably selected from the range of 300 nm to 900 nm, more preferably from the range of 300 nm to 700 nm.

  The FDTD method will be described. In the FDTD method, an analysis model 20 corresponding to the analysis object 10 is arranged in an FDTD area 18 as a virtual space set by a computer. FIG. 2 is a schematic diagram showing an analysis model. In the analysis model 20, a metal fine particle model corresponding to the metal fine particle 12 is referred to as a metal fine particle model 20A, and a medium model corresponding to the medium 14 is referred to as a medium model 20B. In FIG. 2, the metal fine particle model 20 </ b> A is shown by a broken line in order to clearly show it. In FIG. 2, the outer shape of the metal fine particle model 20 </ b> A is represented by a curve, but the FDTD region 18 is spatially discretized as described later. Therefore, the outer shape of the metal fine particle model 20A is represented by a combination of fine straight lines (or cell sides). Such an object representation method is known as staircace approximation. If the discretization of the space becomes finer, the boundary surface of the object becomes closer to an actual curve. For this reason, in order to ensure accuracy when an object is expressed using cells (or grids) described later, fine spatial discretization is preferable. In FIG. 2, the size of the arrangement region of the analysis model 20 and the size of the FDTD region 18 are the same.

The FDTD method is a method of directly solving the Maxwell equation by approximating the space and time derivatives in the Maxwell equation shown in the equations (2a) and (2b) by a finite difference.
r is a position vector from the origin set in the FDTD area 18. μ is the magnetic permeability. In the analysis object 10 of this embodiment, the magnetic permeability μ is usually a constant. E is an electric field vector and H is a magnetic field vector. “∇” means x _u (∂ / ∂x) + y _u (∂ / ∂y) + z _u , where x _u , _yu and z _u are unit vectors in the x-axis direction, y-axis direction and z-axis direction. It is a vector differential operator defined by (∂ / ∂z). “×” is a rotation operator and represents a so-called outer product.

D in the equation (2b) is an electric flux density vector, and D is given by the following equation in the frequency domain.
In equation (3), ε (r, ω) is a dielectric constant.

In order to approximate and solve the Maxwell equation with spatial and temporal finite differences, the FDTD method divides the FDTD region 18 into a plurality of cells 22 as shown in FIG. The number of cells 22 may be set according to the size of the analysis model 20, the calculation time, the memory capacity required for the calculation, and the like. An example of the shape of each cell 22 constituting the FDTD region 18 is a cube. Accordingly, in the cell 22, the lengths of the x-axis direction side (first side), the y-axis direction side (second side), and the z-axis direction side (third side) that are orthogonal to each other are respectively expressed by Δx. When (= Δs ₁ ), Δy (= Δs ₂ ), and Δz (= Δs ₃ ), Δx = Δy = Δz. In the following description, the length of one side of the cubic cell 22 is Δs. In this case, Δx = Δy = Δz = Δs.

  FIG. 3 is an enlarged view of the cell 22 shown in FIG. In FIG. 3, a white arrow indicates an electric field, and a black arrow indicates a magnetic field. As shown in FIG. 3, the electric field strength is calculated at the center of each side constituting the cell 22, and the magnetic field strength is calculated at the center of the surface constituting the cell 22.

Next, formulas for calculating the electric field and the magnetic field will be described separately for calculation formulas for the metal fine particle model 20A and calculation formulas for the medium model 20B. In the following description, the dielectric constant is assumed to be ε. Complex number of the dielectric constant epsilon, when the real part epsilon _r and the imaginary part was epsilon _i, is ε = ε _r -iε _i. First, calculation formulas for the metal fine particle model 20A will be described.

A time step unit for electromagnetic field calculation is expressed as Δt, and n (n is an integer of 0 or more) is the number of time steps. In this case, a certain time t is t = nΔt. The electric field vector and the magnetic field vector at time t is expressed as ^{E n} and ^{H n.}

Since the propagation of electromagnetic waves is calculated in real time using the FDTD method, it is necessary to handle the electric flux density vector represented by Equation (3) in the time domain. The electric flux density vector D in the time domain is expressed by the following equation by performing Fourier transform on both sides of the equation (3) and applying the convolution theorem.
In Expression (4), ε (t) = F {ε (ω)}, and F represents a Fourier transform.

  When the dielectric constant ε (ω) is expressed by a special function having an angular frequency ω, the convolution integral of the equation (4) can be calculated recursively. This is known as an inductive convolution method (RC method).

In the present embodiment, a Drude primary approximation formula (Drude dispersion formula) is adopted as a function representing the dielectric constant ε (ω) of the metal fine particles 12. In this case, the dielectric constant ε (ω) is expressed as Equation (5).
In equation (5), ω _p is the plasma frequency and ν _c is the collision frequency. epsilon _∞ is the dielectric constant when the frequency is infinite, equal to 1. χ (ω) is an electrical specific susceptibility.

When χ (ω) is Fourier transformed, that is, F {χ (ω)} is expressed as χ (t), χ (t) is expressed by the following equation.
U (t) is a unit step function of time. When t = 0, U (t) = 0, and when t> 0, U (t) = 1.

At time t = nΔt (n is an integer of 0 or more) in the discretized time domain, the electric field vector E is a constant vector over a section [mΔt, (m + 1) Δt] (m is an integer of 0 to n), When the convolution integral of Formula (4) is evaluated, the following formula is obtained.
In Equation (7), ε ₀ is the dielectric constant of vacuum.

Equation (7) can be written as:
Ψ is a cumulative variable vector, and the following formula (9) is established.

Expression (9) is expressed as the following Expression (10).
In this formula (10),
It is. From equation (11)
Is established. In equation (12),
It is. In Equation (12), Ψ ⁰ = 0 and Ψ ¹ = E ¹ χ ⁰ .

From the equations (10) to (13b), the following inductive relational expression regarding the cumulative variable Ψ is obtained.

Equation (14) can be solved by setting χ ⁰ = 0 and setting ψ to 0 when n ≦ 0.

In the FDTD method, equations (2a) and (2b) are solved using the central difference method. When a second-order precision finite difference expression is used to approximate the first-order time derivative in equation (2a), equation (2a) is expressed as:

The update equation of the magnetic field vector H is derived from an equation in a differential format equivalent to the equation (2a) and is expressed as follows.
In Expression (16), d is a vector difference operator. d is defined by the following formula.
d _{x, d y} and _{d z} are first derivative ∂ / ∂x, a secondary accuracy finite difference approximation of ∂ / ∂y and ∂ / ∂z. When an arbitrary function is expressed by f (x, y, z), d _x f (x, y, z), _dy f (x, y, z), and d _z f (x, y, z) are , Are given by equations (18a) to (18c).

When the electromagnetic field is three-dimensionally calculated by the FDTD method, the magnetic field vector H has components in the x-axis direction, the y-axis direction, and the z-axis direction, respectively. The update equation for each component of the magnetic field vector H is expressed as follows in consideration of the equations (16), (17), and (18a) to (18c).

Next, the electric field vector E will be described. When a second-order precision finite difference expression is used to approximate the first-order time derivative in Equation (2b), Equation (2b) is expressed as:

The update equation of the electric field vector E is derived from an equation in a differential format equivalent to the equation (2b) and is expressed as follows.

Ψ ^{n + 1} −Ψ ^{n in} Expression (21) is given by Expression (14). D in Formula (21) is defined by Formula (17). d _x f (x, y, z), d _y f (x, y, z) and d _z f (x, y, z) are given by Expression (22a) to Expression (22c).

As in the case of the magnetic field vector H, the electric field vector E has components in the x-axis direction, the y-axis direction, and the z-axis direction, respectively. The update equation for each component of the electric field vector is expressed as follows in consideration of the equations (21), (17), and (22a) to (22c).

  Next, calculation formulas for the medium 14 will be described. When the magnetic permeability μ is a positive value and a constant, the update equation of the magnetic field vector H is the same as the equations (19a), (19b), and (19c).

On the other hand, in the dielectric, both the real part ε _r and the imaginary part ε _i of the dielectric constant ε are positive values. The update equation of each component of the electric field vector E for the medium model is expressed as follows without using the cumulative variable Ψ.

A in equation (24a) ~ formula (24c), when the σ and conductivity can be defined by _{(σΔt) / (2ε r)} = (ε i ωΔt) / (2ε r). However, ε _i ≧ 0 and ε _r > ε _i .

[Stabilization conditions]
Next, stabilization conditions for preventing the FDTD method from diverging will be described. For an algorithm that does not incorporate an inductive convolution method, that is, for the medium model 20B, a stable condition may be calculated based on a so-called Courant condition. However, since the stabilization condition for the algorithm incorporating the inductive convolution method is a stricter condition, the stabilization condition for the algorithm using the inductive convolution method will be described.

In order to stabilize the algorithm of the FDTD method incorporating the inductive convolution method, c ₁ and c ₂ in Equation (14) must satisfy the following conditions.

c ₁ and c ₂ can be calculated using the value of the dielectric constant of the metal microparticles for the wavelength used in the calculation. Here, with respect to the dielectric constant ε (= ε _r −iε _i ) of the metal fine particles 12 corresponding to the wavelength λ (or angular frequency ω) of the incident light 16, the plasma frequency ω _p and the collision frequency ν _c are expressed by the formula ( From 5), it is expressed by the following formula (26) and formula (27).
Therefore, the following expressions are established from the expressions (13a), (13b), (26), and (27).
Considering equations (28a) and (28b), c ₂ is positive and less than 1 as long as ε _i is greater than 0 and ε _r is 0 or less. That is, (Equation 25b) holds for metal fine particles in which ε _i is greater than 0 and ε _r is 0 or less. Furthermore, ωΔt can be selected to satisfy equation (25a).

  Next, a description will be given of stabilization conditions when using the equations (23a) to (23c) applied to the metal fine particle model 20A.

When the finite difference time differential operator is d _t , the equations (21) and (14) can be expressed as follows.
Here, d _t is an operator defined by d _t Ψ ^{n + 1/2} = Ψ ^{n + 1} −Ψ ⁿ .

From equations (29a) and (29b) and equations (19a) to (19c), the electromagnetic field update equation is expressed by the following equation using a matrix.
The matrices M and N are coefficient matrices. From the equation (30), the following equation is obtained.
In Equation (31), the matrix I is an identity matrix. The components of the matrices M and N are obtained as follows when the equations (29a), (29b) and the equations (19a) to (19c) are considered.
Here, _{p E} is represented by (1 / μ) (Δt / Δs) (d ×), a finite difference operator acting on the electric field vector E. p _H is represented by (1 / ε ₀ ) (Δt / Δs) (d ×), and is a finite difference operator acting on the magnetic field vector H.

Here, the matrix C is defined as follows.
When Expression (33) is used, Expression (31) is expressed by the following expression.
By using linear system theory with L _j and q _j (j = 1, 2, 3) as eigenvectors and eigenvalues of matrix C, equations (35) and (36) are obtained.
However, a _j is a constant.

In order to stabilize the algorithm, it is necessary that | q _j | ≦ 1. q _j is an eigenvalue of the matrix C. The eigenvalue is a solution of the eigen equation shown in Expression (37).
In Expression (37), the vector Q is a diagonal matrix defined by Expression (38).
Expression (39) is established from Expression (37).
One solution of equation (39) is zero. Expressing the other two solutions of Equation (39) as q _± , q _± is given by Equation (40).
In order to stabilize the algorithm, it is necessary that | q _± | ≦ 1. Since c ₁ and c ₂ are both 1 or less, if p _E p _H ≦ 2, then | q _± | ≦ 1.

The vector difference operator p _E p _H is given by equation (41).
Therefore, for the stability of the algorithm, equation (42) may be established.

Here, when the light speed of the vacuum (the speed of the electromagnetic wave) is v = 1 / (με ₀ ) ^1/2 , Expression (42) is expressed by Expression (43).

The maximum value of the operator “d × d ×” is evaluated by applying the operator “d × d ×” to the harmonic vector field E. The harmonic vector field E is expressed by the equation (44).
Here, E ₀ is a vector represented by E ₀ = x _u E _0x + y _u E _0y + z _u E _0z . As described above, x _u , _yu , and z _u are unit vectors in the x-axis direction, the y-axis direction, and the z-axis direction, respectively. In Equation (44), k is a wave vector and r is a position vector.

In three dimensions, d × E and d × d × E are given by Equation (45a) and Equation (45b), respectively.

The following equation holds from the definition of the first-order difference operator.
In formula (46), α and β are independent subscripts. α and β are x, y, or z, respectively.

When the x component of d × d × E is represented as (d × d × E) _x , (d × d × E) _x is given by the following equation.

Assuming that the polar angle and azimuth angle with respect to the direction of the vector E ₀ in the vector field of interest are θ and φ, E _0x , E _0y, and E _0z are respectively expressed by equations (48a), (48b), and (48c). expressed.

| sinθcosφ | ≦ 1, | sinθsinφ | ≦ 1 and | cosθ | ≦ 1 a with _{which, | sin (k m Δs /} 2) | since it is ≦ 1, maximum value is 16 right-hand side of equation (47), That is, the x component of d × d × E is (d × d × E). The maximum value of _x can be evaluated as 16. Similarly, the maximum value of the y component and the z component of d × d × E can be evaluated as 16.

Therefore, the following expression is established from Expression (43).
From equation (49), it can be seen that the maximum value of the ratio of discretization of time and space is 0.354 / v (where v is the speed of light in a vacuum) in order to stabilize the algorithm.

  Therefore, when using the equations (23a) to (23c) applied to the metal fine particle model 20A, Δt and Δs may be set so as to satisfy the equation (49). As described above, the stabilization condition for the metal fine particle model 20A is more severe than that for the medium model 20B. Therefore, if Δt and Δs are set so as to satisfy the stabilization conditions for the expressions (23a) to (23c), the medium model 20B can be stably calculated.

Next, the dielectric constant value used in the FDTD method will be described. In the update equation (first equation for electromagnetic field analysis) of the electric field vector E and the magnetic field vector H for the exemplified metal fine particle model 20A, the plasma frequency ω _p and the collision frequency ν _c in the equations (13a) and (13b). Includes the dielectric constant ε of the metal fine particles 12.

More specifically, with respect to the dielectric constant ε (= ε _r −iε _i ) of the metal fine particles 12 corresponding to the wavelength λ (or frequency) of the incident light 16, the plasma frequency ω _p and the collision frequency ν _c are: It represents with Formula (26) and Formula (27).

In the update equation (second electromagnetic field analysis equation) of the electric field vector E and the magnetic field vector H for the medium model, the real part ε of the dielectric constant ε of the medium 14 is represented by “a” in the equations (24a) to (24c). _r and imaginary part ε _i are included.

In the present embodiment, actual values or values obtained by interpolating these actual values (hereinafter referred to as interpolation values) are used as the values of the real part ε _r and the imaginary part ε _i of such a dielectric constant ε. In the present embodiment, the actual measurement values are known actual measurement values such as data values described in books listed as a list such as optical characteristics, or individually measured for the metal fine particles 12 actually used for simulation. It can be a value. Examples of books listed as a list of optical properties are `` CLFoiles, 'Optical properties of pure metals and binary alloys', Chapter 4 of Landolt-Bornstein Numerical Data and Functional Relationships in Science and Technology New Series, Vol. 15, Subvolume b, KHHellwege and JLOlsen, Ed. Springer-Verlag, Berlin 1985, pp. 228 ”. An example of an interpolation method for obtaining an interpolation value is cubic spline interpolation.

  Next, an analysis device 24, a simulation method, and a simulation program for executing the above algorithm will be described.

  FIG. 4 is a block diagram of the analysis device 24 for carrying out an example of the simulation method according to the present invention.

  The analysis device 24 shown in FIG. 4 is a computer device such as a personal computer device or a workstation, and includes an input unit 26, a memory 28, a central processing unit 30, and an output unit 32 as main components.

The input means 26 is a keyboard or the like, and receives input data including parameters for performing a simulation. The parameters include various conditions regarding the incident light 16 including the wavelength λ of the incident light 16, values of ε _r and ε _i of the metal fine particles 12 and the medium 14 with respect to the wavelength λ, various conditions of the analysis target object 10, and boundary conditions, And the maximum time t _max (t _max = nΔt) in the simulation.

The various conditions of the analysis object 10 may be any conditions as long as the analysis model 20 corresponding to the analysis object 10 can be virtually constructed in the analysis device 24. Examples of various conditions of the analysis object 10 are the size (size) of the medium 14 and the metal fine particles 12 and the position of the metal fine particles 12 in the medium 14. The values of the real part ε _r and the imaginary part ε _i of the dielectric constant ε are actually measured values (known values) or interpolation values as described above. As the boundary condition, various absorption boundary conditions or various periodic boundary conditions can be used. Examples of absorbing boundary conditions may include secondary precision Mur boundary conditions, PML (Perfectly matched layer) and UPML (Uniaxial PML) boundary conditions.

  The memory 28 includes various programs for operating the analysis device 24, input data from the input means 26, FDTD program (simulation program) 34 and FDTD program for performing FDTD analysis according to a simulation method described later. 34 stores the calculation result during or after execution. The FDTD program 34 executes an FDTD method and an analysis condition setting module 34A for setting analysis conditions for stably performing the FDTD method by using update equations for application to the metal fine particle model and the medium model, respectively. And an analysis module 34B.

  The central processing unit 30 reads the FDTD program 34 stored in the memory 28 and executes various calculations for simulation. The central processing unit 30 also has a function of reading various programs stored in the memory 28 and controlling components of the analysis device 24.

  The output means 32 is a display etc., for example, and displays the simulation result obtained by being calculated by the central processing unit 30.

  Next, a simulation method according to the FDTD method will be described. FIG. 5 is a flowchart showing a simulation method. The simulation is performed by the central processing unit 30 reading and executing the FDTD program 34 stored in the memory 28 of the analysis device 24.

  As shown in FIG. 5, in step S10, the operator inputs simulation parameters via the input means 26 (input process). When the input means 26 receives a parameter, the memory 28 stores the parameter.

  Next, in step S12, an analysis condition is set (analysis condition setting step). In this analysis condition determination step, the central processing unit 30 sets analysis conditions based on the stabilization conditions. The analysis condition can be set based on the formula (49), for example. Specifically, when the time step unit Δt is constant, Δs from Expression (49), that is, the size of one side of the cell 22 when the analysis model 20 for the analysis target 10 is divided into a plurality of cells 22 is set. obtain. When Δs satisfying equation (49) is set, Δs that can reduce the calculation time is set with higher priority.

  In other embodiments, Δt and Δs may be arbitrarily set so as to satisfy Equation (49).

  In yet another embodiment, Δs and Δt may be set by the following steps.

Step I: Set v in formula (49) to 2.998 × 10 ⁸ m / s.
Step II: Δt is selected so that c ₁ and c ₂ represented by formula (28a) and formula (28b) are 1 or less, respectively.
Step III: Δs is selected so as to satisfy Expression (49).

  In subsequent step S14, the FDTD region 18 defined by the spatial discretization condition set in step S12 is set, and the analysis model 20 is arranged in the FDTD region 18.

  Thereafter, in step S16, the electromagnetic fields of all the cells 22 constituting the FDTD region 18 are initialized. Specifically, the component in each axial direction of the electric field and magnetic field in each cell 22 is set to zero.

  In step S18, the magnetic field and electric field at time t are calculated in all the cells 22, and time t is updated (analysis process). Specifically, for the medium model 20B, the magnetic field and the electric field are calculated based on the equations (19a) to (19c) and (24a) to (24c), and for the metal fine particle model 20A. The magnetic field and the electric field are calculated based on the equations (19a) to (19c) and (23a) to (23c). As shown in Formula (19a) to Formula (19c), Formula (23a) to Formula (23c), and Formula (24a) to Formula (24c), the FDTD method calculates the electric field when calculating the electric field. At the position, the electric field is calculated based on the electric field calculated one unit before the time step and the magnetic field calculated before the half time step unit in the vicinity of the position where the electric field is calculated (a position surrounding the position). . On the other hand, when calculating the magnetic field, at the position where the magnetic field is calculated, the magnetic field calculated one time step before, and the vicinity of the position where the magnetic field is calculated (position surrounding the position) half a time step unit before A magnetic field is calculated based on the calculated electric field.

In step S20, it is determined whether or not the updated time t is the maximum time _tmax . If the time t is smaller than the maximum time t _max (“YES” in step S20), the process returns to step S18.

In step S20, when the time t is equal to or greater than the maximum time _tmax (in the case of “NO” in step S20), the process proceeds to step S22, and the analysis result is output by the output means 32 (output process).

In the simulation method and the FDTD program 34 of the present embodiment described above, Δs and Δt are set based on the stabilization condition. Therefore, the FDTD method can be stably executed even if the analysis object includes the metal fine particles 12 having a negative value of the real part ε _r of the dielectric constant ε.

As the values of the real part ε _r and the imaginary part ε _i of the dielectric constant ε of the metal fine particles 12 included in the expressions (26) and (27), an actual measurement value (known value) with respect to the wavelength λ used in the analysis is used. Therefore, it is possible to simulate more accurately.

  However, the measured value of the dielectric constant ε is normally obtained discretely in the energy or frequency of light. Since the wavelength is inversely proportional to the frequency, even if there are actually measured values of the dielectric constant ε at regular intervals in the energy region, there may be cases where the actual measured values do not exist at regular intervals in the wavelength region. Therefore, there may be a case where there is no measured value of the dielectric constant ε with respect to the wavelength λ used for the analysis. In this case, as described above, the simulation is performed using the interpolation value obtained from the actually measured value using the interpolation method. Thereby, in this embodiment, it is possible to simulate more accurately. An example of the interpolation method is cubic spline interpolation.

When there is no actual measurement value for a certain wavelength λ as the values of the real part ε _r and the imaginary part ε _i of the dielectric constant ε, for example, the dielectric constant ε is calculated using the Drude dispersion equation with optimized parameters. Is also possible.

  However, as shown in FIG. 6, the dielectric constant ε can be calculated more accurately by using the interpolation value.

FIG. 6 is a drawing showing a cubic spline interpolation value of the real part of the dielectric constant of silver. For comparison, FIG. 6 also shows values calculated by Drude's dispersion equation using optimized parameters. The horizontal axis in FIG. 6 indicates the wavelength (nm), and the vertical axis indicates the real part (F / m) of the dielectric constant. A mark “×” in FIG. 6 indicates an actual measurement value. Specifically, “CLFoiles,“ Optical properties of pure metals and binary alloys ”, Chapter 4 of Landolt-Bornstein Numerical Data and Functional Relationships in Science and Technology New Series, Vol. 15, Subvolume b, KHHellwege and JLOlsen, Ed.Springer -Verlag, Berlin 1985, pp. 228 ". A mark “◯” in FIG. 6 is an interpolation value directly calculated from an actual measurement value using cubic spline interpolation. The mark “●” in FIG. 6 is a value obtained from the first-order dispersion equation of Drude using optimized Drude parameters (plasma frequency ω _p and collision frequency ν _c ).

Optimization parameters in the Drude model are genetic algorithms (eg, “S. Banerjee and LN Hazra,“ Experiments with a genetic algorithm for structural design of cemented doublets with prespecified aberration targets ”, App. Opt. Vol. 40, No. 34). , 2001, pp. 6265.). Specifically, the merit function including the real part and the imaginary part of the dielectric constant represented by the plasma frequency ω _p and the collision frequency ν _c that are parameters in the Drude model is optimized using a genetic algorithm. The values indicated by the mark “●” in FIG. 6 are values calculated by setting the plasma frequency ω _p to 1.288 × 10 ¹⁶ s ⁻¹ and the collision frequency ν _c to 6.4713 × 10 ¹³ s ⁻¹ .

  As can be understood from FIG. 6, the tendency obtained from the actual measurement value (for example, the fitting curve of the actual measurement value) from the value calculated by the Drude's dispersion formula including the optimization parameter when the cubic spline interpolation is used. It can be understood that the values are closer. In particular, on the short wavelength side, the value using cubic spline interpolation tends to be closer to the tendency obtained from the actual measurement value.

  Therefore, when there is an actual measurement value corresponding to the wavelength λ, the actual measurement value is adopted. When there is no actual measurement value corresponding to the wavelength λ, the simulation is performed by adopting an interpolation value using cubic spline interpolation. It can be done more accurately.

FIG. 7 is a diagram illustrating an example of a simulation result. FIG. 7 shows the result of simulating light scattering by the metal fine particle model 20A in the vicinity of the metal fine particle model 20A, and shows the simulation result of the xz plane. In FIG. 7, the shading indicates the light intensity, and the brighter the light intensity is. The simulation conditions shown in FIG. 7 are as follows.
-Material of metal fine particles 12: Silver-Particle size of metal fine particles 12: 20 nm
-Magnetic permeability μ of metal fine particles 12: 4π × 10 ⁻⁷ [H / m]
Medium 14: Vacuum (that is, metal fine particles 12 are arranged in a vacuum).
The dielectric constant ε of the medium 14: ε = ε ₀ = 1 / μv ^{2 In} this simulation, v is the speed of light in a vacuum (2.998 × 10 ⁸ m / s).
Magnetic permeability μ of medium 14: 4π × 10 ⁻⁷ [H / m]
The shape of the cell 22: a cube The length Δs of one side of the cell 22: 4 nm
-Wavelength λ of incident light 16: 440 nm
-Number of cells in the x-axis direction: 101
-Number of cells in the y-axis direction: 101
-Number of cells in z-axis direction: 131
Time step unit Δt: 4.72 × 10 ⁻¹⁶ s
-Number of steps n: 18660

  Since the dielectric constant for silver having a wavelength of 440 nm is not described at least in the literature cited in the explanation for the dielectric constant, the dielectric constant was calculated by cubic spline interpolation. The calculated dielectric constant was −6.511 + 0.193i. Under the simulation conditions, Expression (34) that is a stabilization condition is satisfied.

  As shown in FIG. 7, the electromagnetic field in the vicinity of the metal fine particles 12 can be calculated. That is, it can be understood that the FDTD method can be stably executed.

  As mentioned above, although various embodiment of this invention was described, this invention is not limited to the said embodiment, A various deformation | transformation is possible in the range which does not deviate from the meaning of this invention.

  For example, in the form shown in FIG. 2, the size of the cell 22 is the same in the metal fine particle model and the medium model, and satisfies the formula (49). However, the sizes of the cell 22 for the metal fine particle model and the cell 22 for the medium model may be different. Specifically, Δs is set so as to satisfy the equation (49) for the cell 22 for the metal fine particle model, and for the cell 22 for the medium model, the equations (19a) to (19c) and You may set (DELTA) s so that the stabilization conditions (for example, Courant conditions) with respect to Formula (24a)-Formula (24c) may be satisfy | filled. Thus, when the size of the cell 22 is set for each of the metal fine particle model and the medium model, the Δs of the cell 22 with respect to the medium model is usually larger than the Δs of the cell 22 with respect to the metal fine particle model. This is because the stabilization conditions of the cell 22 with respect to the metal fine particle model are more severe. Therefore, when the size of the cell 22 is set for each of the metal fine particle model and the medium model, the number of the cells 22 is reduced as compared with the case where the FDTD region 18 of the analysis model 20 is divided by Δs of the cell 22 with respect to the metal fine particle model. Is done. That is, it is possible to shorten the processing time for calculating the electromagnetic field by the FDTD method.

The shape of the cell 22 is a cube, but the shape of the cell 22 may be a rectangular parallelepiped shape. When the shape of the cell 22 is a rectangular parallelepiped shape, Δx = Δs, which is the length of first to third sides (sides extending in the x-axis, y-axis, and z-axis directions) constituting the cell 22 that are orthogonal to each other. ₁ , Δy = Δs ₂ and Δz = Δs ₃ may be set so as to satisfy the following formulas (50a) to (50c).

  Furthermore, Δt may be a value set in advance or may be determined while changing Δt and Δs so as to satisfy Equation (49) or Equation (50a) to Equation (50c). Moreover, although incident light (electromagnetic wave) was made into a single color, it is not limited to a single color. Furthermore, the metal fine particles are exemplified as the object disposed in the medium. However, the real part of the dielectric constant is 0 or less and the imaginary part is a positive value, and has a three-dimensional shape. If it is.

  DESCRIPTION OF SYMBOLS 10 ... Analysis object, 12 ... Metal fine particle (object), 14 ... Propagation medium (medium), 16 ... Incident light (light), 18 ... FDTD area | region, 20 ... Analysis model, 22 ... Cell, 34 ... Simulation program.

Claims (5)

An electromagnetic field when an electromagnetic wave is incident on an analysis object in which an object having a real part of permittivity of 0 or less and an imaginary part having a positive value is placed in a medium having a positive real part of permittivity, A simulation method for three-dimensional analysis according to the FDTD method,
A FDTD region in which an analysis model of the analysis object is arranged, the size of the cell of the FDTD region composed of a plurality of cells and a time step unit for updating the electromagnetic field in the FDTD method, An analysis condition setting step set by a computer;
An analysis step of calculating an electromagnetic field when the electromagnetic wave is incident on the analysis model by the computer according to the FDTD method;
With
In the analysis step, for the object model corresponding to the object among the analysis models, an electromagnetic field is calculated by the computer using a first electromagnetic field analysis expression derived from Maxwell's equations, and the analysis For the medium model corresponding to the medium among the models, the electromagnetic field is calculated by the computer using the second electromagnetic field analysis formula derived from the Maxwell equation,
The first electromagnetic field analysis equation is derived by incorporating an inductive convolution method into the Maxwell equation,
When the speed of light in vacuum is v, the lengths of at least the first to third sides of the cell in the object model are Δs ₁ , Δs _2, and Δs ₃ , respectively, and the time step unit is Δt In the analysis condition setting step, Δs ₁ , Δs _2, Δs ₃ , and Δt are set based on Expression (1) indicating the following stabilization condition:
Simulation method.
However, j = 1, 2, 3   The simulation method according to claim 1, wherein the electromagnetic wave is a monochromatic electromagnetic wave.   The simulation method according to claim 1, wherein an actual value and an interpolated value obtained by interpolating the actual value are used for the real part and imaginary part values of the dielectric constant in the object and the medium.   The simulation method according to claim 1, wherein the material of the object is at least one of gold, silver, aluminum, and copper.   The simulation program which makes a computer perform the simulation method as described in any one of Claims 1-4.