Z. Angew. Math. Phys. 62 (2011), 1117–1129 c 2011 Springer Basel AG  0044-2275/11/061117-13 published online June 16, 2011 DOI 10.1007/s00033-011-0135-2

Zeitschrift f¨ ur angewandte Mathematik und Physik ZAMP

On the wave propagation in isotropic fractal media Hady Joumaa and Martin Ostoja-Starzewski

Abstract. In this paper, we explore the wave propagation phenomenon in three-dimensional (3D) isotropic fractal media through analytical and computational means. We present the governing scalar wave equation, perform its eigenvalue decomposition, and discuss its corresponding modal solutions. The homogenization through which this fractal wave equation is derived makes its mathematical analysis and consequently the formulation of exact solutions possible if treated in the spherical coordinate system. From the computational perspective, we consider the finite element method and derive the corresponding weak formulation which can be implemented in the numerical scheme. The Newmark time-marching method solves the resulting elastodynamic system and captures the transient response. Two solvers capable of handling problems of arbitrary initial and boundary conditions for arbitrary domains are developed. They are validated in space and time, with particular problems considered on spherical shell domains. The first solver is elementary; it handles problems of purely radial dependence, effectively, 1D. However, the second one deals with general advanced 3D problems of arbitrary spatial dependence. Mathematics Subject Classification (2000). 74H45 Vibration. Keywords. Fractal media · Wave propagation · Spherical shell.

1. Introduction Fractal media are abundant in nature (e.g., rocks, tree leaves, island distributions, clouds) and in living bodies (e.g., brains, cardiovascular systems) [1–3]. Thus, the understanding of their mechanics in the sense of a field theory capable of handling initial-boundary value problems (IBVP) is sorely needed. Fractal bodies are characterized by highly irregular (nonsmooth) topologies; as a result, the Eucledian geometry and consequently the classical continuum mechanics theory fail to provide an expressive model able to simulate the mechanical behavior under scrutiny. The failure of the continuum theory opens the gate for new approaches, the aim of which is the exploration of the laws that govern the general dynamical behavior of fractal media. In fact, fractal mechanics is at a formative stage [4–9] in contradistinction to the extensively explored, advanced, and ubiquitously applied conventional continuum mechanics. A practical and well-known problem that demonstrates the necessity as well as the effectiveness of fractal mechanics is the wave propagation in porous media. Indeed a porous material, cannot be truthfully regarded as a continuum because of the successive void regions spread throughout its volume. The crucial parameter on which the fractal nature of a given body is to be considered or disregarded in the investigation, is the “length scale”. In other words, fractal features must be incorporated in the mechanical analysis only when it is investigated within a scale range limited from above by the overall macroscopic length of the body and from below by the microscopic scale of the smallest fractal feature [10]. In addition, the total mass M of a given porous body of characteristic length L does not scale in power law with the Eucledian dimension of the space that contains that body. It nevertheless scales with the Hausdorff dimension D, which is noninteger and smaller than the Eucledian dimension, i.e., M ∼ kLD where 2 < D < 3 for a fractal body that is contained in a three-dimensional Eucledian space [11].

1118

H. Joumaa and M. Ostoja-Starzewski

ZAMP

In this research, we report a study of wave propagation in 3D isotropic fractal media. This study, although not directly connected with a real application, promotes the general understanding of wave motion and elastodynamics in fractal bodies, the issue that allows scientists and engineers to incorporate this analysis into some mechanics problems where fractal effects are influential and their negligence is erroneous. The 3D fractal wave equation is derived following the approach of homogenizing fractional integrals into the Eucledian space through dimensional regularization [4,5]. The full derivation of the wave equation is provided in the appendix A. This equation, expressed in terms of conventional (nonfractional) derivatives, is given by (Γ being the Euler Gamma function)  2     2D−3 Γ D ∂2p 2−D 2  3  |R| (3 − D) R · ∇p + |R|2 ∇2 p = c (1) 2 ∂t Γ 2 While linear in principle, Eq. (1) involves additional complex terms reflecting the fractal nature of the medium, hindering the chances of obtaining analytical solutions to general problems. In addition to the wave celerity c, the fractal dimension D (a medium-dependent constant) plays a significant role in wave propagation phenomena. We can readily verify that this equation reduces to the classical 3D wave equation in the classical continuum case. Due to the isotropic nature of homogenization, the position vector R, appears in the equation. For this reason, obtaining analytical solutions to this equation in the Cartesian or cylindrical coordinate system is impossible. However, the spherical coordinate system allows one to decompose the equation into three independent modal equations, each of which can be solved analytically. The objective of this work is to first introduce the modal decomposition by explaining the mathematical steps leading to the exact solutions to the reduced modal equations. As an application, we consider a spherically symmetric eigenvalue problem involved with radial wave propagation on a spherical shell, which we treat analytically and asymptotically. Then, we construct the first Finite Element Method (FEM) solver which handles problems of radial dependence only, thus numerically demonstrating this special eigenvalue problem. This solver is validated in space and time with the exact solution. Finally, the general solver which handles full 3D problems is presented, and all applied procedures leading to the final elastodynamic model are explained.

