- Email: [email protected]
Modelling and simulation of electric and magnetic fields in homogeneous non-dispersive anisotropic materials
Available online at www.sciencedirect.com
Computers and Structures 85 (2007) 1623–1634 www.elsevier.com/locate/compstruc
Modelling and simulation of electric and magnetic fields in homogeneous non-dispersive anisotropic materials V.G. Yakhno a
a,*
, T.M. Yakhno
b
Department of Mathematics, Dokuz Eylul University, Fen-Edebiyat Fak., Kaynaklar, Buca, Izmir, Turkey b Department of Computer Engineering, Dokuz Eylul University, Kaynaklar, Buca, Izmir, Turkey Received 5 December 2006; accepted 15 February 2007 Available online 24 April 2007
Abstract The paper describes a new efficient method for finding and simulation of electric and magnetic fields in homogeneous non-dispersive electrically anisotropic materials. Electromagnetic wave propagation in these materials is descried by the time-dependent Maxwell’s system which implies initial value problems (IVPs) for electric and magnetic fields. Our method essentially uses matrix symbolic calculations in MATLAB and allows us to compute explicit formulae of electric and magnetic fields. These formulae are used for generating images of electric and magnetic fields inside different materials. A collection of electric and magnetic field images is presented in the paper. Ó 2007 Elsevier Ltd. All rights reserved. Keywords: Electric and magnetic fields; Anisotropy; Wave propagation; Simulation; Matrix symbolic computations; MATLAB
1. Introduction 1.1. General background information Anisotropic materials attract a great deal of interest in modern engineering applications, in particular, microelectronics [1–3]. In general case a type of anisotropy specifies properties of materials, and electromagnetic waves propagation inside these materials essentially depends on the direction, generating electromagnetic fields of very peculiar forms. Modelling and simulation of waves inside these materials constitute an important interdisciplinary area of research with many cutting-edge scientific and technological applications. Fundamental models for electromagnetic wave propagation used in science and engineering are based on the Maxwell’s system [1–4]. In particular the mathematical models of electric and magnetic waves motions in electrically and magnetically anisotropic materials are very well known *
Corresponding author. Tel.: +90 0232 4128592; fax: +90 0232 4534188. E-mail address: [email protected] (V.G. Yakhno). 0045-7949/$ - see front matter Ó 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.compstruc.2007.02.019
(see, for example, [1–5]). These mathematical models are related to the time-dependent Maxwell’s system containing matrices of the dielectric permittivity E and magnetic permeability M. We note (see, for example, [1]) that a medium is electrically anisotropic if it is described by E ¼ ðeij Þ33 and M ¼ lI; and magnetically anisotropic if it is described by E ¼ eI and M ¼ ðlij Þ33 , where ðeij Þ33 , ðlij Þ33 are arbitrary matrices, e, l are positive constants, I is the identity matrix. A medium can be both electrically and magnetically anisotropic if E ¼ ðeij Þ33 and M ¼ ðlij Þ33 . 1.2. Permittivity, permeability and crystallographic structure of anisotropic materials A smallest block (three-dimensional array of atoms) of anisotropic materials is determined by repeated replication in three dimensions. Its symmetry tells how the constituent atoms are arranged in a regular repeating configuration. The structure of these three-dimensional ‘unit cell’ of atoms in anisotropic materials may have one of seven basic shapes: cubic, hexagonal, tetragonal, trigonal, orthorhombic, monoclinic and triclinic, see [6]. For linear homogeneous
1624
V.G. Yakhno, T.M. Yakhno / Computers and Structures 85 (2007) 1623–1634
non-dispersive electrically and magnetically anisotropic materials the electric and magnetic flux densities D and B have the following relations with the electric and magnetic fields E and H: D ¼ EE, B ¼ MH, where the permittivity E and the magnetic permeability M are symmetric positive definite matrices. In our paper, we will study homogeneous non-dispersive electrically anisotropic materials for which E is a symmetric positive definite matrix and M ¼ lI, l > 0 is a constant. The relations of the crystallographic structures of these materials and their dielectric permittivity tensors E are the following: E ¼ diagðe11 ; e11 ; e11 Þ for the cubic structure; E ¼ diagðe11 ; e11 ; e33 Þ for the hexagonal, tetragonal and trigonal structures; E ¼ ðeij Þ33 , eij ¼ eji , e13 ¼ e23 ¼ 0 for the monoclinic structure; and an arbitrary symmetric positive definite matrix E ¼ ðeij Þ33 for the triclinic structure of materials. 1.3. Basic equations and problems set-up The propagation of electromagnetic waves in homogeneous non-dispersive electrically anisotropic materials is described by the following time-dependent Maxwell’s system with a symmetric positive definite matrix of dielectric permittivity E ¼ ðeij Þ33 and the magnetic permeability M ¼ lI (see, for example, [1,2]). Let x ¼ ðx1 ; x2 ; x3 Þ be a space variable from R3, t be a time variable from R, then Maxwell’s system is given by the following relations: oE þ j; ot oH ; curlx E ¼ M ot
ð1Þ
curlx H ¼ E
ð2Þ
divx ðMHÞ ¼ 0;
ð3Þ
divx ðEEÞ ¼ q;
ð4Þ
where E ¼ ðE1 ; E2 ; E3 Þ, H ¼ ðH 1 ; H 2 ; H 3 Þ are electric and magnetic fields, Ek ¼ Ek ðx; tÞ; H k ¼ H k ðx; tÞ; k ¼ 1; 2; 3; j ¼ ðj1 ; j2 ; j3 Þ is the density of the electric current, jk ¼ jk ðx; tÞ; k ¼ 1; 2; 3; M is the magnetic permeability, E is the dielectric permittivity, q is the density of electric charges. The operators curlx and divx are defined by oE3 oE2 oE1 oE3 oE2 oE1 curlx E ¼ ; ; ; ox2 ox3 ox3 ox1 ox1 ox2 oE1 oE2 oE3 þ þ : divx E ¼ ox1 ox2 ox3 The conservation law of charges is given by (see, for example, [2]) oq þ divx j ¼ 0: ot
ð5Þ
We suppose that E ¼ 0; H ¼ 0; q ¼ 0
for t 6 0;
ð6Þ
this means there is no electric charges at the time t 6 0; electric and magnetic fields vanish for t 6 0.
Differentiating (1) with respect to t and using (2) and (6) we find curlx ðM1 curlx EÞ ¼ E Ejt60 ¼ 0;
o2 E oj þ ; ot2 ot
ð7Þ ð8Þ
where M1 is the inverse matrix to M. In this paper, we assume that M ¼ I, E ¼ ðeij Þ33 is a given symmetric positive definite matrix with constant elements and j is taken in the form jðx; tÞ ¼ edðx1 Þdðx2 Þdðx3 ÞdðtÞ; x ¼ ðx1 ; x2 ; x3 Þ 2 R3 , where e 2 R3 is a given unit vector, dð Þ is the Dirac delta function. This form of the current density is related to a directional impulse electric source concentrated at the origin of coordinates at the time t ¼ 0 in the direction e. As a result of it we conclude that the Maxwell’s system (1)–(4) with conditions (6) implies the following initial value problem for the electric field (IVPEF) and initial value problem for the magnetic field (IVPMF). IVPEF. Let E be a given symmetric positive definite matrix, M ¼ I. Find the vector function E satisfying (7) and (8). IVPMF. Let E be a solution of IVPEF. Find the vector function H satisfying (2) and the condition Hjt60 ¼ 0. 1.4. Information on the previous studies and an assessment of the present results For the electric field equation (7) different problems for isotropic and anisotropic materials have been studied and many applications have been made, see, for example, [1,2,4–17]. In particular, decomposition method for the case of isotropic materials (E is a diagonal matrix of the form E ¼ diagðe11 ; e11 ; e11 Þ) has been suggested in [7]. Analytic methods of Green’s functions constructions have been studied for the case of isotropic materials in [8,9]; for uniaxial anisotropic media ðE ¼ diagðe11 ; e11 ; e33 ÞÞ in [10,11]; for biaxial anisotropic crystals ðE ¼ diagðe11 ; e22 ; e33 ÞÞ in [12]; for arbitrary non-dispersive homogeneous anisotropic dielectrics (E ¼ ðeij Þ33 is a symmetric positive definite matrix) in [13]. Modelling lossy anisotropic dielectric waveguides with the method of lines has been made for inhomogeneous biaxial anisotropic media in [14]. Most of the studies and modelling electromagnetic waves had been made by numerical methods, in particular finite element method [4,5,15,16]. On the other hand nowadays computers can perform very complicated symbolic computations (in addition to numerical calculations) and this opens up new possibilities in modelling and simulation of wave propagation phenomena. Symbolic computations can be considered as a useful tool for analytical methods which can provide exact solutions of problems [18–20]. Unfortunately the exact solutions cannot be found for all complex equations and systems. But when the exact solution can be found it leads to the significant simplification of modelling and simula-
V.G. Yakhno, T.M. Yakhno / Computers and Structures 85 (2007) 1623–1634
tion and allows engineers to focus on the ideas rather than on overcoming calculational difficulties. In this paper, we present a new analytical method for constructing explicit formulae for electric and magnetic fields simultaneously inside different electrically anisotropic non-dispersive homogeneous materials. In general case these explicit formulae are very cumbersome and they have been computed using symbolic computations in MATLAB. Applying these explicit formulae simulations of magnetic waves are obtained. Images of magnetic field components from a pulse current concentrated in a fixed point are presented in the paper. The pulse pointed current was chosen because this source is convenient for modelling fronts of electric and magnetic waves. Presented images of field components give an opportunity to observe the traces of electric and magnetic wave fronts and the distribution of energy in materials with different anisotropic structure. The paper is organized as follows. In Section 2, we describe in details a procedure of finding exact formulae of electric and magnetic fields for arbitrary type of anisotropy. Section 3 contains images of the traces of electric and magnetic fields for different anisotropic materials classified in the introduction. The last section contains conclusions and directions of a future research. 2. Finding exact formulae for electric and magnetic fields
1625
IVPMF can be written in terms of Fourier images e tÞ, Hðm; e tÞ as follows: Eðm; e 1 ðm; tÞ oH e 3 ðm; tÞ m3 E e 2 ; ¼ i½m2 E ot e 2 ðm; tÞ oH e 1 ðm; tÞ m1 E e 3 ðm; tÞ; ¼ i½m3 E ot e 3 ðm; tÞ oH e 2 ðm; tÞ m2 E e 1 ðm; tÞ; ¼ i½m1 E ot fj ðm; tÞjt60 ¼ 0; j ¼ 1; 2; 3: H Integrating (12)–(14) under conditions (15) we find Z t e e 3 ðm; sÞ m3 E e 2 ðm; sÞds; H 1 ðm; tÞ ¼ i ½m2 E 0 Z t e 1 ðm; sÞ m1 E e 3 ðm; sÞds; e H 2 ðm; tÞ ¼ i ½m3 E 0 Z t e 2 ðm; sÞ m2 E e 1 ðm; sÞds: e 3 ðm; tÞ ¼ i ½m1 E H
ð12Þ ð13Þ ð14Þ ð15Þ
ð16Þ ð17Þ ð18Þ
0
One of the mainRpoints of constructing explicit formulae e ðm; sÞds; j ¼ 1; 2; 3 is a procedure of e j ðm; tÞ and t E for E 0 j finding explicitly a non-singular matrix TðmÞ and diagonal matrix DðmÞ with non-negative elements such that TT ðmÞETðmÞ ¼ I;
ð19Þ
T
T ðmÞSðmÞTðmÞ ¼ DðmÞ;
ð20Þ T
2.1. IVPEF and IVPMF in terms of Fourier transform
where I is the identity matrix, T ðmÞ is the transposed matrix to TðmÞ.
e tÞ, Hðm; e tÞ be the Fourier transform images of Let Eðm; the electric and magnetic fields Eðx; tÞ and Hðx; tÞ with respect to x ¼ ðx1 ; x2 ; x3 Þ 2 R3 , i.e.
2.2. Procedure of finding matrices TðmÞ, DðmÞ and its implementation in MATLAB
e tÞ ¼ ð E e 1 ðm; tÞ; E e 2 ðm; tÞ; E e 3 ðm; tÞÞ; Eðm; e tÞ ¼ ð H e 2 ðm; tÞ; H e 3 ðm; tÞÞ; e 1 ðm; tÞ; H Hðm; e j ðm; tÞ ¼ Fx ½Ej ðm; tÞ; E
e j ðm; tÞ ¼ Fx ½H j ðm; tÞ; H
where j ¼ 1; 2; 3; m ¼ ðm1 ; m2 ; m3 Þ 2 R3 ; Z þ1 Z þ1 Z þ1 Fx ½Ej ðm; tÞ ¼ Ej ðx; tÞeimx dx1 dx2 dx3 ; 1
1
xm ¼ x1 m1 þ x2 m2 þ x3 m3 ;
1
i2 ¼ 1:
IVPEF (7) and (8) can be written in terms of the Fourier e tÞ as follows: image Eðm; e o2 E e ¼ ed0 ðtÞ; þ SðmÞ E ot2 e Ej t60 ¼ 0;
E
t 2 R;
where the matrix SðmÞ is given by 0 2 1 m2 þ m23 m1 m2 m1 m3 B C @ m1 m2 m21 þ m23 m2 m3 A: m1 m3
m2 m3
m21 þ m22
ð9Þ ð10Þ
ð11Þ
Let E be an arbitrary symmetric positive definite matrix and the symmetric matrix SðmÞ be defined by (11). Applying MATLAB we can find eigenvalues of SðmÞ. They are k1 ¼ 0, k2 ¼ jmj, k3 ¼ jmj ðjmj ¼ m21 þ m22 þ m23 Þ. Theorem 8.4.2 from [21] (see p. 366) says that for any symmetric matrix the eigenvalues are non-negative if and only if this matrix is positive semi-definite. Using this theorem we conclude that SðmÞ is symmetric positive semi-definite. According to the matrix theory (see, for example, [22], pp. 381–386) for any symmetric positive definite matrix E and any symmetric positive semi-definite matrix SðmÞ there exists a non-singular matrix TðmÞ and a diagonal matrix DðmÞ with non-negative elements which satisfy (19) and (20). The procedure of finding TðmÞ and DðmÞ consists of 1 several steps. The first step is to find E2 . At this step for the given positive definite matrix E ¼ diagðe11 ; e22 ; e33 Þ the 1 matrix E2 is given by 1 1 1 12 E ¼ diag pffiffiffiffiffiffi ; pffiffiffiffiffiffi ; pffiffiffiffiffiffi : e11 e22 e33 For the given positive definite non-diagonal matrix E we need to find a diagonal matrix of its eigenvalues, that is, we need to construct an orthogonal matrix P by eigenfunctions of E such that PT EP ¼ M diagðl1 ; l2 ; l3 Þ, where
1626
V.G. Yakhno, T.M. Yakhno / Computers and Structures 85 (2007) 1623–1634
PT is the transposed matrix to P and lk > 0, k ¼ 1; 2; 3 are 1 eigenvalues of E. The matrix M2 is defined by the rule 1 1 ffiffiffiffi ffi p p ffiffiffiffi ffi p ffiffiffiffi ffi M2 ¼ diagð l1 ; l2 ; l3 Þ. The matrix E2 is found by the 1 2
1 2
T
12
formula E ¼ PM P and E is the inverse matrix to 1 1 E2 . MATLAB code of finding E2 is listed below. INPUT: Eps; [EigVecEps, EigValEps] = eig(Eps); P = EigVecEps; PT = P’; M = EigValEps; Mh = sqrt(M); SqrEps = P * Mh * PT; InvSrtEps = inv(SqrEps); OUTPUT: InvSrtEps. Here Eps, SqrEps, InvSqrEps 1
1
are E, E2 , E2 , respectively. The second step of the procedure is to determine DðmÞ if 1 E2 is found and SðmÞ is given by (11). For this step we 1 1 consider the matrix E2 SðmÞE2 which is symmetric positive semi-definite. This matrix is congruent to a diagonal matrix of its eigenvalues and we construct an orthogonal 1 1 matrix Q by eigenfunctions of E2 SðmÞE2 such that h 1 i 1 QT E2 SðmÞE2 Q ¼ diagðd 1 ðmÞ; d 2 ðmÞ; d 3 ðmÞÞ; where QT ðmÞ is the transposed matrix to QðmÞ and d k ðmÞ P 0, k ¼ 1; 2; 3 are eigenvalues of 1
Here InvSrtEps, D, T, transT 1
are E2 , DðmÞ, TðmÞ, TT ðmÞ, respectively. Applying the above mentioned code for the case E ¼ diagðe11 ; e11 ; e33 Þ we found the following formulae for DðmÞ, TðmÞ: ! 2 jmj e11 m21 þ e11 m22 þ e33 m23 DðmÞ ¼ diag 0; ; ; e11 e11 e33 pffiffiffiffi 0 1 e m1 m 3 m m2 p1ffiffiffiffi 1 c p1ffiffiffiffi b e ðm332 þm 2Þ k e33 m3 e11 m1 11 1 2 B C pffiffiffiffi B C e33 m2 m3 m2 p1ffiffiffiffi b TðmÞ ¼ B p1ffiffiffiffi c k C; 2 2 e11 e11 ðm1 þm2 Þ A @ e33 m3 p1ffiffiffiffi c p1ffiffiffiffi k 0 e33 e33 where 1=2 e11 ðv21 þ v22 Þ þ e33 v23 k¼ ; e11 ðm21 þ m22 Þ 1=2 e11 m21 e11 m22 þ ; c¼ 1þ e33 m23 e33 m23 1=2 m2 ; jmj2 ¼ m21 þ m22 þ m23 : b ¼ 1 þ 22 m1 We note that formulae for matrices DðmÞ, TðmÞ become cumbersome even for the case of the diagonal matrix E with all different elements, and take around 10 printed pages for a general form of symmetric positive definite matrix E. 2.3. Exact formulae for the Fourier transform of electric and magnetic fields
1
E2 SðmÞE2 : At the third step of the procedure we define TðmÞ by the 1 formula TðmÞ ¼ E2 QðmÞ. Since the matrix SðmÞ is given symbolically then it becomes almost impossible to calculate the second and the third steps of our procedure ‘by hand’ even for the case of the diagonal matrix E. At this point we apply the methods of symbolic computations which allow us to automate linear algebra transformations. MATLAB code of these symbolic transformations for finding matrices DðmÞ, TðmÞ, TT ðmÞ is listed below: INPUT: InvSrtEps , S; A=simplify(InvSrtEps * S * InvSrtEps); [EigVecA, EigValA] = eig(A); D = simplify(EigValA); Q = simplify(EigVecA); Q(:,1) = Q(:,1)./sqrt(sum(Q(:,1).^2)); Q(:,2) = Q(:,2)./sqrt(sum(Q(:,2).^2)); Q(:,3) = Q(:,3)./sqrt(sum(Q(:,3).^2)); T = InvSrtEps * Q ; transT = T.’ ; OUTPUT: D, T, transT.
Using matrices TðmÞ, TT ðmÞ, DðmÞ relations (9) and (10) e tÞ as are written in the terms of Yðm; tÞ ¼ TT ðmÞ Eðm; follows: d2 Y þ DðmÞY ¼ TT ðmÞed0 ðtÞ; dt2 Yðm; tÞjt60 ¼ 0:
t 2 R;
ð21Þ ð22Þ
We note that e tÞ ¼ TðmÞYðm; tÞ Eðm;
ð23Þ
and the exact solution of the initial value problem (21) and (22) (see the ordinary differential equations technique, for example, [22]) is given by Yðm; tÞ ¼ columnðY 1 ðm; tÞ; Y 2 ðm; tÞ; Y 3 ðm; tÞÞ;
ð24Þ
where for t < 0 the function Y n ðm; tÞ vanishes and for t P 0 is defined by pffiffiffiffiffiffiffiffiffiffiffi d n ðmÞt ½TT ðmÞen ; ð25Þ Y n ðm; tÞ ¼ cos if d n ðmÞ > 0, and Y n ðm; tÞ ¼ ½TT ðmÞen ;
if d n ðmÞ ¼ 0;
ð26Þ
V.G. Yakhno, T.M. Yakhno / Computers and Structures 85 (2007) 1623–1634
where n ¼ 1; 2; 3; d n ðmÞ are elements of the matrix DðmÞ; ½TT ðmÞen is the nth component of the vector TT ðmÞe. An explicit formula for the Fourier transform of the electric field is given by (23), where Yðm; tÞ is defined by (24)–(26) and elements of matrices DðmÞ, TðmÞ, TT ðmÞ are found explicitly by symbolic computations described above. R t Using formulae (23)–(26) we find explicit formulae for e ðm; sÞds. They are the following E 0 j Z t e j ðm; sÞds ¼ hðtÞ½TðmÞXðm; tÞ ; E ð27Þ j 0
where j ¼ 1; 2; 3; hðtÞ is the Heaviside step function defined by hðtÞ ¼ 1 for t P 0 and hðtÞ ¼ 0 for t < 0; ½TðmÞXðm; tÞj is the jth component of the vector function TðmÞXðm; tÞ; Xðm; tÞ ¼ columnðX 1 ðm; tÞ; X 2 ðm; tÞ; X 3 ðm; tÞÞ is the vector function whose components are defined for t P 0 by pffiffiffiffiffiffiffiffiffiffiffi sinð d n ðmÞtÞ T X n ðm; tÞ ¼ ½T ðmÞen pffiffiffiffiffiffiffiffiffiffiffi ; ð28Þ d n ðmÞ
H 1 ðx; tÞ ¼
hðtÞ 3
ð2pÞ
Z
þ1
1
Z
þ1
1
Z
1627
þ1
1
fm2 ½TðmÞXðm; tÞ3
m3 ½TðmÞXðm; tÞ2 g sinðm1 x1 þ m2 x2 þ m3 x3 Þdm1 dm2 dm3 ;
ð37Þ H 2 ðx; tÞ ¼
hðtÞ 3
ð2pÞ
Z
þ1
1
Z
þ1
1
Z
þ1
1
fm3 ½TðmÞXðm; tÞ1
m1 ½TðmÞXðm; tÞ3 g sinðm1 x1 þ m2 x2 þ m3 x3 Þdm1 dm2 dm3 ;
ð38Þ H 3 ðx; tÞ ¼
hðtÞ 3
ð2pÞ
Z
þ1
1
Z
þ1
1
Z
þ1
1
fm1 ½TðmÞXðm; tÞ2
m2 ½TðmÞXðm; tÞ1 g sinðm1 x1 þ m2 x2 þ m3 x3 Þdm1 dm2 dm3 :
ð39Þ Formulae (34)–(39) give the exact solutions of IVPMF and IVPEF. They are the exact presentations of electric and magnetic fields from the impulse point current of the form jðx; tÞ ¼ edðx1 Þdðx2 Þdðx3 ÞdðtÞ, x ¼ ðx1 ; x2 ; x3 Þ 2 R3 , where e 2 R3 is a given unit vector, dð Þ is the Dirac delta function. 3. Visualization of electric and magnetic fields
if d n ðmÞ > 0, and X n ðm; tÞ ¼ t½TT ðmÞen
if d n ðmÞ ¼ 0:
ð29Þ
Substituting (28), (29) into (27) and then (27) into (12)– (14) we find the formula for the Fourier image of the magnetic field in the form e tÞ ¼ ð H e 1 ðm; tÞ; H e 2 ðm; tÞ; H e 3 ðm; tÞÞ; Hðm;
ð30Þ
where e 1 ðm; tÞ ¼ ihðtÞfm2 ½TðmÞXðm; tÞ m3 ½TðmÞXðm; tÞ g; ð31Þ H 3 2 e H 2 ðm; tÞ ¼ ihðtÞfm3 ½TðmÞXðm; tÞ1 m1 ½TðmÞXðm; tÞ3 g; ð32Þ e 3 ðm; tÞ ¼ ihðtÞfm1 ½TðmÞXðm; tÞ m2 ½TðmÞXðm; tÞ g: ð33Þ H 2
1
Here components of the vector function Xðm; tÞ ¼ ðX 1 ðm; tÞ; X 2 ðm; tÞ; X 3 ðm; tÞÞ are defined by (28) and (29) for m ¼ ðm1 ; m2 ; m3 Þ 2 R3 and t P 0. 2.4. Exact formulae for electric an magnetic fields e j ðm; tÞ, Ej ðx; tÞ, H j ðx; tÞ, We note that functions E j ¼ 1; 2; 3 have real values. Therefore, applying the inverse Fourier transform to (23) and (31)–(33) we obtain Eðx; tÞ ¼ ðE1 ðx; tÞ; E2 ðx; tÞ; E3 ðx; tÞÞ;
ð34Þ
Hðx; tÞ ¼ ðH 1 ðx; tÞ; H 2 ðx; tÞ; H 3 ðx; tÞÞ;
ð35Þ
where Ej ðx; tÞ ¼
hðtÞ
Z
þ1
Z
þ1
Z
þ1
½TðmÞYðm; tÞj 3 ð2pÞ 1 1 1 cosðm1 x1 þ m2 x2 þ m3 x3 Þ dm1 dm2 dm3 ;
j ¼ 1; 2; 3;
ð36Þ
3.1. General characteristics of the visualization For the visualization we consider five different electrically anisotropic non-dispersive homogeneous materials: a positive uniaxial crystal, a negative uniaxial crystal, a biaxial crystal, a crystal with monoclinic structure, and a crystal with a triclinic structure. For all these materials the source of electric and magnetic waves is the current density j taken in the form jðx; tÞ ¼ e1 dðx1 Þdðx2 Þdðx3 ÞdðtÞ; where dðx1 Þdðx2 Þdðx3 ÞdðtÞ is the Dirac delta function concentrated at the origin of the coordinates at the time t ¼ 0 in the direction e1 ¼ ð1; 0; 0Þ. Using the procedure of Section 2.2 the matrices TðmÞ, TT ðmÞ, DðmÞ were computed in the explicit form for every given E. Further, using formulae (34)–(39) we have found explicit formulae for the components of the electric and magnetic fields for each material. Using these formulae the images of E1 ð0; x2 ; x3 ; tÞ, E1 ðx1 ; 0; x3 ; tÞ, E1 ðx1 ; x2 ; 0; tÞ, H 1 ð0; x2 ; x3 ; tÞ, H 1 ðx1 ; 0; x3 ; tÞ, H 1 ðx1 ; x2 ; 0; tÞ were simulated for fixed time t. In this section, we show the following aspects of electric and magnetic fields behavior in different anisotropic materials: dynamics of the electric field propagation; magnitudes of electric and magnetic field components; the behavior of electric and magnetic fields (fronts, distribution of electric and magnetic energy). Images are organized in the following groups. Dynamics of the electric field propagation in a positive uniaxial crystal is shown in Figs. 1, 2 and 8(a). The magnitude of electric and magnetic field components as 3D plots for different types of crystals is shown in Figs. 1, 3–7. Figs. 8 and 9(a)–(c) illustrate shapes of electric and magnetic field
1628
V.G. Yakhno, T.M. Yakhno / Computers and Structures 85 (2007) 1623–1634
Fig. 1. 3D plot of E1 ð0; x2 ; x3 ; 1Þ; E ¼ diagð10; 10; 40Þ.
Fig. 3. 3D plot of H 1 ð0; x2 ; x3 ; 6Þ; E ¼ diagð10; 10; 120Þ.
Fig. 4. 3D plot of H 1 ðx1 ; x2 ; 0; 6Þ; E ¼ diagð120; 120; 10Þ.
Fig. 2. EF in a positive uniaxial crystal: E ¼ diagð10; 10; 40Þ. (a) 2D plot of E1 ð0; x2 ; x3 ; 1Þ; (b) 2D plot of E1 ð0; x2 ; x3 ; 4Þ]
Fig. 5. 3D plot of H 1 ðx1 ; x2 ; 0; 10:5Þ; E ¼ diagð15:15; 20; 30Þ.
V.G. Yakhno, T.M. Yakhno / Computers and Structures 85 (2007) 1623–1634
1629
Fig. 6. 3D plot of H 1 ð0; x2 ; x3 ; 3Þ of the crystal with the monoclinic structure; E ¼ Em .
Fig. 7. 3D plot of E1 ðx1 ; x2 ; 0; 0:09Þ of the crystal with the triclinic structure; E ¼ Et .
components for a positive uniaxial crystal. Figs. 9–14 illustrate the behavior (fronts, energy) of electric and magnetic fields in different types of anisotropic materials. 3.2. Description of input data and corresponding images Example 1 (Positive uniaxial crystal). The matrix E is given by E ¼ diagð10; 10; 40Þ. The result of the simulation E1 ð0; x2 ; x3 ; tÞ for t ¼ 1, t ¼ 4, t ¼ 6 is presented in Figs. 1, 2 and 8(a). Fig. 1 is a 3D plot of E1 ð0; x2 ; x3 ; tÞ for t ¼ 1. The horizontal axes are x2 and x3, respectively. The vertical axis is the magnitude of E1 ð0; x2 ; x3 ; 1Þ. The different colors correspond to different values of E1 ð0; x2 ; x3 ; 1Þ. Fig. 2(a) contains a screen shot of 2D level plot of the same surface E1 ð0; x2 ; x3 ; 1Þ, i.e. a view of the surface z ¼ E1 ð0; x2 ; x3 ; 1Þ presented in Fig. 1 from the top of z-axis. Figs. 2(b), 8(a)– (c) and 9(a)–(c) present screen shorts of 2D level plots
Fig. 8. EF in a positive uniaxial crystal: E ¼ diagð10; 10; 40Þ. (a) 2D plot of E1 ð0; x2 ; x3 ; 6Þ; (b) 2D plot of E1 ðx1 ; 0; x3 ; 6Þ; (c) 2D plot of E1 ðx1 ; x2 ; 0; 6Þ.
1630
V.G. Yakhno, T.M. Yakhno / Computers and Structures 85 (2007) 1623–1634
Fig. 9. MF in a positive uniaxial crystal: E ¼ diagð10; 10; 40Þ. (a) 2D plot of H 1 ð0; x2 ; x3 ; 6Þ; (b) 2D plot of H 1 ðx1 ; 0; x3 ; 6Þ; (c) 2D plot of H 1 ðx1 ; x2 ; 0; 6Þ.
Fig. 10. MF in a positive uniaxial crystal: E ¼ diagð10; 10; 120Þ. (a) 2D plot of H 1 ð0; x2 ; x3 ; 6Þ; (b) 2D plot of H 1 ðx1 ; 0; x3 ; 6Þ; (c) 2D plot of H 1 ðx1 ; x2 ; 0; 6Þ.
V.G. Yakhno, T.M. Yakhno / Computers and Structures 85 (2007) 1623–1634
Fig. 11. MF in a negative uniaxial crystal: E ¼ diagð120; 120; 10Þ. (a) 2D plot of H 1 ð0; x2 ; x3 ; 8Þ; (b) 2D plot of H 1 ðx1 ; 0; x3 ; 8Þ; (c) 2D plot of H 1 ðx1 ; x2 ; 0; 6Þ.
1631
Fig. 12. MF in a biaxial crystal: E ¼ diagð12:5; 20; 30Þ. (a) 2D plot of H 1 ð0; x2 ; x3 ; 10:5Þ; (b) 2D plot of H 1 ðx1 ; 0; x3 ; 10:5Þ; (c) 2D plot of H 1 ðx1 ; x2 ; 0; 10:5Þ.
1632
V.G. Yakhno, T.M. Yakhno / Computers and Structures 85 (2007) 1623–1634
Fig. 13. MF in a crystal with a monoclinic structure. (a) 2D plot of H 1 ð0; x2 ; x3 ; 3Þ; (b) 2D plot of H 1 ðx1 ; 0; x3 ; 3Þ; (c) 2D plot of H 1 ðx1 ; x2 ; 0; 5Þ.
Fig. 14. EF in a crystal with a triclinic structure. (a) 2D plot of E1 ð0; x2 ; x3 ; 0:09Þ; (b) 2D plot of E1 ðx1 ; 0; x3 ; 0:09Þ; (c) 2D plot of E1 ðx1 ; x2 ; 0; 0:09Þ.
V.G. Yakhno, T.M. Yakhno / Computers and Structures 85 (2007) 1623–1634
E1 ð0; x2 ; x3 ; 4Þ, E1 ð0; x2 ; x3 ; 6Þ, E1 ðx1 ; 0; x3 ; 6Þ, E1 ðx1 ; x2 ; 0; 6Þ, H 1 ð0; x2 ; x3 ; 6Þ, H 1 ðx1 ; 0; x3 ; 6Þ, H 1 ðx1 ; x2 ; 0; 6Þ, respectively. Example 2 (Positive uniaxial crystal). The matrix E is given by E ¼ diagð10; 10; 120Þ. Fig. 3 is a 3D plot of H 1 ð0; x2 ; x3 ; tÞ for t ¼ 6. The horizontal axes are x2 and x3, respectively. The vertical axis is the magnitude of H 1 ð0; x2 ; x3 ; 6Þ. Different colors correspond to different values of H 1 ð0; x2 ; x3 ; 6Þ. Fig. 10(a)–(c) contains screen shots of 2D level plots of H 1 ð0; x2 ; x3 ; 6Þ, H 1 ðx1 ; 0; x3 ; 6Þ, H 1 ðx1 ; x2 ; 0; 6Þ. Example 3 (Negative uniaxial crystal). The matrix E is given by E ¼ diagð120; 120; 10Þ. The result of the simulation H 1 ðx1 ; x2 ; 0; 6Þ is presented in Figs. 4 and 11(c). Fig. 4 is a 3D plot of H 1 ðx1 ; x2 ; 0; tÞ for t ¼ 6. The horizontal axes are x1 and x2, respectively. The vertical axis is the magnitude of H 1 ðx1 ; x2 ; 0; 6Þ. Fig. 11(c) contains a screen shot of 2D level plot of the same surface H 1 ðx1 ; x2 ; 0; 6Þ. Fig. 11(a) and (b) presents screen shots of 2D level plots of H 1 ð0; x2 ; x3 ; 8Þ and H 1 ðx1 ; 0; x3 ; 8Þ, respectively. Example 4 (Biaxial crystal). The matrix E is given by E ¼ diagð15:15; 20; 30Þ. Fig. 5 contains a 3D plot of H 1 ðx1 ; x2 ; 0; tÞ for t ¼ 10:5. The horizontal axes are x1 and x2, respectively. The vertical axis is the magnitude of H 1 ðx1 ; x2 ; 0; 10:5Þ. Fig. 12(c) presents a screen shot of 2D level plot of the same surface H 1 ðx1 ; x2 ; 0; 10:5Þ. Fig. 12(a) and (b) contains screen shots of 2D level plots of H 1 ð0; x2 ; x3 ; 10:5Þ and H 1 ðx1 ; 0; x3 ; 10:5Þ, respectively. Example 5 (Monoclinic structure). The matrix E is given by 0 1 17:1598 13:0178 0 B C Em ¼ @ 13:0178 23:6686 0 A: 0 0 44:4444 Fig. 6 contains a 3D plot of H 1 ð0; x2 ; x3 ; tÞ for t ¼ 3. The horizontal axes are x2 and x3, respectively. The vertical axis is the magnitude of H 1 ð0; x2 ; x3 ; 3Þ. Fig. 13(a) presents a screen shot of 2D level plot of the same surface presented in Fig. 6. Fig. 13(b) and (c) contains screen shots of 2D level plots of H 1 ðx1 ; 0; x3 ; 3Þ and H 1 ðx1 ; x2 ; 0; 5Þ, respectively. Example 6 (Triclinic structure). The matrix E is given by 0 1 0:00937 0:01776 0:01477 B C C Et ¼ B @ 0:01776 0:05327 0:03453 A: 0:01477 0:03453 0:08101 3D plot of E1 ðx1 ; x2 ; 0; tÞ for t ¼ 0:09 is shown in Fig. 7. A screen shot of 2D level plot of the same surface is shown in Fig. 14(c). The screen shots of 2D level plots of E1 ð0; x2 ; x3 ; 0:09Þ and E1 ðx1 ; 0; x3 ; 0:09Þ are shown in Fig. 14(a) and (b), respectively.
1633
3.3. Analysis of the visualization The different approaches for modelling of the behavior of electric and magnetic fields may be found in the literature [1–5]. The most popular is a ‘plane wave’ approach, in which the presentation of electric and magnetic fields is based on plane waves. Such presentation is a realistic approximation if, for example, a point pulse source is a great distance away from the observation point. But if a point pulse source is concentrated at a short distance from the observation point then electric and magnetic waves may take very peculiar forms. At the same time electric and magnetic fields are invisible. The simulations of these fields by modern computer tools especially for anisotropic materials allow us to see dependence between the material structure and the behavior of electric and magnetic fields. The method of our paper allows users to observe the exact electric and magnetic fields arising from a pulse point source. The different structures of anisotropic materials produce different responses of electric and magnetic fields inside these anisotropic materials. We can see various shapes of electric and magnetic wave propagations (different form of fronts and magnitude fluctuations) in presented figures.
4. Conclusions In the paper, we have described the new efficient method for modelling electric and magnetic waves inside different electrically anisotropic non-dispersive homogeneous materials. This method is based on constructing explicit formulae for the electric and magnetic fields by matrix symbolic computations and the inverse Fourier transform which is done numerically. We note that the simulation of the wave phenomena based on explicit formulae for exact solutions of partial differential equations is the best one according to accuracy. The robustness of our approach is based on the modern achievements of computational algebra which allow us to make symbolic transformations of cumbersome formulae and matrices. Our approach combines in a rational way analytical techniques and numerical computations. Using this approach we have generated images of electric and magnetic fields components from a pulse pointed electric source (pulse current concentrated at the origin of coordinates at time t ¼ 0). Analyzing these images we have found that electric and magnetic waves in different anisotropic crystals have different shapes of fronts and different distribution of the electric and magnetic energy. Using our method we have created a collection of images of electric and magnetic fields. This collection of images is organized now as a library which contains several hundreds of images classified according to the crystal system. This library can be used as a set of patterns for study and development of anisotropic materials. At the same time it can be used for verification of new numerical methods.
1634
V.G. Yakhno, T.M. Yakhno / Computers and Structures 85 (2007) 1623–1634
As a future application we plan to implement it as a Web application and use media streaming technology to have access to the images from the Internet. This may have considerable impact and potentially wide applications on engineering computational research and engineering education. References [1] Kong JA. Electromagnetic wave theory. New York: John Wiley and Sons; 1986. [2] Ramo S, Whinnery JR, Duzer T. Fields and waves in communication electronics. New York: John Wiley and Sons; 1994. [3] Eom HJ. Electromagnetic wave theory for boundary-value problems. Berlin: Springer; 2004. [4] Monk P. Finite element methods for Maxwell’s equations. Oxford: Clarendon Press; 2003. [5] Cohen GC. Higher-order numerical methods for transient wave equations. Berlin: Springer Verlag; 2002. [6] Dieulesaint E, Royer D. Elastic waves in solids. Chichester: John Wiley and Sons; 1980. [7] Lindell IV. Time-domain TE/TM decomposition of electromagnetic sources. IEEE Trans Anten Propag 1990;38(3):353–8. [8] Haba Z. Green functions and propagation of waves in strongly inhomogeneous media. J Phys A: Math Gen 2004;37:9295–302. [9] Wijnands F, Pendry JB, Garcia-Vidal FJ, Bell PM, Roberts PJ, Moreno LM. Green’s functions for Maxwell’s equations: application to spontaneous emission. Opt Quantum Electron 1997;29:199–216. [10] Li LW, Liu S, Leong MS, Yeo TS. Circular cylindrical waveguide filled with uniaxial anisotropic media-electromagnetic fields and
[11]
[12]
[13]
[14]
[15] [16]
[17]
[18] [19] [20] [21] [22]
dyadic Green’s functions. IEEE Tran Microwave Techn 2001;49(7):1361–4. Gottis PG, Kondylis GD. Properties of the dyadic Green’s function for unbounded anisotropic medium. IEEE Trans Anten Propag 1995;45:154–61. Ortner N, Wagner P. Fundamental matrices of homogeneous hyperbolic systems. Applications to crystal optics, elastodynamics, and piezoelectromagnetism. Z Angew Math Mech 2004;84(5):314–46. Yakhno VG. Constructing Green’s function for time-depending Maxwell system in anisotropic dielectrics. J Phys A: Math Gen 2005;38:2277–87. Berini P, Wu K. Modelling lossy anisotropic dielectric waveguides with the method of lines. IEEE Tran Microwave Techn 1996;44(5):749–59. Zienkiewicz OC, Taylor RL. Finite elements method 1. Oxford: Butterworth; 2000. Cohen GC, Heikkola E, Joly P, Neittaan MP. Mathematical and numerical aspects of waves propagation. Berlin: Springer Verlag; 2003. Yakhno VG, Yakhno TM, Kasap M. A novel approach for modelling and simulation of electromagnetic waves in anisotropic dielectrics. Int J Solids Struct 2006;43:6261–76. Pavlovic MN. Symbolic computation in structural engineering. Comput Struct 2003;81:2121–36. Pavlovic MN, Sapountzakis EJ. Computers and structures: nonnumerical applications. Comput Struct 1986;24:119–27. Beltzer AI. Engineering analysis via simbolic computation – a breakthrough. Appl Mech Rev 1990;43:119–27. Goldberg JL. Matrix theory with applications. New York: McGraw-Hill; 1992. Boyce WE, DiPrima RC. Elementary differential equations and boundary value problems. New York: John Wiley and Sons; 1992.
Computers and Structures 85 (2007) 1623–1634 www.elsevier.com/locate/compstruc
Modelling and simulation of electric and magnetic fields in homogeneous non-dispersive anisotropic materials V.G. Yakhno a
a,*
, T.M. Yakhno
b
Department of Mathematics, Dokuz Eylul University, Fen-Edebiyat Fak., Kaynaklar, Buca, Izmir, Turkey b Department of Computer Engineering, Dokuz Eylul University, Kaynaklar, Buca, Izmir, Turkey Received 5 December 2006; accepted 15 February 2007 Available online 24 April 2007
Abstract The paper describes a new efficient method for finding and simulation of electric and magnetic fields in homogeneous non-dispersive electrically anisotropic materials. Electromagnetic wave propagation in these materials is descried by the time-dependent Maxwell’s system which implies initial value problems (IVPs) for electric and magnetic fields. Our method essentially uses matrix symbolic calculations in MATLAB and allows us to compute explicit formulae of electric and magnetic fields. These formulae are used for generating images of electric and magnetic fields inside different materials. A collection of electric and magnetic field images is presented in the paper. Ó 2007 Elsevier Ltd. All rights reserved. Keywords: Electric and magnetic fields; Anisotropy; Wave propagation; Simulation; Matrix symbolic computations; MATLAB
1. Introduction 1.1. General background information Anisotropic materials attract a great deal of interest in modern engineering applications, in particular, microelectronics [1–3]. In general case a type of anisotropy specifies properties of materials, and electromagnetic waves propagation inside these materials essentially depends on the direction, generating electromagnetic fields of very peculiar forms. Modelling and simulation of waves inside these materials constitute an important interdisciplinary area of research with many cutting-edge scientific and technological applications. Fundamental models for electromagnetic wave propagation used in science and engineering are based on the Maxwell’s system [1–4]. In particular the mathematical models of electric and magnetic waves motions in electrically and magnetically anisotropic materials are very well known *
Corresponding author. Tel.: +90 0232 4128592; fax: +90 0232 4534188. E-mail address: [email protected] (V.G. Yakhno). 0045-7949/$ - see front matter Ó 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.compstruc.2007.02.019
(see, for example, [1–5]). These mathematical models are related to the time-dependent Maxwell’s system containing matrices of the dielectric permittivity E and magnetic permeability M. We note (see, for example, [1]) that a medium is electrically anisotropic if it is described by E ¼ ðeij Þ33 and M ¼ lI; and magnetically anisotropic if it is described by E ¼ eI and M ¼ ðlij Þ33 , where ðeij Þ33 , ðlij Þ33 are arbitrary matrices, e, l are positive constants, I is the identity matrix. A medium can be both electrically and magnetically anisotropic if E ¼ ðeij Þ33 and M ¼ ðlij Þ33 . 1.2. Permittivity, permeability and crystallographic structure of anisotropic materials A smallest block (three-dimensional array of atoms) of anisotropic materials is determined by repeated replication in three dimensions. Its symmetry tells how the constituent atoms are arranged in a regular repeating configuration. The structure of these three-dimensional ‘unit cell’ of atoms in anisotropic materials may have one of seven basic shapes: cubic, hexagonal, tetragonal, trigonal, orthorhombic, monoclinic and triclinic, see [6]. For linear homogeneous
1624
V.G. Yakhno, T.M. Yakhno / Computers and Structures 85 (2007) 1623–1634
non-dispersive electrically and magnetically anisotropic materials the electric and magnetic flux densities D and B have the following relations with the electric and magnetic fields E and H: D ¼ EE, B ¼ MH, where the permittivity E and the magnetic permeability M are symmetric positive definite matrices. In our paper, we will study homogeneous non-dispersive electrically anisotropic materials for which E is a symmetric positive definite matrix and M ¼ lI, l > 0 is a constant. The relations of the crystallographic structures of these materials and their dielectric permittivity tensors E are the following: E ¼ diagðe11 ; e11 ; e11 Þ for the cubic structure; E ¼ diagðe11 ; e11 ; e33 Þ for the hexagonal, tetragonal and trigonal structures; E ¼ ðeij Þ33 , eij ¼ eji , e13 ¼ e23 ¼ 0 for the monoclinic structure; and an arbitrary symmetric positive definite matrix E ¼ ðeij Þ33 for the triclinic structure of materials. 1.3. Basic equations and problems set-up The propagation of electromagnetic waves in homogeneous non-dispersive electrically anisotropic materials is described by the following time-dependent Maxwell’s system with a symmetric positive definite matrix of dielectric permittivity E ¼ ðeij Þ33 and the magnetic permeability M ¼ lI (see, for example, [1,2]). Let x ¼ ðx1 ; x2 ; x3 Þ be a space variable from R3, t be a time variable from R, then Maxwell’s system is given by the following relations: oE þ j; ot oH ; curlx E ¼ M ot
ð1Þ
curlx H ¼ E
ð2Þ
divx ðMHÞ ¼ 0;
ð3Þ
divx ðEEÞ ¼ q;
ð4Þ
where E ¼ ðE1 ; E2 ; E3 Þ, H ¼ ðH 1 ; H 2 ; H 3 Þ are electric and magnetic fields, Ek ¼ Ek ðx; tÞ; H k ¼ H k ðx; tÞ; k ¼ 1; 2; 3; j ¼ ðj1 ; j2 ; j3 Þ is the density of the electric current, jk ¼ jk ðx; tÞ; k ¼ 1; 2; 3; M is the magnetic permeability, E is the dielectric permittivity, q is the density of electric charges. The operators curlx and divx are defined by oE3 oE2 oE1 oE3 oE2 oE1 curlx E ¼ ; ; ; ox2 ox3 ox3 ox1 ox1 ox2 oE1 oE2 oE3 þ þ : divx E ¼ ox1 ox2 ox3 The conservation law of charges is given by (see, for example, [2]) oq þ divx j ¼ 0: ot
ð5Þ
We suppose that E ¼ 0; H ¼ 0; q ¼ 0
for t 6 0;
ð6Þ
this means there is no electric charges at the time t 6 0; electric and magnetic fields vanish for t 6 0.
Differentiating (1) with respect to t and using (2) and (6) we find curlx ðM1 curlx EÞ ¼ E Ejt60 ¼ 0;
o2 E oj þ ; ot2 ot
ð7Þ ð8Þ
where M1 is the inverse matrix to M. In this paper, we assume that M ¼ I, E ¼ ðeij Þ33 is a given symmetric positive definite matrix with constant elements and j is taken in the form jðx; tÞ ¼ edðx1 Þdðx2 Þdðx3 ÞdðtÞ; x ¼ ðx1 ; x2 ; x3 Þ 2 R3 , where e 2 R3 is a given unit vector, dð Þ is the Dirac delta function. This form of the current density is related to a directional impulse electric source concentrated at the origin of coordinates at the time t ¼ 0 in the direction e. As a result of it we conclude that the Maxwell’s system (1)–(4) with conditions (6) implies the following initial value problem for the electric field (IVPEF) and initial value problem for the magnetic field (IVPMF). IVPEF. Let E be a given symmetric positive definite matrix, M ¼ I. Find the vector function E satisfying (7) and (8). IVPMF. Let E be a solution of IVPEF. Find the vector function H satisfying (2) and the condition Hjt60 ¼ 0. 1.4. Information on the previous studies and an assessment of the present results For the electric field equation (7) different problems for isotropic and anisotropic materials have been studied and many applications have been made, see, for example, [1,2,4–17]. In particular, decomposition method for the case of isotropic materials (E is a diagonal matrix of the form E ¼ diagðe11 ; e11 ; e11 Þ) has been suggested in [7]. Analytic methods of Green’s functions constructions have been studied for the case of isotropic materials in [8,9]; for uniaxial anisotropic media ðE ¼ diagðe11 ; e11 ; e33 ÞÞ in [10,11]; for biaxial anisotropic crystals ðE ¼ diagðe11 ; e22 ; e33 ÞÞ in [12]; for arbitrary non-dispersive homogeneous anisotropic dielectrics (E ¼ ðeij Þ33 is a symmetric positive definite matrix) in [13]. Modelling lossy anisotropic dielectric waveguides with the method of lines has been made for inhomogeneous biaxial anisotropic media in [14]. Most of the studies and modelling electromagnetic waves had been made by numerical methods, in particular finite element method [4,5,15,16]. On the other hand nowadays computers can perform very complicated symbolic computations (in addition to numerical calculations) and this opens up new possibilities in modelling and simulation of wave propagation phenomena. Symbolic computations can be considered as a useful tool for analytical methods which can provide exact solutions of problems [18–20]. Unfortunately the exact solutions cannot be found for all complex equations and systems. But when the exact solution can be found it leads to the significant simplification of modelling and simula-
V.G. Yakhno, T.M. Yakhno / Computers and Structures 85 (2007) 1623–1634
tion and allows engineers to focus on the ideas rather than on overcoming calculational difficulties. In this paper, we present a new analytical method for constructing explicit formulae for electric and magnetic fields simultaneously inside different electrically anisotropic non-dispersive homogeneous materials. In general case these explicit formulae are very cumbersome and they have been computed using symbolic computations in MATLAB. Applying these explicit formulae simulations of magnetic waves are obtained. Images of magnetic field components from a pulse current concentrated in a fixed point are presented in the paper. The pulse pointed current was chosen because this source is convenient for modelling fronts of electric and magnetic waves. Presented images of field components give an opportunity to observe the traces of electric and magnetic wave fronts and the distribution of energy in materials with different anisotropic structure. The paper is organized as follows. In Section 2, we describe in details a procedure of finding exact formulae of electric and magnetic fields for arbitrary type of anisotropy. Section 3 contains images of the traces of electric and magnetic fields for different anisotropic materials classified in the introduction. The last section contains conclusions and directions of a future research. 2. Finding exact formulae for electric and magnetic fields
1625
IVPMF can be written in terms of Fourier images e tÞ, Hðm; e tÞ as follows: Eðm; e 1 ðm; tÞ oH e 3 ðm; tÞ m3 E e 2 ; ¼ i½m2 E ot e 2 ðm; tÞ oH e 1 ðm; tÞ m1 E e 3 ðm; tÞ; ¼ i½m3 E ot e 3 ðm; tÞ oH e 2 ðm; tÞ m2 E e 1 ðm; tÞ; ¼ i½m1 E ot fj ðm; tÞjt60 ¼ 0; j ¼ 1; 2; 3: H Integrating (12)–(14) under conditions (15) we find Z t e e 3 ðm; sÞ m3 E e 2 ðm; sÞds; H 1 ðm; tÞ ¼ i ½m2 E 0 Z t e 1 ðm; sÞ m1 E e 3 ðm; sÞds; e H 2 ðm; tÞ ¼ i ½m3 E 0 Z t e 2 ðm; sÞ m2 E e 1 ðm; sÞds: e 3 ðm; tÞ ¼ i ½m1 E H
ð12Þ ð13Þ ð14Þ ð15Þ
ð16Þ ð17Þ ð18Þ
0
One of the mainRpoints of constructing explicit formulae e ðm; sÞds; j ¼ 1; 2; 3 is a procedure of e j ðm; tÞ and t E for E 0 j finding explicitly a non-singular matrix TðmÞ and diagonal matrix DðmÞ with non-negative elements such that TT ðmÞETðmÞ ¼ I;
ð19Þ
T
T ðmÞSðmÞTðmÞ ¼ DðmÞ;
ð20Þ T
2.1. IVPEF and IVPMF in terms of Fourier transform
where I is the identity matrix, T ðmÞ is the transposed matrix to TðmÞ.
e tÞ, Hðm; e tÞ be the Fourier transform images of Let Eðm; the electric and magnetic fields Eðx; tÞ and Hðx; tÞ with respect to x ¼ ðx1 ; x2 ; x3 Þ 2 R3 , i.e.
2.2. Procedure of finding matrices TðmÞ, DðmÞ and its implementation in MATLAB
e tÞ ¼ ð E e 1 ðm; tÞ; E e 2 ðm; tÞ; E e 3 ðm; tÞÞ; Eðm; e tÞ ¼ ð H e 2 ðm; tÞ; H e 3 ðm; tÞÞ; e 1 ðm; tÞ; H Hðm; e j ðm; tÞ ¼ Fx ½Ej ðm; tÞ; E
e j ðm; tÞ ¼ Fx ½H j ðm; tÞ; H
where j ¼ 1; 2; 3; m ¼ ðm1 ; m2 ; m3 Þ 2 R3 ; Z þ1 Z þ1 Z þ1 Fx ½Ej ðm; tÞ ¼ Ej ðx; tÞeimx dx1 dx2 dx3 ; 1
1
xm ¼ x1 m1 þ x2 m2 þ x3 m3 ;
1
i2 ¼ 1:
IVPEF (7) and (8) can be written in terms of the Fourier e tÞ as follows: image Eðm; e o2 E e ¼ ed0 ðtÞ; þ SðmÞ E ot2 e Ej t60 ¼ 0;
E
t 2 R;
where the matrix SðmÞ is given by 0 2 1 m2 þ m23 m1 m2 m1 m3 B C @ m1 m2 m21 þ m23 m2 m3 A: m1 m3
m2 m3
m21 þ m22
ð9Þ ð10Þ
ð11Þ
Let E be an arbitrary symmetric positive definite matrix and the symmetric matrix SðmÞ be defined by (11). Applying MATLAB we can find eigenvalues of SðmÞ. They are k1 ¼ 0, k2 ¼ jmj, k3 ¼ jmj ðjmj ¼ m21 þ m22 þ m23 Þ. Theorem 8.4.2 from [21] (see p. 366) says that for any symmetric matrix the eigenvalues are non-negative if and only if this matrix is positive semi-definite. Using this theorem we conclude that SðmÞ is symmetric positive semi-definite. According to the matrix theory (see, for example, [22], pp. 381–386) for any symmetric positive definite matrix E and any symmetric positive semi-definite matrix SðmÞ there exists a non-singular matrix TðmÞ and a diagonal matrix DðmÞ with non-negative elements which satisfy (19) and (20). The procedure of finding TðmÞ and DðmÞ consists of 1 several steps. The first step is to find E2 . At this step for the given positive definite matrix E ¼ diagðe11 ; e22 ; e33 Þ the 1 matrix E2 is given by 1 1 1 12 E ¼ diag pffiffiffiffiffiffi ; pffiffiffiffiffiffi ; pffiffiffiffiffiffi : e11 e22 e33 For the given positive definite non-diagonal matrix E we need to find a diagonal matrix of its eigenvalues, that is, we need to construct an orthogonal matrix P by eigenfunctions of E such that PT EP ¼ M diagðl1 ; l2 ; l3 Þ, where
1626
V.G. Yakhno, T.M. Yakhno / Computers and Structures 85 (2007) 1623–1634
PT is the transposed matrix to P and lk > 0, k ¼ 1; 2; 3 are 1 eigenvalues of E. The matrix M2 is defined by the rule 1 1 ffiffiffiffi ffi p p ffiffiffiffi ffi p ffiffiffiffi ffi M2 ¼ diagð l1 ; l2 ; l3 Þ. The matrix E2 is found by the 1 2
1 2
T
12
formula E ¼ PM P and E is the inverse matrix to 1 1 E2 . MATLAB code of finding E2 is listed below. INPUT: Eps; [EigVecEps, EigValEps] = eig(Eps); P = EigVecEps; PT = P’; M = EigValEps; Mh = sqrt(M); SqrEps = P * Mh * PT; InvSrtEps = inv(SqrEps); OUTPUT: InvSrtEps. Here Eps, SqrEps, InvSqrEps 1
1
are E, E2 , E2 , respectively. The second step of the procedure is to determine DðmÞ if 1 E2 is found and SðmÞ is given by (11). For this step we 1 1 consider the matrix E2 SðmÞE2 which is symmetric positive semi-definite. This matrix is congruent to a diagonal matrix of its eigenvalues and we construct an orthogonal 1 1 matrix Q by eigenfunctions of E2 SðmÞE2 such that h 1 i 1 QT E2 SðmÞE2 Q ¼ diagðd 1 ðmÞ; d 2 ðmÞ; d 3 ðmÞÞ; where QT ðmÞ is the transposed matrix to QðmÞ and d k ðmÞ P 0, k ¼ 1; 2; 3 are eigenvalues of 1
Here InvSrtEps, D, T, transT 1
are E2 , DðmÞ, TðmÞ, TT ðmÞ, respectively. Applying the above mentioned code for the case E ¼ diagðe11 ; e11 ; e33 Þ we found the following formulae for DðmÞ, TðmÞ: ! 2 jmj e11 m21 þ e11 m22 þ e33 m23 DðmÞ ¼ diag 0; ; ; e11 e11 e33 pffiffiffiffi 0 1 e m1 m 3 m m2 p1ffiffiffiffi 1 c p1ffiffiffiffi b e ðm332 þm 2Þ k e33 m3 e11 m1 11 1 2 B C pffiffiffiffi B C e33 m2 m3 m2 p1ffiffiffiffi b TðmÞ ¼ B p1ffiffiffiffi c k C; 2 2 e11 e11 ðm1 þm2 Þ A @ e33 m3 p1ffiffiffiffi c p1ffiffiffiffi k 0 e33 e33 where 1=2 e11 ðv21 þ v22 Þ þ e33 v23 k¼ ; e11 ðm21 þ m22 Þ 1=2 e11 m21 e11 m22 þ ; c¼ 1þ e33 m23 e33 m23 1=2 m2 ; jmj2 ¼ m21 þ m22 þ m23 : b ¼ 1 þ 22 m1 We note that formulae for matrices DðmÞ, TðmÞ become cumbersome even for the case of the diagonal matrix E with all different elements, and take around 10 printed pages for a general form of symmetric positive definite matrix E. 2.3. Exact formulae for the Fourier transform of electric and magnetic fields
1
E2 SðmÞE2 : At the third step of the procedure we define TðmÞ by the 1 formula TðmÞ ¼ E2 QðmÞ. Since the matrix SðmÞ is given symbolically then it becomes almost impossible to calculate the second and the third steps of our procedure ‘by hand’ even for the case of the diagonal matrix E. At this point we apply the methods of symbolic computations which allow us to automate linear algebra transformations. MATLAB code of these symbolic transformations for finding matrices DðmÞ, TðmÞ, TT ðmÞ is listed below: INPUT: InvSrtEps , S; A=simplify(InvSrtEps * S * InvSrtEps); [EigVecA, EigValA] = eig(A); D = simplify(EigValA); Q = simplify(EigVecA); Q(:,1) = Q(:,1)./sqrt(sum(Q(:,1).^2)); Q(:,2) = Q(:,2)./sqrt(sum(Q(:,2).^2)); Q(:,3) = Q(:,3)./sqrt(sum(Q(:,3).^2)); T = InvSrtEps * Q ; transT = T.’ ; OUTPUT: D, T, transT.
Using matrices TðmÞ, TT ðmÞ, DðmÞ relations (9) and (10) e tÞ as are written in the terms of Yðm; tÞ ¼ TT ðmÞ Eðm; follows: d2 Y þ DðmÞY ¼ TT ðmÞed0 ðtÞ; dt2 Yðm; tÞjt60 ¼ 0:
t 2 R;
ð21Þ ð22Þ
We note that e tÞ ¼ TðmÞYðm; tÞ Eðm;
ð23Þ
and the exact solution of the initial value problem (21) and (22) (see the ordinary differential equations technique, for example, [22]) is given by Yðm; tÞ ¼ columnðY 1 ðm; tÞ; Y 2 ðm; tÞ; Y 3 ðm; tÞÞ;
ð24Þ
where for t < 0 the function Y n ðm; tÞ vanishes and for t P 0 is defined by pffiffiffiffiffiffiffiffiffiffiffi d n ðmÞt ½TT ðmÞen ; ð25Þ Y n ðm; tÞ ¼ cos if d n ðmÞ > 0, and Y n ðm; tÞ ¼ ½TT ðmÞen ;
if d n ðmÞ ¼ 0;
ð26Þ
V.G. Yakhno, T.M. Yakhno / Computers and Structures 85 (2007) 1623–1634
where n ¼ 1; 2; 3; d n ðmÞ are elements of the matrix DðmÞ; ½TT ðmÞen is the nth component of the vector TT ðmÞe. An explicit formula for the Fourier transform of the electric field is given by (23), where Yðm; tÞ is defined by (24)–(26) and elements of matrices DðmÞ, TðmÞ, TT ðmÞ are found explicitly by symbolic computations described above. R t Using formulae (23)–(26) we find explicit formulae for e ðm; sÞds. They are the following E 0 j Z t e j ðm; sÞds ¼ hðtÞ½TðmÞXðm; tÞ ; E ð27Þ j 0
where j ¼ 1; 2; 3; hðtÞ is the Heaviside step function defined by hðtÞ ¼ 1 for t P 0 and hðtÞ ¼ 0 for t < 0; ½TðmÞXðm; tÞj is the jth component of the vector function TðmÞXðm; tÞ; Xðm; tÞ ¼ columnðX 1 ðm; tÞ; X 2 ðm; tÞ; X 3 ðm; tÞÞ is the vector function whose components are defined for t P 0 by pffiffiffiffiffiffiffiffiffiffiffi sinð d n ðmÞtÞ T X n ðm; tÞ ¼ ½T ðmÞen pffiffiffiffiffiffiffiffiffiffiffi ; ð28Þ d n ðmÞ
H 1 ðx; tÞ ¼
hðtÞ 3
ð2pÞ
Z
þ1
1
Z
þ1
1
Z
1627
þ1
1
fm2 ½TðmÞXðm; tÞ3
m3 ½TðmÞXðm; tÞ2 g sinðm1 x1 þ m2 x2 þ m3 x3 Þdm1 dm2 dm3 ;
ð37Þ H 2 ðx; tÞ ¼
hðtÞ 3
ð2pÞ
Z
þ1
1
Z
þ1
1
Z
þ1
1
fm3 ½TðmÞXðm; tÞ1
m1 ½TðmÞXðm; tÞ3 g sinðm1 x1 þ m2 x2 þ m3 x3 Þdm1 dm2 dm3 ;
ð38Þ H 3 ðx; tÞ ¼
hðtÞ 3
ð2pÞ
Z
þ1
1
Z
þ1
1
Z
þ1
1
fm1 ½TðmÞXðm; tÞ2
m2 ½TðmÞXðm; tÞ1 g sinðm1 x1 þ m2 x2 þ m3 x3 Þdm1 dm2 dm3 :
ð39Þ Formulae (34)–(39) give the exact solutions of IVPMF and IVPEF. They are the exact presentations of electric and magnetic fields from the impulse point current of the form jðx; tÞ ¼ edðx1 Þdðx2 Þdðx3 ÞdðtÞ, x ¼ ðx1 ; x2 ; x3 Þ 2 R3 , where e 2 R3 is a given unit vector, dð Þ is the Dirac delta function. 3. Visualization of electric and magnetic fields
if d n ðmÞ > 0, and X n ðm; tÞ ¼ t½TT ðmÞen
if d n ðmÞ ¼ 0:
ð29Þ
Substituting (28), (29) into (27) and then (27) into (12)– (14) we find the formula for the Fourier image of the magnetic field in the form e tÞ ¼ ð H e 1 ðm; tÞ; H e 2 ðm; tÞ; H e 3 ðm; tÞÞ; Hðm;
ð30Þ
where e 1 ðm; tÞ ¼ ihðtÞfm2 ½TðmÞXðm; tÞ m3 ½TðmÞXðm; tÞ g; ð31Þ H 3 2 e H 2 ðm; tÞ ¼ ihðtÞfm3 ½TðmÞXðm; tÞ1 m1 ½TðmÞXðm; tÞ3 g; ð32Þ e 3 ðm; tÞ ¼ ihðtÞfm1 ½TðmÞXðm; tÞ m2 ½TðmÞXðm; tÞ g: ð33Þ H 2
1
Here components of the vector function Xðm; tÞ ¼ ðX 1 ðm; tÞ; X 2 ðm; tÞ; X 3 ðm; tÞÞ are defined by (28) and (29) for m ¼ ðm1 ; m2 ; m3 Þ 2 R3 and t P 0. 2.4. Exact formulae for electric an magnetic fields e j ðm; tÞ, Ej ðx; tÞ, H j ðx; tÞ, We note that functions E j ¼ 1; 2; 3 have real values. Therefore, applying the inverse Fourier transform to (23) and (31)–(33) we obtain Eðx; tÞ ¼ ðE1 ðx; tÞ; E2 ðx; tÞ; E3 ðx; tÞÞ;
ð34Þ
Hðx; tÞ ¼ ðH 1 ðx; tÞ; H 2 ðx; tÞ; H 3 ðx; tÞÞ;
ð35Þ
where Ej ðx; tÞ ¼
hðtÞ
Z
þ1
Z
þ1
Z
þ1
½TðmÞYðm; tÞj 3 ð2pÞ 1 1 1 cosðm1 x1 þ m2 x2 þ m3 x3 Þ dm1 dm2 dm3 ;
j ¼ 1; 2; 3;
ð36Þ
3.1. General characteristics of the visualization For the visualization we consider five different electrically anisotropic non-dispersive homogeneous materials: a positive uniaxial crystal, a negative uniaxial crystal, a biaxial crystal, a crystal with monoclinic structure, and a crystal with a triclinic structure. For all these materials the source of electric and magnetic waves is the current density j taken in the form jðx; tÞ ¼ e1 dðx1 Þdðx2 Þdðx3 ÞdðtÞ; where dðx1 Þdðx2 Þdðx3 ÞdðtÞ is the Dirac delta function concentrated at the origin of the coordinates at the time t ¼ 0 in the direction e1 ¼ ð1; 0; 0Þ. Using the procedure of Section 2.2 the matrices TðmÞ, TT ðmÞ, DðmÞ were computed in the explicit form for every given E. Further, using formulae (34)–(39) we have found explicit formulae for the components of the electric and magnetic fields for each material. Using these formulae the images of E1 ð0; x2 ; x3 ; tÞ, E1 ðx1 ; 0; x3 ; tÞ, E1 ðx1 ; x2 ; 0; tÞ, H 1 ð0; x2 ; x3 ; tÞ, H 1 ðx1 ; 0; x3 ; tÞ, H 1 ðx1 ; x2 ; 0; tÞ were simulated for fixed time t. In this section, we show the following aspects of electric and magnetic fields behavior in different anisotropic materials: dynamics of the electric field propagation; magnitudes of electric and magnetic field components; the behavior of electric and magnetic fields (fronts, distribution of electric and magnetic energy). Images are organized in the following groups. Dynamics of the electric field propagation in a positive uniaxial crystal is shown in Figs. 1, 2 and 8(a). The magnitude of electric and magnetic field components as 3D plots for different types of crystals is shown in Figs. 1, 3–7. Figs. 8 and 9(a)–(c) illustrate shapes of electric and magnetic field
1628
V.G. Yakhno, T.M. Yakhno / Computers and Structures 85 (2007) 1623–1634
Fig. 1. 3D plot of E1 ð0; x2 ; x3 ; 1Þ; E ¼ diagð10; 10; 40Þ.
Fig. 3. 3D plot of H 1 ð0; x2 ; x3 ; 6Þ; E ¼ diagð10; 10; 120Þ.
Fig. 4. 3D plot of H 1 ðx1 ; x2 ; 0; 6Þ; E ¼ diagð120; 120; 10Þ.
Fig. 2. EF in a positive uniaxial crystal: E ¼ diagð10; 10; 40Þ. (a) 2D plot of E1 ð0; x2 ; x3 ; 1Þ; (b) 2D plot of E1 ð0; x2 ; x3 ; 4Þ]
Fig. 5. 3D plot of H 1 ðx1 ; x2 ; 0; 10:5Þ; E ¼ diagð15:15; 20; 30Þ.
V.G. Yakhno, T.M. Yakhno / Computers and Structures 85 (2007) 1623–1634
1629
Fig. 6. 3D plot of H 1 ð0; x2 ; x3 ; 3Þ of the crystal with the monoclinic structure; E ¼ Em .
Fig. 7. 3D plot of E1 ðx1 ; x2 ; 0; 0:09Þ of the crystal with the triclinic structure; E ¼ Et .
components for a positive uniaxial crystal. Figs. 9–14 illustrate the behavior (fronts, energy) of electric and magnetic fields in different types of anisotropic materials. 3.2. Description of input data and corresponding images Example 1 (Positive uniaxial crystal). The matrix E is given by E ¼ diagð10; 10; 40Þ. The result of the simulation E1 ð0; x2 ; x3 ; tÞ for t ¼ 1, t ¼ 4, t ¼ 6 is presented in Figs. 1, 2 and 8(a). Fig. 1 is a 3D plot of E1 ð0; x2 ; x3 ; tÞ for t ¼ 1. The horizontal axes are x2 and x3, respectively. The vertical axis is the magnitude of E1 ð0; x2 ; x3 ; 1Þ. The different colors correspond to different values of E1 ð0; x2 ; x3 ; 1Þ. Fig. 2(a) contains a screen shot of 2D level plot of the same surface E1 ð0; x2 ; x3 ; 1Þ, i.e. a view of the surface z ¼ E1 ð0; x2 ; x3 ; 1Þ presented in Fig. 1 from the top of z-axis. Figs. 2(b), 8(a)– (c) and 9(a)–(c) present screen shorts of 2D level plots
Fig. 8. EF in a positive uniaxial crystal: E ¼ diagð10; 10; 40Þ. (a) 2D plot of E1 ð0; x2 ; x3 ; 6Þ; (b) 2D plot of E1 ðx1 ; 0; x3 ; 6Þ; (c) 2D plot of E1 ðx1 ; x2 ; 0; 6Þ.
1630
V.G. Yakhno, T.M. Yakhno / Computers and Structures 85 (2007) 1623–1634
Fig. 9. MF in a positive uniaxial crystal: E ¼ diagð10; 10; 40Þ. (a) 2D plot of H 1 ð0; x2 ; x3 ; 6Þ; (b) 2D plot of H 1 ðx1 ; 0; x3 ; 6Þ; (c) 2D plot of H 1 ðx1 ; x2 ; 0; 6Þ.
Fig. 10. MF in a positive uniaxial crystal: E ¼ diagð10; 10; 120Þ. (a) 2D plot of H 1 ð0; x2 ; x3 ; 6Þ; (b) 2D plot of H 1 ðx1 ; 0; x3 ; 6Þ; (c) 2D plot of H 1 ðx1 ; x2 ; 0; 6Þ.
V.G. Yakhno, T.M. Yakhno / Computers and Structures 85 (2007) 1623–1634
Fig. 11. MF in a negative uniaxial crystal: E ¼ diagð120; 120; 10Þ. (a) 2D plot of H 1 ð0; x2 ; x3 ; 8Þ; (b) 2D plot of H 1 ðx1 ; 0; x3 ; 8Þ; (c) 2D plot of H 1 ðx1 ; x2 ; 0; 6Þ.
1631
Fig. 12. MF in a biaxial crystal: E ¼ diagð12:5; 20; 30Þ. (a) 2D plot of H 1 ð0; x2 ; x3 ; 10:5Þ; (b) 2D plot of H 1 ðx1 ; 0; x3 ; 10:5Þ; (c) 2D plot of H 1 ðx1 ; x2 ; 0; 10:5Þ.
1632
V.G. Yakhno, T.M. Yakhno / Computers and Structures 85 (2007) 1623–1634
Fig. 13. MF in a crystal with a monoclinic structure. (a) 2D plot of H 1 ð0; x2 ; x3 ; 3Þ; (b) 2D plot of H 1 ðx1 ; 0; x3 ; 3Þ; (c) 2D plot of H 1 ðx1 ; x2 ; 0; 5Þ.
Fig. 14. EF in a crystal with a triclinic structure. (a) 2D plot of E1 ð0; x2 ; x3 ; 0:09Þ; (b) 2D plot of E1 ðx1 ; 0; x3 ; 0:09Þ; (c) 2D plot of E1 ðx1 ; x2 ; 0; 0:09Þ.
V.G. Yakhno, T.M. Yakhno / Computers and Structures 85 (2007) 1623–1634
E1 ð0; x2 ; x3 ; 4Þ, E1 ð0; x2 ; x3 ; 6Þ, E1 ðx1 ; 0; x3 ; 6Þ, E1 ðx1 ; x2 ; 0; 6Þ, H 1 ð0; x2 ; x3 ; 6Þ, H 1 ðx1 ; 0; x3 ; 6Þ, H 1 ðx1 ; x2 ; 0; 6Þ, respectively. Example 2 (Positive uniaxial crystal). The matrix E is given by E ¼ diagð10; 10; 120Þ. Fig. 3 is a 3D plot of H 1 ð0; x2 ; x3 ; tÞ for t ¼ 6. The horizontal axes are x2 and x3, respectively. The vertical axis is the magnitude of H 1 ð0; x2 ; x3 ; 6Þ. Different colors correspond to different values of H 1 ð0; x2 ; x3 ; 6Þ. Fig. 10(a)–(c) contains screen shots of 2D level plots of H 1 ð0; x2 ; x3 ; 6Þ, H 1 ðx1 ; 0; x3 ; 6Þ, H 1 ðx1 ; x2 ; 0; 6Þ. Example 3 (Negative uniaxial crystal). The matrix E is given by E ¼ diagð120; 120; 10Þ. The result of the simulation H 1 ðx1 ; x2 ; 0; 6Þ is presented in Figs. 4 and 11(c). Fig. 4 is a 3D plot of H 1 ðx1 ; x2 ; 0; tÞ for t ¼ 6. The horizontal axes are x1 and x2, respectively. The vertical axis is the magnitude of H 1 ðx1 ; x2 ; 0; 6Þ. Fig. 11(c) contains a screen shot of 2D level plot of the same surface H 1 ðx1 ; x2 ; 0; 6Þ. Fig. 11(a) and (b) presents screen shots of 2D level plots of H 1 ð0; x2 ; x3 ; 8Þ and H 1 ðx1 ; 0; x3 ; 8Þ, respectively. Example 4 (Biaxial crystal). The matrix E is given by E ¼ diagð15:15; 20; 30Þ. Fig. 5 contains a 3D plot of H 1 ðx1 ; x2 ; 0; tÞ for t ¼ 10:5. The horizontal axes are x1 and x2, respectively. The vertical axis is the magnitude of H 1 ðx1 ; x2 ; 0; 10:5Þ. Fig. 12(c) presents a screen shot of 2D level plot of the same surface H 1 ðx1 ; x2 ; 0; 10:5Þ. Fig. 12(a) and (b) contains screen shots of 2D level plots of H 1 ð0; x2 ; x3 ; 10:5Þ and H 1 ðx1 ; 0; x3 ; 10:5Þ, respectively. Example 5 (Monoclinic structure). The matrix E is given by 0 1 17:1598 13:0178 0 B C Em ¼ @ 13:0178 23:6686 0 A: 0 0 44:4444 Fig. 6 contains a 3D plot of H 1 ð0; x2 ; x3 ; tÞ for t ¼ 3. The horizontal axes are x2 and x3, respectively. The vertical axis is the magnitude of H 1 ð0; x2 ; x3 ; 3Þ. Fig. 13(a) presents a screen shot of 2D level plot of the same surface presented in Fig. 6. Fig. 13(b) and (c) contains screen shots of 2D level plots of H 1 ðx1 ; 0; x3 ; 3Þ and H 1 ðx1 ; x2 ; 0; 5Þ, respectively. Example 6 (Triclinic structure). The matrix E is given by 0 1 0:00937 0:01776 0:01477 B C C Et ¼ B @ 0:01776 0:05327 0:03453 A: 0:01477 0:03453 0:08101 3D plot of E1 ðx1 ; x2 ; 0; tÞ for t ¼ 0:09 is shown in Fig. 7. A screen shot of 2D level plot of the same surface is shown in Fig. 14(c). The screen shots of 2D level plots of E1 ð0; x2 ; x3 ; 0:09Þ and E1 ðx1 ; 0; x3 ; 0:09Þ are shown in Fig. 14(a) and (b), respectively.
1633
3.3. Analysis of the visualization The different approaches for modelling of the behavior of electric and magnetic fields may be found in the literature [1–5]. The most popular is a ‘plane wave’ approach, in which the presentation of electric and magnetic fields is based on plane waves. Such presentation is a realistic approximation if, for example, a point pulse source is a great distance away from the observation point. But if a point pulse source is concentrated at a short distance from the observation point then electric and magnetic waves may take very peculiar forms. At the same time electric and magnetic fields are invisible. The simulations of these fields by modern computer tools especially for anisotropic materials allow us to see dependence between the material structure and the behavior of electric and magnetic fields. The method of our paper allows users to observe the exact electric and magnetic fields arising from a pulse point source. The different structures of anisotropic materials produce different responses of electric and magnetic fields inside these anisotropic materials. We can see various shapes of electric and magnetic wave propagations (different form of fronts and magnitude fluctuations) in presented figures.
4. Conclusions In the paper, we have described the new efficient method for modelling electric and magnetic waves inside different electrically anisotropic non-dispersive homogeneous materials. This method is based on constructing explicit formulae for the electric and magnetic fields by matrix symbolic computations and the inverse Fourier transform which is done numerically. We note that the simulation of the wave phenomena based on explicit formulae for exact solutions of partial differential equations is the best one according to accuracy. The robustness of our approach is based on the modern achievements of computational algebra which allow us to make symbolic transformations of cumbersome formulae and matrices. Our approach combines in a rational way analytical techniques and numerical computations. Using this approach we have generated images of electric and magnetic fields components from a pulse pointed electric source (pulse current concentrated at the origin of coordinates at time t ¼ 0). Analyzing these images we have found that electric and magnetic waves in different anisotropic crystals have different shapes of fronts and different distribution of the electric and magnetic energy. Using our method we have created a collection of images of electric and magnetic fields. This collection of images is organized now as a library which contains several hundreds of images classified according to the crystal system. This library can be used as a set of patterns for study and development of anisotropic materials. At the same time it can be used for verification of new numerical methods.
1634
V.G. Yakhno, T.M. Yakhno / Computers and Structures 85 (2007) 1623–1634
As a future application we plan to implement it as a Web application and use media streaming technology to have access to the images from the Internet. This may have considerable impact and potentially wide applications on engineering computational research and engineering education. References [1] Kong JA. Electromagnetic wave theory. New York: John Wiley and Sons; 1986. [2] Ramo S, Whinnery JR, Duzer T. Fields and waves in communication electronics. New York: John Wiley and Sons; 1994. [3] Eom HJ. Electromagnetic wave theory for boundary-value problems. Berlin: Springer; 2004. [4] Monk P. Finite element methods for Maxwell’s equations. Oxford: Clarendon Press; 2003. [5] Cohen GC. Higher-order numerical methods for transient wave equations. Berlin: Springer Verlag; 2002. [6] Dieulesaint E, Royer D. Elastic waves in solids. Chichester: John Wiley and Sons; 1980. [7] Lindell IV. Time-domain TE/TM decomposition of electromagnetic sources. IEEE Trans Anten Propag 1990;38(3):353–8. [8] Haba Z. Green functions and propagation of waves in strongly inhomogeneous media. J Phys A: Math Gen 2004;37:9295–302. [9] Wijnands F, Pendry JB, Garcia-Vidal FJ, Bell PM, Roberts PJ, Moreno LM. Green’s functions for Maxwell’s equations: application to spontaneous emission. Opt Quantum Electron 1997;29:199–216. [10] Li LW, Liu S, Leong MS, Yeo TS. Circular cylindrical waveguide filled with uniaxial anisotropic media-electromagnetic fields and
[11]
[12]
[13]
[14]
[15] [16]
[17]
[18] [19] [20] [21] [22]
dyadic Green’s functions. IEEE Tran Microwave Techn 2001;49(7):1361–4. Gottis PG, Kondylis GD. Properties of the dyadic Green’s function for unbounded anisotropic medium. IEEE Trans Anten Propag 1995;45:154–61. Ortner N, Wagner P. Fundamental matrices of homogeneous hyperbolic systems. Applications to crystal optics, elastodynamics, and piezoelectromagnetism. Z Angew Math Mech 2004;84(5):314–46. Yakhno VG. Constructing Green’s function for time-depending Maxwell system in anisotropic dielectrics. J Phys A: Math Gen 2005;38:2277–87. Berini P, Wu K. Modelling lossy anisotropic dielectric waveguides with the method of lines. IEEE Tran Microwave Techn 1996;44(5):749–59. Zienkiewicz OC, Taylor RL. Finite elements method 1. Oxford: Butterworth; 2000. Cohen GC, Heikkola E, Joly P, Neittaan MP. Mathematical and numerical aspects of waves propagation. Berlin: Springer Verlag; 2003. Yakhno VG, Yakhno TM, Kasap M. A novel approach for modelling and simulation of electromagnetic waves in anisotropic dielectrics. Int J Solids Struct 2006;43:6261–76. Pavlovic MN. Symbolic computation in structural engineering. Comput Struct 2003;81:2121–36. Pavlovic MN, Sapountzakis EJ. Computers and structures: nonnumerical applications. Comput Struct 1986;24:119–27. Beltzer AI. Engineering analysis via simbolic computation – a breakthrough. Appl Mech Rev 1990;43:119–27. Goldberg JL. Matrix theory with applications. New York: McGraw-Hill; 1992. Boyce WE, DiPrima RC. Elementary differential equations and boundary value problems. New York: John Wiley and Sons; 1992.