2. Modal decomposition The aim of the modal decomposition is to be able to formulate analytical solutions to general problems using the technique of modal superposition. For the continuum case, modal decomposition in Cartesian (x, y, z), cylindrical (r, θ, z), and spherical (r, θ, φ) coordinates generates three decoupled ordinary differential equations [12, Chap. 9]. In the fractal domain, the first two frames of references are futile; they fail to produce independent equations. Fortunately, the spherical system produces the desirable decoupled equations. We start our modal analysis by separating time and space variables and introducing the frequency ω, thus p (R, t) = P (R) eiωt

(2)

Substituting (2) into (1), we obtain what we will refer to as the fractal Helmholtz equation expressed in compact and expanded forms, respectively, ∇2D P + k 2 P = 0  2     2D−3 Γ D  3  2 |R|2−D c (3 − D) R · ∇P + |R|2 ∇2 P + k 2 P = 0 Γ 2

(3a) (3b)

where the wavenumber k = ω/c. Expressing the dependent variable P in (2) as a product of three independent modal functions, P (r, θ, φ) ≡ F (r) · G (θ) · H (φ), and substituting back into (3b), we obtain

Vol. 62 (2011)

On the wave propagation in isotropic fractal media

1119

three independent Sturm–Liouville equations for a 3D eigenvalue problem expressed as φ : H  + m2 H = 0  m2 G + l (l + 1) − θ : G + G=0 tan θ sin2 θ

 2 r : r2 F  + (5 − D) rF  + (λk) r2D−4 − l (l + 1) F = 0

(4a) (4b) (4c)

Recall that, in spherical coordinates, we have r = |R|

(5a)

∂P R · ∇P = r ∂r



∂2P ∂ ∂P ∂ 1 1 ∂P 1 ∇2 P = 2 r2 + 2 sin θ + 2 2 r ∂r ∂r r sin θ ∂θ ∂θ r sin θ ∂φ2

(5b) (5c)

The constant λ showing in (4c) will always be pre-multiplied by the wave number k, it is defined as follows   Γ 32 3−D 2 (6) λ = D Γ 2 and l and m are two independent integers characterizing the modes, so-called mode numbers. Let us physically interpret these equations. First, we realize that the Hausdorff dimension D appears in the radial equation only. This is because Tarasov’s approach to homogenization relies solely on the radial distance from the origin to the point of interest, so that the azimuthal and transcendental modes are free of fractal effects [4]. Consequently, the azimuthal and transcendental equations are exactly identical to those of the continuum problem, see [12, p. 251]. In addition, upon setting D = 3, the radial mode equation reduces to the spherical Bessel equation, which offers an additional verification of this analysis. Having realized that the radial mode equation is similar to the spherical Bessel equation triggers an important hint in forming its solution as will be later explained. The analytical solutions to the modal equations are now discussed. Concerning (4a), it has the wellknown harmonic solution H (φ) = h1 sin (mφ) + h2 cos (mφ)

(7)

The transcendental equation (4b) can be transformed into a Legendre differential equation if we apply the transformation η = cos θ. The transformed equation is    d2 G m2 dG 1 − η2 + l (l + 1) − − 2η G=0 (8) dη 2 dη 1 − η2 |m|

One solution to this equation is the Legendre function Pl (η). Note that, if we set m = 0, we obtain the Legendre polynomial, denoted as Pl (η). The other homogeneous solution to (8) is not analytic at η = ±1. This renders it physically meaningless, for this problem, therefore the Legendre function is the solely considered solution for the spherical wave problem. The radial mode equation (4c) is the most challenging to solve—no closed form solution is obvious at a glance. We propose a series form solution for this equation so as to obtain two linearly independent homogeneous solutions given as

D−4 λkrD−2 (1) 2 Jν Fν (r, k, D) = r (9a) D−2

D−4 λkrD−2 Fν(2) (r, k, D) = r 2 Yν (9b) D−2 F (r) = f1 Fν(1) (r) + f2 Fν(2) (r)

(9c)

1120

H. Joumaa and M. Ostoja-Starzewski

ZAMP

v v v

v v v

(a)

(b)

Fig. 1. The fractal radial harmonic function of first kind, shown for different D, k = 1. a l = 0 b l = 1

v v v

v v v

(a)

(b)

Fig. 2. The fractal radial harmonic function of second kind, shown for different D, k = 1. a l = 0 b l = 1

√ 4l(l+1)+(D−4)2 where the constant ν = > 0, while Jν and Yν are, respectively, Bessel functions of the 2(D−2) first and second kind of order ν. Effectively, in (9a) and (9b) we have introduced the fractal radial harmonic functions of the first and second kind. Figures 1 and 2 show some plots of these two functions for several values of D and l. The unnormalized sinc function is reproduced for the case D = 3, l = 0. Having solved all three Eq. (4a–4c), we can now express the solution to wave propagation problems in terms of these functions while satisfying the necessary initial and boundary conditions as shown in the next section.

Vol. 62 (2011)

On the wave propagation in isotropic fractal media

(a)

1121

(b)

Fig. 3. The spherical shell domain in two views along with the spatial meshing along the radial direction for the one-dimensional reduced problem. a Outer view of the shell b Cross-section view of the shell

3. Fractal IBVP—case study We demonstrated in the modal decoupling analysis that fractal effects are observable in the radial harmonics only. Therefore, 3D problems can be reduced to 1D (in the radial coordinate) while preserving the fractal effects that we intend to explore. As an application, a Sturm–Liouville (eigenvalue) problem involved in the wave propagation inside a spherical shell is considered. The analytical solution is presented in terms of the fractal radial harmonic functions described before. In addition, the WKB asymptotic method is applied, providing an approximate but meaningful solution to the problem. 3.1. Exact solution Consider the spherical shell centered at the origin with inner radius Rin = 1 and outer radius Rout = 2, as shown in Fig. 3. Suppressing the dependence of θ and φ on p, the resulting wave equation for this reduced problem simplifies to  2 ∂ 2 p  c 2 4−2D ∂p 2∂ p +r = r (5 − D) r (10) ∂t2 λ ∂r ∂r2 The corresponding reduced fractal Helmholtz equation is r2 P  + (5 − D) rP  + λ2 k 2 r2D−4 P = 0

(11)

The conditions are of Dirichlet type at both ends of the domain. For simplicity, we choose  boundary  P Rin = P (Rout ) = 0 The general homogeneous solution for P (r) is obtained in terms of harmonic functions as Pn (r) = an Fν(1) (r, kn , D) + bn Fν(2) (r, kn , D) 4−D 2(D−2)

(12)

The nontrivial values (eigenvalues) of kn , along with the corresponding relation between where ν = constants an and bn , are determined through the application of boundary conditions which result in the

1122

H. Joumaa and M. Ostoja-Starzewski

ZAMP

Fig. 4. First three orthonormal modal functions. Note the remarkable match between the exact and asymptotic solutions

following relations Fν(1) (1, λkn , D) Fν(2) (2, λkn , D) − Fν(1) (2, λkn , D) Fν(2) (1, λkn , D) = 0 an = bn

(2) Fν − (1) Fν

(1, λkn , D)

(13a) (13b)

(1, λkn , D)

Equation (13a) is of transcendental nature; it admits infinite solutions in k. A general root-finding algorithm for nonlinear algebraic equations (bisection) is applied to solve for the eigenvalues. Rewriting (11) in the Sturm–Liouville form [13, p. 29] as  d 2 5−D dP (14) r + (λk) rD−1 P = 0 dr dr leads to the orthonormality condition expressed as 2

Pn Pm rD−1 dr = δmn

(15)

1

The first three orthonormal modal functions are shown in Fig. 4.

3.2. Asymptotic modal analysis The determination of the higher frequency modes can be performed effectively through an asymptotic approach. The WKB method is a well-known asymptotic procedure that predicts higher modes in eigenvalue problems particularly when the formulation of the exact solution represents a mathematical challenge. The WKB method is thoroughly explained in reference [13, Chap. 10]. We will briefly present the steps followed to reach the asymptotic solution. The WKB method expresses the solution P (r) as an infinite series of the form P (r) = eiλkσ(r)

∞  Ψj (r) j=0

j

(λk)

(16)

Vol. 62 (2011)

On the wave propagation in isotropic fractal media

1123

Table 1. Exact and asymptotic (WKB) solution of the eigenvalue problem for the first five modes, D = 2.5 n 1 2 3 4 5

λkn 3.8377 7.6076 11.392 15.181 18.971

˜n λk 3.7922 7.5844 11.376 15.169 18.961

Rel. error 0.0118 0.0030 0.0013 0.0007 0.0004

an bn

−20.4082 0.4551 −1.1084 1.6348 −0.2749

An Bn

−3.6202 0.5981 −0.9548 1.8624 −0.2271

Rel. error 0.8226 0.3142 0.1386 0.1222 0.1739

Substituting this form into (11), we obtain a series of decreasing power of k. Satisfying each and every order of k, the following relations are obtained   rD−2 O k 2 : −r6−2D σ 2 + 1 = 0 ⇒ σ (r) = ± (17a) D−2 1 O (k) : 2rΨ0 + 2Ψ0 = 0 ⇒ Ψ0 (r) = (17b) r  5−D   Ψj d (rΨj+1 ) i d r = (17c) dr 2r dr Equation (17c) provides the recursive relation to solve for Ψj+1 knowing Ψj . Truncating the series to one term only (Ψ0 ), the resulting nth mode corresponding to the nth wavenumber kn is approximated by: 

λkn rD−2 λkn rD−2 1  Pn (r) = An cos + Bn sin (18) r D−2 D−2 Applying the boundary conditions on the asymptotic solution Pn , a transcendental equation for the wavenumber kn and the (An , Bn ) relation are obtained as follows nπ (D − 2) k˜n = λ (2D−2 − 1)

nπ An = − tan Bn 2D−2 − 1

(19a) (19b)

If the asymptotic modes were to be orthonormalized by the same condition of (15), the resulting form is obtained 

rD−2 − 1 2(D − 2) 1  sin nπ D−2 Pn (r) = (20) 2D−2 − 1 r 2 −1 The exact and asymptotic solutions for the case D = 2.5 are listed in Table 1. We performed the analysis of the first five modes. Note the remarkable matching between both solutions at higher modes. Even though the WKB method targets higher modes, the lower modes are still well approximated. The plot of the modal functions in Fig. 4 illustrates the significance of the asymptotic approximation. Interestingly, the components of the WKB homogeneous solution shown in (18) are nothing but the asymptotic expansion of the harmonic functions F (1) and F (2) at large k. References [13, p. 572] and [12, App. A4] provide a detailed asymptotic expansion to the Bessel functions leading to a first order approximation for F (1) and F (2) given as   √  1 1 (1) (1) F1.5 ≈ F1.5 (r, k, D = 2.5) = − cos 2λk r as k → +∞ (21a) r πλk   √  1 1 (2) (2) sin 2λk r as k → +∞ F1.5 ≈ F1.5 (r, k, D = 2.5) = − (21b) r πλk As such, we can perform the asymptotic analysis by applying the boundary conditions on Pn (r) = An F (1) + Bn F (2) and reproduce the same results of the WKB work. The inconvenience of this approach

1124

H. Joumaa and M. Ostoja-Starzewski

ZAMP

lies in its reliance on the exact solution, nevertheless, it remains significant through its strong approximation of the higher modes while avoiding the differential equation solving dictated by the WKB and the implementation of the root finding algorithm required by the exact analysis. In conclusion, the two asymptotic approaches are consistent, both providing a meaningful prediction to higher modes.

4. FEM for reduced problem In this section, we introduce the FEM as a numerical alternative to solve the reduced problem of wave propagation. We discuss the FEM scheme construction through the weak formulation. We next solve transient problems based on modal excitation. The numerical results are verified and convergence is confirmed through error analysis.

4.1. FEM construction The reduced wave equation is presented in (10). We seek an FEM transient solution to the reduced problem whose domain and Dirichlet boundary conditions are already presented in the modal analysis. On the temporal side, the system is excited with a particular mode, i.e., p0 = Pn , a modal function. As such, the expected response is a single frequency harmonic, corresponding to the excited mode. The weak formulation is obtained following standard FEM procedures explained in [14,15]. We introduce an admissible function pˆ, continuous throughout the domain of interest and satisfying the boundary conditions. Consider the strong (differential) form of (10), then multiplying both sides by pˆ and integrating over the domain Ω, the weak (integral) form is generated. 2  2   2 ∂ p λ 5−2D ∂p 6−2D ∂ p p ˆ dr + p ˆ dr = (5 − D) r r pˆ dr (22) c ∂t2 ∂r ∂r2 Ω

Ω

Ω

Integrating by parts and performing additional calculus, we obtain 2  2 Rout   λ ∂ p ∂p ∂p 6−2D 5−2D ∂p pˆ dr + r pˆ dr = (D − 1) r pˆ − 2 c ∂t ∂r ∂r ∂r Rin Ω

Ω

Ω

∂ pˆ 6−2D r dr ∂r

(23)

Upon applying the boundary conditions, the middle  term on the right-hand side vanishes. The integration is then performed on an elemental basis Ωe , where e Ωe = Ω. The weak form is thus rewritten as 2  Ne  Ne  Ne    λ ∂2p ∂p ∂ pˆ 6−2D 5−2D ∂p pˆ dr + r pˆ dr + (1 − D) r dr = 0 (24) 2 c ∂t ∂r ∂r ∂r e=1 e=1 e=1 Ωe

Ωe

Ωe

For simplicity, the domain is meshed with a 2-node linear element. The integration can be performed either analytically or by implementing the Gaussian quadrature rule [14, Chap. 3]. The former method is chosen because the latter one fails to produce the exact value due to the noninteger powers present in the kernels. The elemental integration of the first term of (24) generates the elemental mass matrix Me , the second term generates the fractal elastic matrix Le , while the third term generates the regular elastic matrix Ke . The nomenclature of these matrices will be better appreciated when discussing the 3D formulation in the upcoming section. Summing all the resulting elemental matrices, the discrete formulation is produced N  Ne Ne e    Me p¨ + Le + Ke p = 0 ⇒ M¨ p + (L + K) p = 0 (25) e=1

e=1

e=1

Vol. 62 (2011)

On the wave propagation in isotropic fractal media

1125

The obtained model is constitutively similar to its continuum elastodynamic counterpart which balances inertial forces with elastic ones. As such, the numerical treatments of the continuous system remain applicable in the fractal field. In this problem, the mass matrix M and the regular elastic matrix K are symmetric positive definite (SPD), while the fractal elastic matrix L is not. Thus, the total stiffness tensor L+K is in general not SPD. These remarks can be proved by observing the elemental forms of these matrices. Having formulated the discrete elastodyanmic model, eigenvalue (frequency) analysis in addition to transient analysis can be conducted. The latter topic is discussed next.

4.2. Time-marching solution We explain the time-marching scheme we applied to solve (25) and obtain the transient response. The numerical analysis literature ([14, Chap. 9],[15, Chap. 9]) describes a variety of methods to solve this linear hyperbolic system, whereby we note the Newmark method. It can be implemented either implicitly or explicitly depending on its parameters adjustments. For our analysis, we have chosen the trapezoidal (implicit) case of Newmark’s method; it is unconditionally stable, does not generate any dissipative errors, and has a good accuracy. In conclusion, the numerical algorithm of the trapezoidal method can be implemented for fractal problems; the matching between numerical and analytical solutions is a clear proof of the method’s validity. The first three modal excitations are reported in Fig. 5. The error analysis in time and space is conducted; convergence plots are shown in Fig. 6 where the L2 norm of the error is plotted with respect to the element size and time step.

5. FEM for 3D problems In this section, we extend the application of the FEM to general 3D problems. The formulation is conceptually similar to that of the reduced problem but mathematically more complex. We adhere to the same problem solved on the spherical shell to ease the perception of the obtained solutions. Needless to say, the designed solver can handle general problems of arbitrary domain and boundary conditions.

5.1. FEM construction The weak form for this problem is obtained by considering an admissible function pˆ, multiplying both sides of the strong form shown in (1) by this function, and then integrating both sides of the equation over the domain of interest, Ω. The corresponding weak form is given as 2  2   ∂ p λ 4−2D p ˆ dΩ = (3 − D) p ˆ |R| R · ∇pdΩ + pˆ|R|6−2D ∇2 pdΩ (26) c ∂t2 Ω

Ω

Ω

Calculus rearrangements are made for the kernel of the third integral to engage the admissible function into the gradient operator. Applying the Green–Gauss theorem to this integral, we obtain   ˆ · ∇p (27a) pˆr6−2D ∇2 p = ∇ · pˆr6−2D ∇p − r6−2D ∇ˆ p · ∇p − (6 − 2D) pˆr4−2D R     pˆr6−2D ∇2 p dΩ = pˆr6−2D ∇p dS − r6−2D ∇ˆ p · ∇p dΩ + (2D − 6) pˆr4−2D R · ∇p dΩ (27b) Ω

S

Ω

Ω

1126

H. Joumaa and M. Ostoja-Starzewski

ZAMP

(b)

(a)

(c) Fig. 5. The system numerical (1D FEM) and exact solutions for the first three modal excitations. a First mode excitation b Second mode excitation c Third mode excitation

Applying the Dirichlet boundary conditions on the inner and outer spherical boundaries, the surface integral vanishes. Incorporating (27b) into (26), the weak form becomes 2  2   λ ∂ p 4−2D p ˆ dΩ = (D − 3) p ˆ r R · ∇p dΩ − r6−2D ∇ˆ p · ∇p dΩ (28) c ∂t2 Ω

Ω

Ω

The integration is conducted on the elemental level. We thus have 2  Ne  Ne  Ne    λ ∂2p 4−2D p ˆ dΩ = (D − 3) p ˆ r R · ∇p dΩ − r6−2D ∇ˆ p · ∇p dΩe e e 2 c ∂t e=1 e=1 e=1 Ωe

Ωe

(29)

Ωe

The above equation is similar to (24). The three governing matrices (M, L, and K) introduced in the reduced problem are easily identified. For the continuum case (D = 3), the middle term, which produces the fractal elastic matrix L vanishes. It is clear by now why L is named as such: it disappears

Vol. 62 (2011)

On the wave propagation in isotropic fractal media

1127

(a)

(b)

Fig. 6. Convergence plots in time and space. The order of accuracy is shown for each mode. a Time convergence b Space convergence

(a)

(b)

Fig. 7. The system numerical (3D FEM) and exact solutions for the first two modal excitations. a First mode excitation. b Second mode excitation

from the formulation when fractal effects are absent. On the other hand, the symmetric regular elastic matrix K regains its conventional form and reflects regular elastic effects in the case of continuous medium. The mass matrix M is indifferent to fractal effects; this is a direct consequence of the homogeneity (constant density) property of the fractal medium. The procedure to derive the elemental matrices requires integrating over 3D elements, 4-node tetrahedral being the simplest choice. The integrals cannot be evaluated exactly; the use of a quadrature method becomes a must. A useful mathematical algorithm for the derivation of quadrature points and corresponding weights over tetrahedral elements is provided in [16]. Assembling the governing matrices, the discrete elastodynamic form presented in (25) is obtained but with a much larger number of degrees of freedom. The Newmark method is applied to solve for the transient solution. The first two modal excitations are simulated and the response is shown in Fig. 7.

1128

H. Joumaa and M. Ostoja-Starzewski

ZAMP

6. Conclusions The applicability of analytical and numerical methods in understanding the wave propagation phenomena in isotopic fractal media is demonstrated in this research work. The remarkable consistency observed in all the results validates all these diverse approaches to solving the problem. In pursuing an analytical solution of the reduced 3D problem, we introduce the fractal radial harmonic functions of the first and second kind, a generalization of Bessel functions of the first and second kind. The elastodynamic analysis embedding modal decomposition, solving eigenvalue system, and predicting transient response is achieved for this noncontinuum type of material. The development of a general 3D solver, implementing the robust FEM, and applying the stable and convergent trapezoidal time-stepping scheme, allows solution of IBVPs in media with fractal geometries. Advanced fractal materials of greater complexity, encountered in real life applications, can now be modeled and investigated by further building upon the work presented here.

Acknowledgments This research was made possible by the support of Sandia-DTRA (grant HDTRA1-08-10-BRCWMD) and the NSF (grant CMMI-1030940).

Appendix A: derivation of the wave equation The wave equation in isotropic fractal media is derived by considering small perturbations to the governing balance laws describing the fractal medium dynamics. We strongly encourage the reader to refer to reference [4], where all the dynamic laws for the fractional continuous medium model are developed. Two relevant mathematical operators are first defined. The fractional spatial derivative (gradient) operator ∇D k is given as   2D−3 Γ D ∂A D  3  2 |R|3−D ∇k [A] = (30) ∂xk Γ 2 and the fractional material derivative is defined as d ∂A + uk ∇D [A] = k A dt D ∂t For the three-dimensional isotropic fractal media, the mass balance law is presented as follows d ρ = −ρ∇D k uk dt D while the momentum balance law is given as d 1 uk = fk − ∇D p dt D ρ k

(31)

(32)

(33)

Assuming isentropic (reversible and adiabatic) processes for small amplitude oscillation (small disturbance variables: uk , p∗ , and ρ∗ ) about an equilibrium point (¯ p,¯ ρ), we get, for a fluid of compressibility κ p − p¯ = κ

ρ − ρ¯ ρ∗ ⇒ p∗ = κ ρ¯ ρ¯

Introducing the gradient operator into (33), we have  d D D ∇k ρ uk = −∇D k ∇k p dt D

(34)

(35)

Vol. 62 (2011)

On the wave propagation in isotropic fractal media

1129

Neglecting the higher order terms (convective ones and ρ∗ uk product), we simplify the above equation to  D ∂uk D ∗ (36) ρ¯∇k = −∇D k ∇k p ∂t Doing the same for the mass balance in (32) we obtain, ∂ρ∗ = −¯ ρ∇D k uk ∂t Differentiating the above equation with respect to time,

∂uk ∂ 2 ρ∗ ∂  D  D ∇k uk = −¯ ρ∇k = −¯ ρ ∂t2 ∂t ∂t

(37)

(38)

Combining (35) and (38), we finally obtain ∂ 2 ρ∗ ∂ 2 p∗ κ D ∗ = ∇D = ∇D ∇D p∗ k ∇k p ⇒ 2 2 ∂t ∂t ρ¯ k k If we assign the isotropic wave celerity c2 = the fractal wave equation presented in (1)

κ ρ¯ ,

(39)

∗ and expand the ∇D k operator twice on p , we reproduce

References 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15. 16.

Mandelbrot, B.B.: The Fractal Geometry of Nature. W.H. Freeman & Co, San Francisco (1982) Barnsely, M.: Fractals Everywhere. Morgan Kaufmann, Los Altos (1993) Le Mehaute, A.: Fractal Geometries Theory and Applications. CRC Press, Boca Raton (1991) Tarasov, V.E.: Fractional hydrodynamic equations for fractal media. Ann. Phys. 318/2, 286–307 (2005) Tarasov, V.E.: Wave equation for fractal solid string. Mod. Phys. Lett. B 19(15), 721–728 (2005) Ostoja-Starzewski, M.: Towards thermomechanics of fractal media. ZAMP 58(6), 1085–1096 (2007) Ostoja-Starzewski, M.: Extremum and variational principles for elastic and inelastic media with fractal geometries. Acta Mech. 205, 161–170 (2009) Ostoja-Starzewski, M.: On turbulence in fractal porous media. ZAMP 59(6), 1111–1117 (2008) Ostoja-Starzewski, M., Li, J.: Fractal materials, beams and fracture mechanics. ZAMP 60, 1–12 (2009) Falconer, K.: Fractal Geometry: Mathematical Foundations and Applications. Wiley, England (2003) Hastings, H.M., Sugihara, G.: Fractals: A User’s Guide for the Natural Sciences. Oxford Science Publications, Oxford (1993) Kinsler, L.E., Frey, A.R.: Fundamentals of Acoustics. Wiley, New York (2000) Bender, C.M., Orszag, Steven A.: Advanced Mathematical Methods for Scientists and Engineers. Springer, Berlin (1999) Hughes, T.J.R.: The Finite Element Method. Dover Publications, New York (2000) Bathe, K.J.: Finite Element Procedures. Prentice Hall, NJ (1982) Rathold, H.T., Venkatesudu, B., Nagaraja, K.V., Shafiqul Islam, Md.: Gauss Legendre-Gauss Jacobi quadrature rules over a tetrahedral region. Appl. Math. Comp. 190, 186–194 (2007)

Hady Joumaa and Martin Ostoja-Starzewski Department of Mechanical Science and Engineering University of Illinois at Urbana-Champaign 1206 W. Green St. Urbana IL, 61801 USA e-mail: [email protected] Martin Ostoja-Starzewski e-mail: [email protected] (Received: December 30, 2010)