- Email: [email protected]
On the influence of material properties on the wave propagation in Mindlin-type microstructured solids
Wave Motion 50 (2013) 1127–1139
Contents lists available at ScienceDirect
Wave Motion journal homepage: www.elsevier.com/locate/wavemoti
On the influence of material properties on the wave propagation in Mindlin-type microstructured solids Andrus Salupere ∗ , Kert Tamm Centre for Nonlinear Studies, Institute of Cybernetics at Tallinn University of Technology, Akadeemia tee 21, 12618, Tallinn, Estonia
highlights • • • • •
Wave propagation in Mindlin-type microstructured solids is simulated numerically. The model equation includes nonlinearities in the macro- and microscale. The influence of material parameters on the wave propagation is demonstrated. Solitonic behaviour of emerged wave structures is demonstrated. The results could be applied in the nondestructive testing of material properties.
article
info
Article history: Received 30 December 2012 Accepted 30 May 2013 Available online 17 June 2013 Keywords: Microstructure Nonlinearity Solitary waves Dispersion
abstract A mathematical model describing 1D wave propagation in Mindlin-type microstructured solids with nonlinearities in the macro- and microscale is used for studying propagation of solitary waves in such media. The results could be used for the stress analysis as well as for the nondestructive testing of material properties. The model equations are solved numerically under the localized initial conditions and periodic boundary conditions by the pseudospectral method. It is demonstrated how the values of the model parameters influence the wave propagation, the evolution and the interaction of waves under the framework of considered models. For this reason the solutions of the model equations are compared under different parameter combinations against one fixed combination of material parameters which is called ‘the reference case’. © 2013 Elsevier B.V. All rights reserved.
1. Introduction The conventional theories of continua describe usually statics and dynamics of homogeneous materials, however, in the technological applications this can be inadequate for two reasons. Firstly, materials where homogeneous approach is overlooking important properties of the material are already widely used in the applications, like, for example, composites, functionally graded materials, granular materials and polycrystalline solids, which all have inherent microstructure at different scales. Secondly, high-frequency excitations can occur in the modern technology meaning that wavelengths of excitations are comparable with the internal scales in materials. This in turn means that in some applications the microstructural effects must be definitely taken into account, especially when modelling and analysing dynamical phenomena. As far as the materials with microstructure are already widely used in the applications, the benefits of increasing the knowledge of these materials can be significant. The key benefits from the research in this direction can be: (1) applications in the nondestructive testing of materials; (2) designing of new materials with predefined properties.
∗
Corresponding author. Tel.: +372 5249941. E-mail addresses: [email protected], [email protected] (A. Salupere), [email protected] (K. Tamm).
0165-2125/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.wavemoti.2013.05.004
1128
A. Salupere, K. Tamm / Wave Motion 50 (2013) 1127–1139
The influence of internal structure of solids on wave propagation in such media clearly manifests itself in the wave dispersion. The modelling of such effects can be traced to the Born–von Karman lattice model [1]. Following such an approach, the basic idea is the continualization of discrete systems which brings the higher-order terms in the governing equations into account. An excellent overview on corresponding mathematical models is given by Maugin [2]. The clear modelling of microstructured continua is presented by Mindlin [3] and Eringen [4] where the balance laws for macro- and microstructures are derived from a variational principles. Another possibility for describing microstructural effects is to use gradient elasticity models [5,6] which clearly includes both inertial and elastic properties of the microstructure, see also [7]. The Mindlin-type model is studied by Engelbrecht et al. [8] by using the Euler–Lagrange formalism for deriving the governing equations. Berezovski et al. [9] show that the same result can be achieved by using the formalism of internal variables. The comparative studies by Askes and Aifantis [10], Metrikine and Askes [5], Berezovski et al. [11], Papargyri-Beskou et al. [7] have demonstrated how the different mathematical models are related to each other. All these are characterized by the leading higher order terms UXXXX and UXXTT which are responsible for dispersive effects with the accuracy of their coefficients in different models. Actually the governing equations are of the Boussinesq type [12] grasping bi-directionality of waves, nonlinearity and dispersion. Clearly such models are able to support solitonic solutions. However, contrary to the classical one-wave models like the Korteweg–de Vries (KdV) equation, bi-directionality is an important factor which might introduce difficulties in the analysis. The Boussinesq solitons and their interaction processes are analysed by Christov et al. [12–14] by using numerical simulations (finite difference methods). Tamm [15] has analysed the structure, interaction and emergence of solitonic solutions in the Mindlin-type model. Engelbrecht et al. [16] have analysed the emergence of soliton trains and (inelastic) collision of solitons based on the same model. It must be stressed that asymmetric localized travelling wave solutions can exist in such models if the coefficients (combination of material parameters) of the governing equation(s) satisfy certain conditions [17]. That is why the influence of material parameters on solitonic solutions must be analysed in detail. Furthermore, we have demonstrated in our previous studies (see[15,18]) that the initial symmetric pulses can be morphed to these of asymmetric shape during propagation and interactions in the case of Mindlin-type model. Similar phenomena (reshaping and transformation of initial pulses and formation of different wave structures) are recently studied for different physical problems (see for example [19–21]). In the present paper the Mindlin-type mathematical model is used for describing the propagation of 1D longitudinal waves in microstructured solids. The influence of microstructure is described by higher-order terms in governing equations, nonlinearities are taken into account at both the macro- and microscale. The influence of material parameters on wave propagation in such materials is studied. For this end the model equations are integrated numerically making use of the pseudospectral method (PSM). In Section 2 of the present paper the governing equations are presented, the problem is stated and a numerical scheme is described briefly. In Section 3 the character of the combined parameter space is described and points of interest are highlighted therein. In Section 4 results are presented and discussed and in Section 5 conclusions are drawn. 2. Model equations, statement of the problem and the numerical method As mentioned in the Introduction, the 1D Mindlin-type model is used which is based on Mindlin’s and Eringen’s earlier works [3,4]. It should be noted that similar approach is also taken by Achenbach et al. [22] when describing wave propagation in laminated medium. The model equations are derived by Engelbrecht and Pastrone (see [8,23]) and here the derivation is shortly summarized. 2.1. Model equations The basic steps of deriving the model equations start from the Lagrangian L
L = K − W,
1
K =
2
1
ρ u2t + I ϕt2 ,
W = W (ux , ϕ, ϕx ),
2
(1)
where K is the kinetic energy and W is the free energy. The subscripts denote partial derivatives with the respect to subscripted variables, i.e., ux = ∂ u(x, t )/∂ x, while ϕ denotes microdeformation and u denotes macro-displacement. The free energy W is determined as W =
A 2
u2x +
B 2
C
ϕ2 +
2
ϕx2 + Dϕ ux +
N 6
u3x +
M 6
ϕx3 .
(2)
Here A, B, C , D are material parameters responsible for the linear part of the model and N , M are related to the nonlinearity in the macro- and microscale, respectively [8,17]. The Euler–Lagrange equations are used for deriving governing equations. Taking relevant derivatives from the Lagrangian and taking into account the free energy (2) one arrives at equations of motion. For further analysis dimensionless variables and parameters
√ X =
x Lo
,
T = √
At
ρ Lo
,
U =
u Uo
,
δ=
l2o L2o
,
ε=
Uo Lo
(3)
A. Salupere, K. Tamm / Wave Motion 50 (2013) 1127–1139
1129
are introduced [17]. Here Uo and Lo are the amplitude and the wavelength of the initial excitation, and lo is the characteristic scale of the microstructure. The equations of motion can be thus written in the dimensionless form as UTT = UXX + s1 ϕX + µUX UXX , s1 =
D
,
s2 =
Cρ
,
s=
ϕTT = s2 ϕXX − sϕ − α UX + νϕX ϕXX , Bρ L2o
,
Dρ Uo Lo
α=
,
ν=
Mρ
,
µ=
Nε
.
(4)
Aε IA IA IA IALo A Eqs. (4) are referred to as the full system of equations (the FSE for short) below. By making use of the slaving principle (see [8,24] for details) it is possible to derive a single equation in terms of dimensionless macrodisplacement U from the governing equations: UTT − bUXX −
b=1−
D2 AB
,
µ
UX2 X
2
= δ β UTT − γ UXX + ID2
β=
ρ l2o B2
,
CD2
γ =
AB2 l2o
√ λ δ 2
,
2 UXX
(5)
XX
,
µ=
Nε A
,
λ=
D3 M ε AB3 l3o
.
Eq. (5) can be considered as an approximation of full system of Eqs. (4) and the wave operator at the right hand side characterizes the influence of the microstructure on wave motion. Eq. (5) is referred to as the hierarchical equation (HE) below as it is hierarchical in Whitham’s sense [25]. 2.2. Statement of the problem The main goal of the present study is to investigate how material parameters influence the wave propagation for the considered models. More specifically, the goal of the present study is to demonstrate how the free energy (2) parameters influence the wave evolution (shape, propagation speed, interaction characteristics, etc.). Parameters B and C are purely related to the (linear) properties of the microstructure. Parameter D is related to the interaction between (linear) macroand microstructural properties and parameters M and N are related to the nonlinear properties of the micro- and as macrostructure, respectively. Parameter A is not varied as it is related the bulk medium on the macroscopic scale and is thus assumed to be fixed. 2.3. Initial and boundary conditions, numerical method Eqs. (4) and (5) are solved using sech2 -type localized initial conditions and periodic boundary conditions U (X , 0) = Uo sech2 Bo X ,
U (X , T ) = U (X + 2Kmπ , T ),
m = 1, 2, . . . ,
(6)
where K = 24, i.e., the total length of the spatial period is 48π . The initial phase speed is taken to be zero, which can be interpreted as starting from the peak of the interaction of two waves propagating in opposite directions. For the FSE two more initial conditions are needed for the microdeformation. We assume that at T = 0 the microdeformation and the corresponding velocity are zero, i.e., ϕ(X , 0) = 0 and ϕT (X , 0) = 0. The integration interval is from zero to Tf = 480. In all considered cases two solitary waves that propagate in opposite directions emerge from the initial pulse (6). In the present paper the discrete Fourier transform (DFT) based PSM is used [26,27]. The version of the DFT used is:
U (k, T ) = F [U] =
n −1
U (j∆X , T ) exp −
j =0
2π ijk
n
,
(7)
where n is the number of space-grid points, ∆X = 2π /n is the space step, k = 0, ±1, ±2, . . . , ±(n/2 − 1), −n/2; i is the imaginary unit, F denotes the DFT and F−1 denotes the inverse DFT. Basically, the idea of the PSM is to approximate space derivatives by making use of the DFT
∂ mU = F−1 (ik)m F(U ) , m ∂X
(8)
reducing therefore PDE to ODE and then to use standard ODE solvers for integration with respect to time. The regular PSM algorithm is derived for ut = Φ (u, ux , u2x , . . . , umx ) type equations. In our case, however, we have also a mixed partial derivative term δβ UTTXX in the HE (5) and thus the standard PSM has to be modified [27–29]. Therefore we rewrite the HE (5) so that all partial derivatives with respect to time are in the left-hand side of the HE UTT − δβ UTTXX = bUXX +
µ 2
UX2 X
− δ γ UXX −
√ λ δ 2
,
2 UXX
(9)
XX
and introduce a new variable Φ = U − δβ UXX . After that, making use of properties of the DFT, one can express the variable U and its spatial derivatives in terms of the new variable Φ : −1
U =F
F(Φ ) 1 + δβ k2
,
m ∂ mU −1 (ik) F(Φ ) =F . ∂Xm 1 + δβ k2
(10)
1130
A. Salupere, K. Tamm / Wave Motion 50 (2013) 1127–1139
Finally, Eq. (5) can be rewritten in terms of the variable Φ
ΦTT = bUXX +
µ 2
UX2 X
− δ γ UXX −
√ λ δ 2
.
2 UXX
(11)
XX
In Eq. (11) all partial derivatives of U with respect to X are calculated in terms of Φ by using expression (10) and therefore one can apply the PSM for numerical integration of Eq. (11). The FSE (4) is reduced to the system of first-order differential equations which are solved by the standard PSM without any further modifications. The calculations are carried out with the Python package SciPy [30], using the FFTW library [31] for the DFT and the F2PY [32] generated Python interface to the ODEPACK Fortran code [33] for the ODE solver. 3. Material parameters and dispersion In order to reduce the dimension of the material parameters space the material parameters in the free energy expression (2) which are responsible for the linear part of the model are collected into two combined parameters γA2 and γ12 which have clear physical interpretation. 3.1. Combined parameters Parameters γA2 and γ12 combine the linear parts of the material parameters (see [8,24]) describing macro- and microstructure and interaction between those
γA2 = 1 − b =
D2 AB
,
γ12 = s2 =
γ ρC = . β AI
(12)
Parameter γA2 is related to the speed of long waves in the system and γ12 to the speed of short waves (Which feel the influence of the microstructure) in the system. The parameters describing nonlinearity are used separately. Combined parameter γA2 has a physical meaning in the domain from 0 up to 1. In the limiting cases, if γA2 = 0 then the HE (5) is reduced into the classical wave equation and if γA2 = 1 then the normalized speed of the long waves would be zero. 3.2. Dispersion Dispersive properties of a material provide insight into the internal structure and degrees of freedom in the system. We have the HE (5) which is derived from the FSE (4) through slaving principle. The HE has only a single dispersion curve (branch) while the dispersion curve of the FSE maintains two branches — so called acoustic and optical (higher-order) branches, where the former provides insight about scale of the underlying microstructure and the latter is related to the internal degrees of freedom (see [34–36]). The dispersion type for the HE (and for the acoustic dispersion branch of the FSE) can be determined by the sign of quantity (see [8] for details):
Γ = 1 − γA2 − γ12 .
(13)
If Γ is positive, we have the normal dispersion case, if negative, we have the anomalous dispersion case and if it is equal to zero, we have the dispersionless case. For the optical dispersion branch of the FSE the dispersion type is always anomalous. It has been shown by Metrikine in [37] that the higher order dispersion branch is essential in such models for achieving causality. In essence, without the higher order dispersion branch an infinite propagation speed for some disturbances in such models may appear. It should be noted that in the dispersionless case the acoustic dispersion branches of the FSE and HE coincide. For the sake of clarity — the classical way of dealing with dispersion types is to compare the group and phase speeds in the system. If these are equal to each other then we have dispersionless case. If the group speed is larger than the phase speed (cgr > cph ) then we have the anomalous dispersion case and if the group speed is smaller than the phase speed (cgr < cph ) then the normal dispersion case. In terms of wave propagation the essence of a dispersion is that the phase speed of the wave depends on the wave number of the harmonics. This is due to the length scale that the underlying microstructure brings into the models. In Fig. 1 a combined parameter domain and points of interest therein are presented. The dashed line represents parameter combinations corresponding to the dispersionless case, above that line the dispersion type is anomalous while below the dashed line the dispersion type is normal. The crosses, a ring and the triangles in Fig. 1 represent the investigated points in the combined parameter space for which the particular parameter combinations are presented later in the paper. The domain between thick solid lines corresponds to parameter combinations where concurrency between solutions of the HE and FSE is found to be good according to the dispersion analysis and numerical simulations [15,24]. In the area between the thick solid lines in Fig. 1 the difference between acoustical dispersion curves for the FSE and HE is 5% or less at dimensionless wavenumber of 1.5 (see Fig. 2 for examples of dispersion curves). In Fig. 2 some dispersion curves are presented for various points in the plane of combined parameters. Two upper subfigures give examples for two points in the horizontal line of crosses and the two lower subfigures give examples for
A. Salupere, K. Tamm / Wave Motion 50 (2013) 1127–1139
1131
Fig. 1. The combined parameter plane and points of interest corresponding to different values of parameters B, C and D.
end-points of the vertical line of triangles in Fig. 1. Top left and bottom right subfigures correspond to the normal dispersion case and top right and bottom left subfigures to these of the anomalous dispersion (with regard to the acoustic branch). 4. Results and discussion Material parameters are taken as A = 20,
B=
N = 50 ± 25,
D2 AγA2
,
C = 11 ± 5,
D = 10 ± 5, (14)
M = 5000 ± 2500,
ρ = 1,
I =
4 5
,
where parameter B is picked up to visit the points of interests defined by the parameter D in the combined parameter space i.e., parameter B has values of 60,
125 3
, 1500 , 375 , 500 , 15, 1500 , 125 , 1500 , 375 , 20 49 16 27 121 12 169 49 3
. Parameters C and D are investigated
in the given range with the step size 1 while parameter N is investigated with the step size 5 and parameter M with the step size 500. The points of interest corresponding to these parameter combinations in the combined parameter space can be seen in Fig. 1. In Fig. 1 the ‘reference case’ corresponds to the parameters combination against what other cases are compared to spot the effects which can caused by increasing or decreasing of material parameters. In the reference case B = 15, C = 11, D = 10, N = 50, M = 5000. The geometrical parameters are taken lo = 1,
Lo = 50,
Uo = 1,
Bo =
1 4
,
Ao = 1,
(15)
where lo is the characteristic scale of the microstructure, Lo — the characteristic wavelength of the initial excitation and Uo — the amplitude of the initial excitation while Bo — the parameter controlling the width of the initial pulse and Ao — the normalized amplitude of the initial pulse. The values above are taken for a certain hypotetic functionally graded material [15]. A reference case of formation, propagation and interaction of two solitary waves from a localized initial pulse in case of the FSE is shown in Fig. 3. The initial bell-shaped pulse splits into two solitary pulses propagating in the opposite directions which evolve differently in time — the pulse propagating to the left evolves towards a peakon-like shape [18,38,39] while the pulse propagating to the right maintains its bell-like shape which can be spotted also in the example solution in Fig. 3 at closer scrutiny. Such a asymmetry is caused by the internal structure of the model equations. In Fig. 4 one can see the phase plot corresponding to the waveprofiles at the end of integration (dashed and dotted lines) and the initial solitary wave
1132
A. Salupere, K. Tamm / Wave Motion 50 (2013) 1127–1139
Fig. 2. Dispersion curves for the FSE and HE at γA2 = 0.168, γ12 = 0.688 (top left), γA2 = 0.48, γ12 = 0.688 (top right), γA2 = 0.3333, γ12 = 1 (bottom left), and γA2 = 0.333, γ12 = 0.375 (bottom right). Variable η is the dimensionless frequency and variable ξ the dimensionless wavenumber.
Fig. 3. The typical solution of Eq. (4).
(solid line) — on horizontal axis there is a spatial derivative of the displacement (UX ) and on vertical axis the displacement amplitude of the wave (U ). A term ‘normalized initial wave’ is used for the waveprofile at a short time after the start where the initial pulse has divided into the two solitary pulses propagating in opposite directions but still having almost exactly sech2 -type profile, just with an amplitude of 0.5 instead of the initial amplitude of one. The main characteristics under investigation are the wave speeds which are tracked by finding the trajectory of the pulse peak and calculating the speed from these and the wave shape which is investigated mainly by making use of the phase-plots where even relatively small changes in the pulse shape are clearly visible.
A. Salupere, K. Tamm / Wave Motion 50 (2013) 1127–1139
1133
Fig. 4. The phase plot corresponding to the reference case at T = 480. Solid line — normalized initial wave, (blue) dotted line — peakon-like pulse, (red) dashed line — bell-like pulse.
Table 1 Calculated average speeds cB and cD and characteristic speed co against parameters B and D. D
cD
B
cB
γA2
co
cD - co
5 6 7 8 9 10 11 12 13 14 15
0.960 0.939 0.916 0.888 0.857 0.819 0.775 0.725 0.662 0.593 0.502
60 125/3 1500/49 375/16 500/27 15 1500/121 125/12 1500/169 375/49 20/3
0.960 0.939 0.917 0.888 0.856 0.818 0.775 0.723 0.662 0.590 0.501
0.083 0.120 0.163 0.213 0.270 0.333 0.403 0.480 0.563 0.653 0.750
0.917 0.880 0.837 0.787 0.730 0.667 0.597 0.520 0.437 0.347 0.250
0.043 0.059 0.080 0.102 0.127 0.152 0.179 0.205 0.226 0.246 0.252
4.1. Parameters B and D The average speeds in Table 1 have been calculated for the pulse propagating to the right in the interval from separating from the initial pulse until the first interaction for the FSE. The speeds for the HE are practically the same and are within ±0.02 which is the accuracy of the speeds in the table. In Table 1, cB denotes the speed calculated for the cases where parameter B was used to control the value of the combined parameter γA2 and cD denotes the speed calculated for the cases where parameter D was used to control the value of the combined parameter γA2 . The speed co is the characteristic speed for √ the bulk medium and can be found as co = A/ρ . The speed differences between the solutions of the FSE (4) and HE (5) as well as between speeds corresponding to the parameters B and D resulting in the same combined parameter γA2 values are negligible. What can be noted, however, is that with the increasing parameter γA2 the gap between the speeds cB , cD and the characteristic speed co increases. While the parameters D and B have similar influence on the wave propagation speed their effect on the pulse shape is slightly different. At lower values of the combined parameter γA2 the differences are hard to spot but at highest value of D or lowest value of B (D = 15 or B = 20/3) one can see from the phase plots in Fig. 5 that in the case of high D the oscillatory structure forming in the peakon like profile has higher frequency oscillations than in the case of low values of parameter B. Regardless of the observed differences qualitatively the oscillatory structure is the same as for both high D and low B values altogether 5 oscillations can be distinguished in Fig. 5 on the peakon-like pulse. No notable differences can be detected in the bell-like pulse for the high values of parameter D or low values of parameter B.
1134
A. Salupere, K. Tamm / Wave Motion 50 (2013) 1127–1139
84 83 82 81 80 79 78 77 76 75 74 47π
T
156 155 154 153 152 151 150 149 148 147 146 47π
B = 15; D = 15 B = 20/3; D = 10 47.5π
48π
48.5π
49π
T
T
T
, D = 10, T = 480 (top right); B = 15, D = 5, T = 460 (bottom left); B = 15, Fig. 5. Phase plots for B = 60, D = 10, T = 460 (top left); B = 20 3 D = 15, T = 480 (bottom right). Solid line — normalized initial wave, (blue) dotted line — peakon-like pulse, (red) dashed line — bell-like pulse, and parameters C , M and N have reference values.
B = 15; D =5 B = 60; D = 10 47.5π
48π 47π ≤ X < 49π
48.5π
49π
307 306 305 304 303 302 301 300 299 298 297 23π 163 162 161 160 159 158 157 156 155 154 153 23π
B = 15; D = 15 B = 20/3; D = 10 23.5π
24π
25.5π
25π
B = 15; D =5 B = 60; D = 10 23.5π
24π 23π ≤ X < 25π
25.5π
25π
Fig. 6. Trajectories of the first (left column) and the second (right column) interaction of solitary waves for different parameter B and D values. Top row — γA2 = 0.75, bottom row — γA2 = 0.083.
In order to clarify the influence of parameters B and D on the behaviour of the solution during interactions we trace wave-trajectories for γA2 = 0.083 and γA2 = 0.75 (see Fig. 6). According to Eq. (12) these values can be reached by different combinations of parameters B and D. For both values of γA2 two combinations of B and D are considered. It is clear that for smaller values of D compared trajectories practically coincide (lower subfigures), but for higher values of D small differences
A. Salupere, K. Tamm / Wave Motion 50 (2013) 1127–1139
1135
Fig. 7. Phase plots for C = 6 (left) and C = 16 (right) at T = 480. Solid line — normalized initial wave, (blue) dotted line — peakon-like pulse, (red) dashed line — bell-like pulse, and parameters B, D, M and N have reference values.
can be detected (higher subfigures). It should be mentioned that in here the solitary pulses do not experience phase shifts during the interactions due to the normalization procedure used which removes the usual speed–amplitude dependence which exists in the Boussinesq- and KdV-type equations (see, for example [12,40]). 4.2. Parameter C The effect of parameter C on the calculated average speeds is negligible and it has no influence on the observed interaction trajectories of the solitary pulses. In Fig. 7 one can see the two phase-plots representing waveprofile shape at the end of integration at T = 480 corresponding to the different values of the parameter C . The most notable effect on the phase plots is that the oscillatory structure which develops by the end of the integration interval switched from the ‘front’ (in the direction of propagation, anomalous dispersion) to the ‘back’ (opposite to the direction of propagation, normal dispersion) of the sharper pulse. The closer the point is in the γA2 − γ12 plane to the dispersionless case the sharper the peakonlike pulse. Strictly speaking a peakon should have a discontinuity at its peak but our numerical scheme does not handle specifically discontinuities so the integration interval has been picked short enough to avoid reaching a true discontinuity at the sharper profile peak in the dispersionless case. Parameter C in the present study is used to control the values of combined parameter γ12 which is interpreted as a speed of short waves which feel the influence of the microstructure. Under the used parameter combinations, this translates into controlling the dispersion type for the acoustic dispersion branch of the governing equations. More specifically, looking at Eqs. (12) and (13) it is clear that parameter C has a significant impact on the speed of harmonics with a shorter wavelength and dispersion type which would explain the observed differences in the waveprofile shapes at different parameter C values. 4.3. Nonlinear parameters N and M The effect of parameters N and M on the calculated average speeds is negligible. In Fig. 8 the highest presented parameter N value is N = 55 as at higher parameter N values the wave starts breaking because of the peakon-like profile peak gets too sharp as to interfere with the numerical scheme used. At the highest presented value of parameter N one can already see emerging high frequency oscillations from the discontinuity at the waveprofile peak. The lower the parameter N the smaller the difference between waveprofiles propagating in the opposite directions and vice versa — the higher the parameter N the greater the difference (Fig. 8) — the sharper pulse gets thinner and sharper while the smoother pulse gets even thicker with increasing N. Overall the parameter M influences on the waveprofile evolution under the used parameter combinations can be considered weak (see Fig. 8). However, one can note that even at the lowest considered parameter M value the beginnings of the oscillatory structure can be spotted in the phase-plot corresponding to the peakon-like pulse. Compared to the parameter N the influence of parameter M on the waveprofile shape is negligible under the investigated parameter combinations. It is known from previous studies (see, for example [41] where its shown for the HE (5)) that under some parameter combinations the parameter M can affect the asymmetry of a waveprofile so that it ‘tilts’ in the direction of propagation or opposite of the direction of propagation (according to on dispersion type for such parameter combination). However, it must be added that in here the conditions described in [41] are not fulfilled. A change in the sharper pulse shape at highest investigated value of parameter M can be noted in Fig. 8 where the emerging oscillations appear to be slightly stronger than in the case of weakest considered parameter M (Fig. 8) or the reference case (Fig. 4).
1136
A. Salupere, K. Tamm / Wave Motion 50 (2013) 1127–1139
Fig. 8. Phase plots for M = 5000, N = 25 (top left), M = 5000, N = 55 (top right) and M = 2500, N = 50 (bottom left) and M = 7500, N = 50 (bottom right) at T = 480. Solid line — normalized initial wave, (blue) dotted line — peakon-like pulse, (red) dashed line — bell-like pulse, and parameters B, C and D have reference values.
Another notable difference between the influence of nonlinear parameters on the waveprofile evolution emerges when we investigate the apparent trajectories of the waveprofiles by tracking their maxima. In the case of parameter N it is possible to observe a clear difference between trajectories drawn by the waveprofile peak during interactions at different parameter N values (Fig. 9). Outside the interactions the values of parameter N have no visible effect on the trajectories. At closer inspection (see Fig. 10) an interesting phenomenon appears: the pulse itself is practically stationary during the interactions, but the peak of the waveprofile moves. At the beginning of the interaction a combined waveprofile is asymmetric and tilted to the left. During the interaction it reaches symmetric shape at the peak of interaction and is tilted to the right after the peak of interaction until the waveprofiles separate again after the interaction. In addition, in Fig. 10 it is possible to spot that during the fourth interaction the waveprofile has sharper peak than during the second interaction. Of course, it should be useful to perform rigorous analysis of interactions for the present model, like it is done for the KdV equation [42]. 4.4. Concurrency of the solutions In order to calculate the difference between the solutions of the HE and FSE the quantity ∆S is introduced as
∆S =
n 1
n i =1
∆i ,
∆i = U HE (Xi , Tf ) − U FSE (Xi , Tf ) ,
(16)
where n is the number of grid points. In the sense of the quantity ∆S for the parameters C , N and M the concurrency of the solutions of the HE (5) and FSE (4) does not depend on the parameter values, however, for the parameters B and D the parameter value also affects the concurrency of the solutions. The concurrency of the solutions gets worse with decreasing
190 189 188 187 186 185 184 183 182 181 180 179 22π
T
Ν=25 Ν=50 Ν=75
23π
190 189 188 187 186 185 184 183 182 181 180 179 22π
24π
25π
26π
Μ=2500 Μ=5000 Μ=7500
T
T
T
A. Salupere, K. Tamm / Wave Motion 50 (2013) 1127–1139
23π
24π
25π
26π
374 373 372 371 370 369 368 367 366 365 364 22π
1137
Ν=25 Ν=50 Ν=75
23π
24π
25π
374 373 372 371 370 369 368 367 366 365 364
26π
Μ=2500 Μ=5000 Μ=7500
22π
23π
24π
25π
26π
1.0 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 18π
30π
1.0 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 18π
30π
1.0 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 18π
U, T = 183.20 U, T = 184.00 U, T = 184.80 U, T = 185.60 U, T = 186.40
U
1.0 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 18π
20π
22π
24π
26π
28π
U, T = 183.20 U, T = 184.00 U, T = 184.80 U, T = 185.60 U, T = 186.40
U
U
U
Fig. 9. Trajectories of the second (left) and the fourth (right) interactions for the nonlinear parameter values N = 25, 50, 75 and M = 5000 (top) and M = 2500, 5000, 7500 and N = 50 (bottom). Parameters B, C and D have reference values.
20π
22π
24π X
26π
28π
U, T = 367.80 U, T = 368.60 U, T = 369.40 U, T = 370.20 U, T = 371.00
20π
22π
24π
26π
28π
30π
24π X
26π
28π
30π
U, T = 367.80 U, T = 368.60 U, T = 369.40 U, T = 370.20 U, T = 371.00
20π
22π
Fig. 10. The second (left) and the fourth (right) interactions for the nonlinear parameter values N = 75, M = 5000 (top) and N = 50, M = 7500 (bottom). Parameters B, C and D have reference values.
of the parameter B or increasing the parameter D. Under the investigated parameter combinations the quantity log10 (∆S ) is between −5 and −7 which we consider to be a good concurrency between the solutions of the FSE (4) and HE (5) because this accuracy is sufficient to investigate in detail do the solutions of the FSE and HE behave in a similar way or not. If quantity log10 (∆S ) is close to −3 then the concurrency of the solutions cannot be considered anymore satisfactory as there can be already substantial differences between the solutions of the FSE (4) and HE (5) (see [15] for details) after sufficiently short evolution in time.
1138
A. Salupere, K. Tamm / Wave Motion 50 (2013) 1127–1139
5. Conclusions In the present paper we studied wave propagation in solids having a Mindlin-type microstructure. The considered free energy function (2) includes material parameters of different kind. Parameter A characterizes linear properties of the material at the macroscale and parameters B and C characterize these at the microscale. Interaction between macro- and microscale effects is described by the parameter D. In their turn, parameters M and N characterize the influence of the nonlinear properties of the micro- and macroscale. We were particularly interested how parameters B, C , D, M and N influence the propagation and evolution of solitary waves. To summarize our findings following conclusions can be drawn: (i) The speed of emerged solitary waves is directly related to values of parameters B and D — the higher the value of B or the lower value of D the higher the speed of waves. The influence of parameters C , M and N on the wave speed can be neglected. (ii) Parameters C and M influence the character of asymmetry of solitary pulses (including the location of oscillatory structures with respect to the direction of propagation of main pulses). In the considered domain of parameters the influence of M is weak, but according to earlier studies [17,41] the character of asymmetry of the travelling wave solution of the HE depends directly on the value of M. (iii) Parameters D and N influence the shape of single pulses in the following way: the greater the parameter D or N the sharper the peakon-like pulse and the wider the bell-like pulse. Parameter B has the opposite effect on the pulse-shapes. (iv) Parameter N influences the rate of the change of asymmetry during the interaction process — the higher the value of N the higher the rate of the change of asymmetry (i.e. the faster the peak of the waveprofile moves during the interaction). It influences also the trajectories within the interaction process. Last but not least, if the propagation of solitary pulses in nonlinear dispersive media is studied then the question about solitonic character of the wave structure appears. Notwithstanding that the interactions between solitary pulses are not completely elastic (the peakon-like pulse gets sharper and the bell-like pulse gets wider and oscillatory structures can be generated due to interactions), the main pulses propagate at a constant speed and almost restore their shape and speed after interactions. Therefore we conclude that the emerged wave structures modelled by a Boussinesq-type equation (5) or its full form (4) have solitonic behaviour (at least throughout several first interactions). Acknowledgements The research was supported by the Estonian Science Foundation Grant No. 8658, the Estonian Ministry of Education and Research (block grant SF0140077s08) and by the European Union through the Regional Development Fund (projects CENS and SCDM). The authors thank Prof. J. Engelbrecht for helpful discussions and suggestions related to this paper. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27]
M. Born, T. von Kármán, Uber Schwingungen in Raumgittern, Phys. Z. 13 (1) (1912) 297–309. G.A. Maugin, Nonlinear Waves in Elastic Crystals, Oxford Univercity Press, Oxford, 1999. R. Mindlin, Micro-structure in linear elasticity, Arch. Ration. Mech. Anal. 16 (1) (1964) 51–78. A. Eringen, Linear theory of micropolar elasticity, J. Math. Mech. 15 (1966) 909–923. A. Metrikine, H. Askes, One-dimensional dynamically consistent gradient elasticity models derived from a discrete microstructure, part 1: generic formulation, Eur. J. Mech. A Solids 21 (2002) 555–572. V. Erofeev, Wave Processes in Solids with Microstructure, World Scientific, Singapore, 2003. S. Papargyri-Beskou, D. Polyzos, D. Beskos, Wave dispersion in gradient elastic solids and structures: a unified treatment, Int. J. Solids Struct. 46 (21) (2009) 3751–3759. J. Engelbrecht, A. Berezovski, F. Pastrone, M. Braun, Waves in microstructured materials and dispersion, Phil. Mag. 85 (33–35) (2005) 4127–4141. A. Berezovski, J. Engelbrecht, G.A. Maugin, Generalized thermomechanics with dual internal variables, Arch. Appl. Mech. 81 (2) (2010) 229–240. H. Askes, E.C. Aifantis, Gradient elasticity theories in statics and dynamics — a unification of approaches, Int. J. Fract. 139 (2) (2006) 297–304. A. Berezovski, J. Engelbrecht, M. Berezovski, Waves in microstructured solids: a unified viewpoint of modeling, Acta Mech. 220 (1–4) (2011) 349–363. C.I. Christov, G.A. Maugin, A.V. Porubov, On Boussinesq’s paradigm in nonlinear wave propagation, C. R. Mecanique 335 (9–10) (2007) 521–535. C.I. Christov, M. Velarde, Inelastic interaction of Boussinesq solitons, Internat. J. Bifur. Chaos Appl. Sci. Engrg. 4 (5) (1994) 1095–1112. C.I. Christov, G.A. Maugin, M. Velarde, Well-posed Boussinesq paradigm with purely spatial higher-order derivatives, Phys. Rev. E 54 (4) (1996) 3621–3638. K. Tamm, Wave propagation and interaction in mindlin-type microstructured solids: numerical simulation, Ph.D. Thesis, Tallinn University of Technology, 2011. J. Engelbrecht, A. Salupere, K. Tamm, Waves in microstructured solids and the Boussinesq paradigm, Wave Motion 48 (2011) 717–726. J. Janno, J. Engelbrecht, Solitary waves in nonlinear microstructured materials, J. Phys. A: Math. Gen. 38 (2005) 5159–5172. K. Tamm, A. Salupere, On the propagation of 1D solitary waves in Mindlin-type microstructured solids, Math. Comput. Simul. 82 (7) (2012) 1308–1320. A. Sergeeva, E. Pelinovsky, T. Talipova, Nonlinear random wave field in shallow water: variable Korteweg–de Vries framework, Nat. Hazards Earth Syst. Sci. 11 (2) (2011) 323–330. A.V. Porubov, G.A. Maugin, B.R. Andrievsky, Solitary wave interactions and reshaping in coupled systems, Wave Motion 48 (8) (2011) 773–781. A.V. Porubov, I.V. Andrianov, V.V. Danishevskyy, Nonlinear strain wave localization in periodic composites, Int. J. Solids Struct. 49 (2012) 3381–3387. C. Sun, J.D. Achenbach, G. Herrmann, Continuum theory for a laminated medium, J. Appl. Mech. 35 (3) (1968) 467. J. Engelbrecht, F. Pastrone, Waves in microstructured solids with nonlinearities in microscale, Proc. Est. Acad. Sci., Phys. Math. 52 (1) (2003) 12–21. T. Peets, M. Randrüüt, J. Engelbrecht, On modelling dispersion in microstructured solids, Wave Motion 45 (4) (2008) 471–480. G. Whitham, Linear and Nonlinear Waves, Wiley-Interscience, New York, 1974. B. Fornberg, A Practical Guide to Pseudospectral Methods, Cambridge University Press, Cambridge, 1998. A. Salupere, The pseudospectral method and discrete spectral analysis, in: E. Quak, T. Soomere (Eds.), Applied Wave Mathematics, Springer, Berlin, 2009, pp. 301–334.
A. Salupere, K. Tamm / Wave Motion 50 (2013) 1127–1139
1139
[28] L. Ilison, A. Salupere, Propagation of sech2 -type solitary waves in hierarchical KdV-type systems, Math. Comput. Simul. 79 (11) (2009) 3314–3327. [29] L. Ilison, A. Salupere, P. Peterson, On the propagation of localized perturbations in media with microstructure, Proc. Est. Acad. Sci., Phys. Math. 56 (2) (2007) 84–92. [30] E. Jones, T. Oliphant, P. Peterson, SciPy: open source scientific tools for Python, 2007. http://www.scipy.org. [31] M. Frigo, S. Johnson, The design and implementation of FFTW 3, Proc. IEEE 93 (2) (2005) 216–231. [32] P. Peterson, F2PY: Fortran to Python interface generator, 2005. http://cens.ioc.ee/projects/f2py2e/. [33] A. Hindmarsh, ODEPACK, A Systematized Collection of ODE Solvers, Vol. 1, North-Holland, Amsterdam, 1983. [34] L. Brillouin, Wave Propagation in Periodic Structures, Dover Publications Inc., New York, 1953. [35] A. Eringen, Microcontinuum Field Theories I: Foundations and Solids, Springer, New York, 1999. [36] K. Tamm, T. Peets, On the influence of internal degrees of freedom on dispersion in microstructured solids, Mech. Res. Commun. 47 (2013) 106–111. [37] A. Metrikine, On causality of the gradient elasticity models, J. Sound Vib. 297 (3–5) (2006) 727–742. [38] R. Camassa, D. Holm, An integrable shallow water equation with peaked solitons, Phys. Rev. Lett. 71 (11) (1993) 1661–1664. [39] H. Kalisch, J. Lenells, Numerical study of traveling-wave solutions for the Camassa–Holm equation, Chaos Solitons Fractals 25 (2) (2005) 287–298. [40] A. Salupere, On hidden solitons in KdV related systems, Math. Comput. Simulat. (2013) submitted for publication. [41] J. Janno, J. Engelbrecht, An inverse solitary wave problem related to microstructured materials, Inverse Problems 21 (2005) 2019–2034. [42] A. Salupere, J. Engelbrecht, P. Peterson, Long-time behaviour of soliton ensembles, part I—emergence of ensembles, Chaos Solitons Fractals 14 (2002) 1413–1424.
Contents lists available at ScienceDirect
Wave Motion journal homepage: www.elsevier.com/locate/wavemoti
On the influence of material properties on the wave propagation in Mindlin-type microstructured solids Andrus Salupere ∗ , Kert Tamm Centre for Nonlinear Studies, Institute of Cybernetics at Tallinn University of Technology, Akadeemia tee 21, 12618, Tallinn, Estonia
highlights • • • • •
Wave propagation in Mindlin-type microstructured solids is simulated numerically. The model equation includes nonlinearities in the macro- and microscale. The influence of material parameters on the wave propagation is demonstrated. Solitonic behaviour of emerged wave structures is demonstrated. The results could be applied in the nondestructive testing of material properties.
article
info
Article history: Received 30 December 2012 Accepted 30 May 2013 Available online 17 June 2013 Keywords: Microstructure Nonlinearity Solitary waves Dispersion
abstract A mathematical model describing 1D wave propagation in Mindlin-type microstructured solids with nonlinearities in the macro- and microscale is used for studying propagation of solitary waves in such media. The results could be used for the stress analysis as well as for the nondestructive testing of material properties. The model equations are solved numerically under the localized initial conditions and periodic boundary conditions by the pseudospectral method. It is demonstrated how the values of the model parameters influence the wave propagation, the evolution and the interaction of waves under the framework of considered models. For this reason the solutions of the model equations are compared under different parameter combinations against one fixed combination of material parameters which is called ‘the reference case’. © 2013 Elsevier B.V. All rights reserved.
1. Introduction The conventional theories of continua describe usually statics and dynamics of homogeneous materials, however, in the technological applications this can be inadequate for two reasons. Firstly, materials where homogeneous approach is overlooking important properties of the material are already widely used in the applications, like, for example, composites, functionally graded materials, granular materials and polycrystalline solids, which all have inherent microstructure at different scales. Secondly, high-frequency excitations can occur in the modern technology meaning that wavelengths of excitations are comparable with the internal scales in materials. This in turn means that in some applications the microstructural effects must be definitely taken into account, especially when modelling and analysing dynamical phenomena. As far as the materials with microstructure are already widely used in the applications, the benefits of increasing the knowledge of these materials can be significant. The key benefits from the research in this direction can be: (1) applications in the nondestructive testing of materials; (2) designing of new materials with predefined properties.
∗
Corresponding author. Tel.: +372 5249941. E-mail addresses: [email protected], [email protected] (A. Salupere), [email protected] (K. Tamm).
0165-2125/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.wavemoti.2013.05.004
1128
A. Salupere, K. Tamm / Wave Motion 50 (2013) 1127–1139
The influence of internal structure of solids on wave propagation in such media clearly manifests itself in the wave dispersion. The modelling of such effects can be traced to the Born–von Karman lattice model [1]. Following such an approach, the basic idea is the continualization of discrete systems which brings the higher-order terms in the governing equations into account. An excellent overview on corresponding mathematical models is given by Maugin [2]. The clear modelling of microstructured continua is presented by Mindlin [3] and Eringen [4] where the balance laws for macro- and microstructures are derived from a variational principles. Another possibility for describing microstructural effects is to use gradient elasticity models [5,6] which clearly includes both inertial and elastic properties of the microstructure, see also [7]. The Mindlin-type model is studied by Engelbrecht et al. [8] by using the Euler–Lagrange formalism for deriving the governing equations. Berezovski et al. [9] show that the same result can be achieved by using the formalism of internal variables. The comparative studies by Askes and Aifantis [10], Metrikine and Askes [5], Berezovski et al. [11], Papargyri-Beskou et al. [7] have demonstrated how the different mathematical models are related to each other. All these are characterized by the leading higher order terms UXXXX and UXXTT which are responsible for dispersive effects with the accuracy of their coefficients in different models. Actually the governing equations are of the Boussinesq type [12] grasping bi-directionality of waves, nonlinearity and dispersion. Clearly such models are able to support solitonic solutions. However, contrary to the classical one-wave models like the Korteweg–de Vries (KdV) equation, bi-directionality is an important factor which might introduce difficulties in the analysis. The Boussinesq solitons and their interaction processes are analysed by Christov et al. [12–14] by using numerical simulations (finite difference methods). Tamm [15] has analysed the structure, interaction and emergence of solitonic solutions in the Mindlin-type model. Engelbrecht et al. [16] have analysed the emergence of soliton trains and (inelastic) collision of solitons based on the same model. It must be stressed that asymmetric localized travelling wave solutions can exist in such models if the coefficients (combination of material parameters) of the governing equation(s) satisfy certain conditions [17]. That is why the influence of material parameters on solitonic solutions must be analysed in detail. Furthermore, we have demonstrated in our previous studies (see[15,18]) that the initial symmetric pulses can be morphed to these of asymmetric shape during propagation and interactions in the case of Mindlin-type model. Similar phenomena (reshaping and transformation of initial pulses and formation of different wave structures) are recently studied for different physical problems (see for example [19–21]). In the present paper the Mindlin-type mathematical model is used for describing the propagation of 1D longitudinal waves in microstructured solids. The influence of microstructure is described by higher-order terms in governing equations, nonlinearities are taken into account at both the macro- and microscale. The influence of material parameters on wave propagation in such materials is studied. For this end the model equations are integrated numerically making use of the pseudospectral method (PSM). In Section 2 of the present paper the governing equations are presented, the problem is stated and a numerical scheme is described briefly. In Section 3 the character of the combined parameter space is described and points of interest are highlighted therein. In Section 4 results are presented and discussed and in Section 5 conclusions are drawn. 2. Model equations, statement of the problem and the numerical method As mentioned in the Introduction, the 1D Mindlin-type model is used which is based on Mindlin’s and Eringen’s earlier works [3,4]. It should be noted that similar approach is also taken by Achenbach et al. [22] when describing wave propagation in laminated medium. The model equations are derived by Engelbrecht and Pastrone (see [8,23]) and here the derivation is shortly summarized. 2.1. Model equations The basic steps of deriving the model equations start from the Lagrangian L
L = K − W,
1
K =
2
1
ρ u2t + I ϕt2 ,
W = W (ux , ϕ, ϕx ),
2
(1)
where K is the kinetic energy and W is the free energy. The subscripts denote partial derivatives with the respect to subscripted variables, i.e., ux = ∂ u(x, t )/∂ x, while ϕ denotes microdeformation and u denotes macro-displacement. The free energy W is determined as W =
A 2
u2x +
B 2
C
ϕ2 +
2
ϕx2 + Dϕ ux +
N 6
u3x +
M 6
ϕx3 .
(2)
Here A, B, C , D are material parameters responsible for the linear part of the model and N , M are related to the nonlinearity in the macro- and microscale, respectively [8,17]. The Euler–Lagrange equations are used for deriving governing equations. Taking relevant derivatives from the Lagrangian and taking into account the free energy (2) one arrives at equations of motion. For further analysis dimensionless variables and parameters
√ X =
x Lo
,
T = √
At
ρ Lo
,
U =
u Uo
,
δ=
l2o L2o
,
ε=
Uo Lo
(3)
A. Salupere, K. Tamm / Wave Motion 50 (2013) 1127–1139
1129
are introduced [17]. Here Uo and Lo are the amplitude and the wavelength of the initial excitation, and lo is the characteristic scale of the microstructure. The equations of motion can be thus written in the dimensionless form as UTT = UXX + s1 ϕX + µUX UXX , s1 =
D
,
s2 =
Cρ
,
s=
ϕTT = s2 ϕXX − sϕ − α UX + νϕX ϕXX , Bρ L2o
,
Dρ Uo Lo
α=
,
ν=
Mρ
,
µ=
Nε
.
(4)
Aε IA IA IA IALo A Eqs. (4) are referred to as the full system of equations (the FSE for short) below. By making use of the slaving principle (see [8,24] for details) it is possible to derive a single equation in terms of dimensionless macrodisplacement U from the governing equations: UTT − bUXX −
b=1−
D2 AB
,
µ
UX2 X
2
= δ β UTT − γ UXX + ID2
β=
ρ l2o B2
,
CD2
γ =
AB2 l2o
√ λ δ 2
,
2 UXX
(5)
XX
,
µ=
Nε A
,
λ=
D3 M ε AB3 l3o
.
Eq. (5) can be considered as an approximation of full system of Eqs. (4) and the wave operator at the right hand side characterizes the influence of the microstructure on wave motion. Eq. (5) is referred to as the hierarchical equation (HE) below as it is hierarchical in Whitham’s sense [25]. 2.2. Statement of the problem The main goal of the present study is to investigate how material parameters influence the wave propagation for the considered models. More specifically, the goal of the present study is to demonstrate how the free energy (2) parameters influence the wave evolution (shape, propagation speed, interaction characteristics, etc.). Parameters B and C are purely related to the (linear) properties of the microstructure. Parameter D is related to the interaction between (linear) macroand microstructural properties and parameters M and N are related to the nonlinear properties of the micro- and as macrostructure, respectively. Parameter A is not varied as it is related the bulk medium on the macroscopic scale and is thus assumed to be fixed. 2.3. Initial and boundary conditions, numerical method Eqs. (4) and (5) are solved using sech2 -type localized initial conditions and periodic boundary conditions U (X , 0) = Uo sech2 Bo X ,
U (X , T ) = U (X + 2Kmπ , T ),
m = 1, 2, . . . ,
(6)
where K = 24, i.e., the total length of the spatial period is 48π . The initial phase speed is taken to be zero, which can be interpreted as starting from the peak of the interaction of two waves propagating in opposite directions. For the FSE two more initial conditions are needed for the microdeformation. We assume that at T = 0 the microdeformation and the corresponding velocity are zero, i.e., ϕ(X , 0) = 0 and ϕT (X , 0) = 0. The integration interval is from zero to Tf = 480. In all considered cases two solitary waves that propagate in opposite directions emerge from the initial pulse (6). In the present paper the discrete Fourier transform (DFT) based PSM is used [26,27]. The version of the DFT used is:
U (k, T ) = F [U] =
n −1
U (j∆X , T ) exp −
j =0
2π ijk
n
,
(7)
where n is the number of space-grid points, ∆X = 2π /n is the space step, k = 0, ±1, ±2, . . . , ±(n/2 − 1), −n/2; i is the imaginary unit, F denotes the DFT and F−1 denotes the inverse DFT. Basically, the idea of the PSM is to approximate space derivatives by making use of the DFT
∂ mU = F−1 (ik)m F(U ) , m ∂X
(8)
reducing therefore PDE to ODE and then to use standard ODE solvers for integration with respect to time. The regular PSM algorithm is derived for ut = Φ (u, ux , u2x , . . . , umx ) type equations. In our case, however, we have also a mixed partial derivative term δβ UTTXX in the HE (5) and thus the standard PSM has to be modified [27–29]. Therefore we rewrite the HE (5) so that all partial derivatives with respect to time are in the left-hand side of the HE UTT − δβ UTTXX = bUXX +
µ 2
UX2 X
− δ γ UXX −
√ λ δ 2
,
2 UXX
(9)
XX
and introduce a new variable Φ = U − δβ UXX . After that, making use of properties of the DFT, one can express the variable U and its spatial derivatives in terms of the new variable Φ : −1
U =F
F(Φ ) 1 + δβ k2
,
m ∂ mU −1 (ik) F(Φ ) =F . ∂Xm 1 + δβ k2
(10)
1130
A. Salupere, K. Tamm / Wave Motion 50 (2013) 1127–1139
Finally, Eq. (5) can be rewritten in terms of the variable Φ
ΦTT = bUXX +
µ 2
UX2 X
− δ γ UXX −
√ λ δ 2
.
2 UXX
(11)
XX
In Eq. (11) all partial derivatives of U with respect to X are calculated in terms of Φ by using expression (10) and therefore one can apply the PSM for numerical integration of Eq. (11). The FSE (4) is reduced to the system of first-order differential equations which are solved by the standard PSM without any further modifications. The calculations are carried out with the Python package SciPy [30], using the FFTW library [31] for the DFT and the F2PY [32] generated Python interface to the ODEPACK Fortran code [33] for the ODE solver. 3. Material parameters and dispersion In order to reduce the dimension of the material parameters space the material parameters in the free energy expression (2) which are responsible for the linear part of the model are collected into two combined parameters γA2 and γ12 which have clear physical interpretation. 3.1. Combined parameters Parameters γA2 and γ12 combine the linear parts of the material parameters (see [8,24]) describing macro- and microstructure and interaction between those
γA2 = 1 − b =
D2 AB
,
γ12 = s2 =
γ ρC = . β AI
(12)
Parameter γA2 is related to the speed of long waves in the system and γ12 to the speed of short waves (Which feel the influence of the microstructure) in the system. The parameters describing nonlinearity are used separately. Combined parameter γA2 has a physical meaning in the domain from 0 up to 1. In the limiting cases, if γA2 = 0 then the HE (5) is reduced into the classical wave equation and if γA2 = 1 then the normalized speed of the long waves would be zero. 3.2. Dispersion Dispersive properties of a material provide insight into the internal structure and degrees of freedom in the system. We have the HE (5) which is derived from the FSE (4) through slaving principle. The HE has only a single dispersion curve (branch) while the dispersion curve of the FSE maintains two branches — so called acoustic and optical (higher-order) branches, where the former provides insight about scale of the underlying microstructure and the latter is related to the internal degrees of freedom (see [34–36]). The dispersion type for the HE (and for the acoustic dispersion branch of the FSE) can be determined by the sign of quantity (see [8] for details):
Γ = 1 − γA2 − γ12 .
(13)
If Γ is positive, we have the normal dispersion case, if negative, we have the anomalous dispersion case and if it is equal to zero, we have the dispersionless case. For the optical dispersion branch of the FSE the dispersion type is always anomalous. It has been shown by Metrikine in [37] that the higher order dispersion branch is essential in such models for achieving causality. In essence, without the higher order dispersion branch an infinite propagation speed for some disturbances in such models may appear. It should be noted that in the dispersionless case the acoustic dispersion branches of the FSE and HE coincide. For the sake of clarity — the classical way of dealing with dispersion types is to compare the group and phase speeds in the system. If these are equal to each other then we have dispersionless case. If the group speed is larger than the phase speed (cgr > cph ) then we have the anomalous dispersion case and if the group speed is smaller than the phase speed (cgr < cph ) then the normal dispersion case. In terms of wave propagation the essence of a dispersion is that the phase speed of the wave depends on the wave number of the harmonics. This is due to the length scale that the underlying microstructure brings into the models. In Fig. 1 a combined parameter domain and points of interest therein are presented. The dashed line represents parameter combinations corresponding to the dispersionless case, above that line the dispersion type is anomalous while below the dashed line the dispersion type is normal. The crosses, a ring and the triangles in Fig. 1 represent the investigated points in the combined parameter space for which the particular parameter combinations are presented later in the paper. The domain between thick solid lines corresponds to parameter combinations where concurrency between solutions of the HE and FSE is found to be good according to the dispersion analysis and numerical simulations [15,24]. In the area between the thick solid lines in Fig. 1 the difference between acoustical dispersion curves for the FSE and HE is 5% or less at dimensionless wavenumber of 1.5 (see Fig. 2 for examples of dispersion curves). In Fig. 2 some dispersion curves are presented for various points in the plane of combined parameters. Two upper subfigures give examples for two points in the horizontal line of crosses and the two lower subfigures give examples for
A. Salupere, K. Tamm / Wave Motion 50 (2013) 1127–1139
1131
Fig. 1. The combined parameter plane and points of interest corresponding to different values of parameters B, C and D.
end-points of the vertical line of triangles in Fig. 1. Top left and bottom right subfigures correspond to the normal dispersion case and top right and bottom left subfigures to these of the anomalous dispersion (with regard to the acoustic branch). 4. Results and discussion Material parameters are taken as A = 20,
B=
N = 50 ± 25,
D2 AγA2
,
C = 11 ± 5,
D = 10 ± 5, (14)
M = 5000 ± 2500,
ρ = 1,
I =
4 5
,
where parameter B is picked up to visit the points of interests defined by the parameter D in the combined parameter space i.e., parameter B has values of 60,
125 3
, 1500 , 375 , 500 , 15, 1500 , 125 , 1500 , 375 , 20 49 16 27 121 12 169 49 3
. Parameters C and D are investigated
in the given range with the step size 1 while parameter N is investigated with the step size 5 and parameter M with the step size 500. The points of interest corresponding to these parameter combinations in the combined parameter space can be seen in Fig. 1. In Fig. 1 the ‘reference case’ corresponds to the parameters combination against what other cases are compared to spot the effects which can caused by increasing or decreasing of material parameters. In the reference case B = 15, C = 11, D = 10, N = 50, M = 5000. The geometrical parameters are taken lo = 1,
Lo = 50,
Uo = 1,
Bo =
1 4
,
Ao = 1,
(15)
where lo is the characteristic scale of the microstructure, Lo — the characteristic wavelength of the initial excitation and Uo — the amplitude of the initial excitation while Bo — the parameter controlling the width of the initial pulse and Ao — the normalized amplitude of the initial pulse. The values above are taken for a certain hypotetic functionally graded material [15]. A reference case of formation, propagation and interaction of two solitary waves from a localized initial pulse in case of the FSE is shown in Fig. 3. The initial bell-shaped pulse splits into two solitary pulses propagating in the opposite directions which evolve differently in time — the pulse propagating to the left evolves towards a peakon-like shape [18,38,39] while the pulse propagating to the right maintains its bell-like shape which can be spotted also in the example solution in Fig. 3 at closer scrutiny. Such a asymmetry is caused by the internal structure of the model equations. In Fig. 4 one can see the phase plot corresponding to the waveprofiles at the end of integration (dashed and dotted lines) and the initial solitary wave
1132
A. Salupere, K. Tamm / Wave Motion 50 (2013) 1127–1139
Fig. 2. Dispersion curves for the FSE and HE at γA2 = 0.168, γ12 = 0.688 (top left), γA2 = 0.48, γ12 = 0.688 (top right), γA2 = 0.3333, γ12 = 1 (bottom left), and γA2 = 0.333, γ12 = 0.375 (bottom right). Variable η is the dimensionless frequency and variable ξ the dimensionless wavenumber.
Fig. 3. The typical solution of Eq. (4).
(solid line) — on horizontal axis there is a spatial derivative of the displacement (UX ) and on vertical axis the displacement amplitude of the wave (U ). A term ‘normalized initial wave’ is used for the waveprofile at a short time after the start where the initial pulse has divided into the two solitary pulses propagating in opposite directions but still having almost exactly sech2 -type profile, just with an amplitude of 0.5 instead of the initial amplitude of one. The main characteristics under investigation are the wave speeds which are tracked by finding the trajectory of the pulse peak and calculating the speed from these and the wave shape which is investigated mainly by making use of the phase-plots where even relatively small changes in the pulse shape are clearly visible.
A. Salupere, K. Tamm / Wave Motion 50 (2013) 1127–1139
1133
Fig. 4. The phase plot corresponding to the reference case at T = 480. Solid line — normalized initial wave, (blue) dotted line — peakon-like pulse, (red) dashed line — bell-like pulse.
Table 1 Calculated average speeds cB and cD and characteristic speed co against parameters B and D. D
cD
B
cB
γA2
co
cD - co
5 6 7 8 9 10 11 12 13 14 15
0.960 0.939 0.916 0.888 0.857 0.819 0.775 0.725 0.662 0.593 0.502
60 125/3 1500/49 375/16 500/27 15 1500/121 125/12 1500/169 375/49 20/3
0.960 0.939 0.917 0.888 0.856 0.818 0.775 0.723 0.662 0.590 0.501
0.083 0.120 0.163 0.213 0.270 0.333 0.403 0.480 0.563 0.653 0.750
0.917 0.880 0.837 0.787 0.730 0.667 0.597 0.520 0.437 0.347 0.250
0.043 0.059 0.080 0.102 0.127 0.152 0.179 0.205 0.226 0.246 0.252
4.1. Parameters B and D The average speeds in Table 1 have been calculated for the pulse propagating to the right in the interval from separating from the initial pulse until the first interaction for the FSE. The speeds for the HE are practically the same and are within ±0.02 which is the accuracy of the speeds in the table. In Table 1, cB denotes the speed calculated for the cases where parameter B was used to control the value of the combined parameter γA2 and cD denotes the speed calculated for the cases where parameter D was used to control the value of the combined parameter γA2 . The speed co is the characteristic speed for √ the bulk medium and can be found as co = A/ρ . The speed differences between the solutions of the FSE (4) and HE (5) as well as between speeds corresponding to the parameters B and D resulting in the same combined parameter γA2 values are negligible. What can be noted, however, is that with the increasing parameter γA2 the gap between the speeds cB , cD and the characteristic speed co increases. While the parameters D and B have similar influence on the wave propagation speed their effect on the pulse shape is slightly different. At lower values of the combined parameter γA2 the differences are hard to spot but at highest value of D or lowest value of B (D = 15 or B = 20/3) one can see from the phase plots in Fig. 5 that in the case of high D the oscillatory structure forming in the peakon like profile has higher frequency oscillations than in the case of low values of parameter B. Regardless of the observed differences qualitatively the oscillatory structure is the same as for both high D and low B values altogether 5 oscillations can be distinguished in Fig. 5 on the peakon-like pulse. No notable differences can be detected in the bell-like pulse for the high values of parameter D or low values of parameter B.
1134
A. Salupere, K. Tamm / Wave Motion 50 (2013) 1127–1139
84 83 82 81 80 79 78 77 76 75 74 47π
T
156 155 154 153 152 151 150 149 148 147 146 47π
B = 15; D = 15 B = 20/3; D = 10 47.5π
48π
48.5π
49π
T
T
T
, D = 10, T = 480 (top right); B = 15, D = 5, T = 460 (bottom left); B = 15, Fig. 5. Phase plots for B = 60, D = 10, T = 460 (top left); B = 20 3 D = 15, T = 480 (bottom right). Solid line — normalized initial wave, (blue) dotted line — peakon-like pulse, (red) dashed line — bell-like pulse, and parameters C , M and N have reference values.
B = 15; D =5 B = 60; D = 10 47.5π
48π 47π ≤ X < 49π
48.5π
49π
307 306 305 304 303 302 301 300 299 298 297 23π 163 162 161 160 159 158 157 156 155 154 153 23π
B = 15; D = 15 B = 20/3; D = 10 23.5π
24π
25.5π
25π
B = 15; D =5 B = 60; D = 10 23.5π
24π 23π ≤ X < 25π
25.5π
25π
Fig. 6. Trajectories of the first (left column) and the second (right column) interaction of solitary waves for different parameter B and D values. Top row — γA2 = 0.75, bottom row — γA2 = 0.083.
In order to clarify the influence of parameters B and D on the behaviour of the solution during interactions we trace wave-trajectories for γA2 = 0.083 and γA2 = 0.75 (see Fig. 6). According to Eq. (12) these values can be reached by different combinations of parameters B and D. For both values of γA2 two combinations of B and D are considered. It is clear that for smaller values of D compared trajectories practically coincide (lower subfigures), but for higher values of D small differences
A. Salupere, K. Tamm / Wave Motion 50 (2013) 1127–1139
1135
Fig. 7. Phase plots for C = 6 (left) and C = 16 (right) at T = 480. Solid line — normalized initial wave, (blue) dotted line — peakon-like pulse, (red) dashed line — bell-like pulse, and parameters B, D, M and N have reference values.
can be detected (higher subfigures). It should be mentioned that in here the solitary pulses do not experience phase shifts during the interactions due to the normalization procedure used which removes the usual speed–amplitude dependence which exists in the Boussinesq- and KdV-type equations (see, for example [12,40]). 4.2. Parameter C The effect of parameter C on the calculated average speeds is negligible and it has no influence on the observed interaction trajectories of the solitary pulses. In Fig. 7 one can see the two phase-plots representing waveprofile shape at the end of integration at T = 480 corresponding to the different values of the parameter C . The most notable effect on the phase plots is that the oscillatory structure which develops by the end of the integration interval switched from the ‘front’ (in the direction of propagation, anomalous dispersion) to the ‘back’ (opposite to the direction of propagation, normal dispersion) of the sharper pulse. The closer the point is in the γA2 − γ12 plane to the dispersionless case the sharper the peakonlike pulse. Strictly speaking a peakon should have a discontinuity at its peak but our numerical scheme does not handle specifically discontinuities so the integration interval has been picked short enough to avoid reaching a true discontinuity at the sharper profile peak in the dispersionless case. Parameter C in the present study is used to control the values of combined parameter γ12 which is interpreted as a speed of short waves which feel the influence of the microstructure. Under the used parameter combinations, this translates into controlling the dispersion type for the acoustic dispersion branch of the governing equations. More specifically, looking at Eqs. (12) and (13) it is clear that parameter C has a significant impact on the speed of harmonics with a shorter wavelength and dispersion type which would explain the observed differences in the waveprofile shapes at different parameter C values. 4.3. Nonlinear parameters N and M The effect of parameters N and M on the calculated average speeds is negligible. In Fig. 8 the highest presented parameter N value is N = 55 as at higher parameter N values the wave starts breaking because of the peakon-like profile peak gets too sharp as to interfere with the numerical scheme used. At the highest presented value of parameter N one can already see emerging high frequency oscillations from the discontinuity at the waveprofile peak. The lower the parameter N the smaller the difference between waveprofiles propagating in the opposite directions and vice versa — the higher the parameter N the greater the difference (Fig. 8) — the sharper pulse gets thinner and sharper while the smoother pulse gets even thicker with increasing N. Overall the parameter M influences on the waveprofile evolution under the used parameter combinations can be considered weak (see Fig. 8). However, one can note that even at the lowest considered parameter M value the beginnings of the oscillatory structure can be spotted in the phase-plot corresponding to the peakon-like pulse. Compared to the parameter N the influence of parameter M on the waveprofile shape is negligible under the investigated parameter combinations. It is known from previous studies (see, for example [41] where its shown for the HE (5)) that under some parameter combinations the parameter M can affect the asymmetry of a waveprofile so that it ‘tilts’ in the direction of propagation or opposite of the direction of propagation (according to on dispersion type for such parameter combination). However, it must be added that in here the conditions described in [41] are not fulfilled. A change in the sharper pulse shape at highest investigated value of parameter M can be noted in Fig. 8 where the emerging oscillations appear to be slightly stronger than in the case of weakest considered parameter M (Fig. 8) or the reference case (Fig. 4).
1136
A. Salupere, K. Tamm / Wave Motion 50 (2013) 1127–1139
Fig. 8. Phase plots for M = 5000, N = 25 (top left), M = 5000, N = 55 (top right) and M = 2500, N = 50 (bottom left) and M = 7500, N = 50 (bottom right) at T = 480. Solid line — normalized initial wave, (blue) dotted line — peakon-like pulse, (red) dashed line — bell-like pulse, and parameters B, C and D have reference values.
Another notable difference between the influence of nonlinear parameters on the waveprofile evolution emerges when we investigate the apparent trajectories of the waveprofiles by tracking their maxima. In the case of parameter N it is possible to observe a clear difference between trajectories drawn by the waveprofile peak during interactions at different parameter N values (Fig. 9). Outside the interactions the values of parameter N have no visible effect on the trajectories. At closer inspection (see Fig. 10) an interesting phenomenon appears: the pulse itself is practically stationary during the interactions, but the peak of the waveprofile moves. At the beginning of the interaction a combined waveprofile is asymmetric and tilted to the left. During the interaction it reaches symmetric shape at the peak of interaction and is tilted to the right after the peak of interaction until the waveprofiles separate again after the interaction. In addition, in Fig. 10 it is possible to spot that during the fourth interaction the waveprofile has sharper peak than during the second interaction. Of course, it should be useful to perform rigorous analysis of interactions for the present model, like it is done for the KdV equation [42]. 4.4. Concurrency of the solutions In order to calculate the difference between the solutions of the HE and FSE the quantity ∆S is introduced as
∆S =
n 1
n i =1
∆i ,
∆i = U HE (Xi , Tf ) − U FSE (Xi , Tf ) ,
(16)
where n is the number of grid points. In the sense of the quantity ∆S for the parameters C , N and M the concurrency of the solutions of the HE (5) and FSE (4) does not depend on the parameter values, however, for the parameters B and D the parameter value also affects the concurrency of the solutions. The concurrency of the solutions gets worse with decreasing
190 189 188 187 186 185 184 183 182 181 180 179 22π
T
Ν=25 Ν=50 Ν=75
23π
190 189 188 187 186 185 184 183 182 181 180 179 22π
24π
25π
26π
Μ=2500 Μ=5000 Μ=7500
T
T
T
A. Salupere, K. Tamm / Wave Motion 50 (2013) 1127–1139
23π
24π
25π
26π
374 373 372 371 370 369 368 367 366 365 364 22π
1137
Ν=25 Ν=50 Ν=75
23π
24π
25π
374 373 372 371 370 369 368 367 366 365 364
26π
Μ=2500 Μ=5000 Μ=7500
22π
23π
24π
25π
26π
1.0 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 18π
30π
1.0 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 18π
30π
1.0 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 18π
U, T = 183.20 U, T = 184.00 U, T = 184.80 U, T = 185.60 U, T = 186.40
U
1.0 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 18π
20π
22π
24π
26π
28π
U, T = 183.20 U, T = 184.00 U, T = 184.80 U, T = 185.60 U, T = 186.40
U
U
U
Fig. 9. Trajectories of the second (left) and the fourth (right) interactions for the nonlinear parameter values N = 25, 50, 75 and M = 5000 (top) and M = 2500, 5000, 7500 and N = 50 (bottom). Parameters B, C and D have reference values.
20π
22π
24π X
26π
28π
U, T = 367.80 U, T = 368.60 U, T = 369.40 U, T = 370.20 U, T = 371.00
20π
22π
24π
26π
28π
30π
24π X
26π
28π
30π
U, T = 367.80 U, T = 368.60 U, T = 369.40 U, T = 370.20 U, T = 371.00
20π
22π
Fig. 10. The second (left) and the fourth (right) interactions for the nonlinear parameter values N = 75, M = 5000 (top) and N = 50, M = 7500 (bottom). Parameters B, C and D have reference values.
of the parameter B or increasing the parameter D. Under the investigated parameter combinations the quantity log10 (∆S ) is between −5 and −7 which we consider to be a good concurrency between the solutions of the FSE (4) and HE (5) because this accuracy is sufficient to investigate in detail do the solutions of the FSE and HE behave in a similar way or not. If quantity log10 (∆S ) is close to −3 then the concurrency of the solutions cannot be considered anymore satisfactory as there can be already substantial differences between the solutions of the FSE (4) and HE (5) (see [15] for details) after sufficiently short evolution in time.
1138
A. Salupere, K. Tamm / Wave Motion 50 (2013) 1127–1139
5. Conclusions In the present paper we studied wave propagation in solids having a Mindlin-type microstructure. The considered free energy function (2) includes material parameters of different kind. Parameter A characterizes linear properties of the material at the macroscale and parameters B and C characterize these at the microscale. Interaction between macro- and microscale effects is described by the parameter D. In their turn, parameters M and N characterize the influence of the nonlinear properties of the micro- and macroscale. We were particularly interested how parameters B, C , D, M and N influence the propagation and evolution of solitary waves. To summarize our findings following conclusions can be drawn: (i) The speed of emerged solitary waves is directly related to values of parameters B and D — the higher the value of B or the lower value of D the higher the speed of waves. The influence of parameters C , M and N on the wave speed can be neglected. (ii) Parameters C and M influence the character of asymmetry of solitary pulses (including the location of oscillatory structures with respect to the direction of propagation of main pulses). In the considered domain of parameters the influence of M is weak, but according to earlier studies [17,41] the character of asymmetry of the travelling wave solution of the HE depends directly on the value of M. (iii) Parameters D and N influence the shape of single pulses in the following way: the greater the parameter D or N the sharper the peakon-like pulse and the wider the bell-like pulse. Parameter B has the opposite effect on the pulse-shapes. (iv) Parameter N influences the rate of the change of asymmetry during the interaction process — the higher the value of N the higher the rate of the change of asymmetry (i.e. the faster the peak of the waveprofile moves during the interaction). It influences also the trajectories within the interaction process. Last but not least, if the propagation of solitary pulses in nonlinear dispersive media is studied then the question about solitonic character of the wave structure appears. Notwithstanding that the interactions between solitary pulses are not completely elastic (the peakon-like pulse gets sharper and the bell-like pulse gets wider and oscillatory structures can be generated due to interactions), the main pulses propagate at a constant speed and almost restore their shape and speed after interactions. Therefore we conclude that the emerged wave structures modelled by a Boussinesq-type equation (5) or its full form (4) have solitonic behaviour (at least throughout several first interactions). Acknowledgements The research was supported by the Estonian Science Foundation Grant No. 8658, the Estonian Ministry of Education and Research (block grant SF0140077s08) and by the European Union through the Regional Development Fund (projects CENS and SCDM). The authors thank Prof. J. Engelbrecht for helpful discussions and suggestions related to this paper. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27]
M. Born, T. von Kármán, Uber Schwingungen in Raumgittern, Phys. Z. 13 (1) (1912) 297–309. G.A. Maugin, Nonlinear Waves in Elastic Crystals, Oxford Univercity Press, Oxford, 1999. R. Mindlin, Micro-structure in linear elasticity, Arch. Ration. Mech. Anal. 16 (1) (1964) 51–78. A. Eringen, Linear theory of micropolar elasticity, J. Math. Mech. 15 (1966) 909–923. A. Metrikine, H. Askes, One-dimensional dynamically consistent gradient elasticity models derived from a discrete microstructure, part 1: generic formulation, Eur. J. Mech. A Solids 21 (2002) 555–572. V. Erofeev, Wave Processes in Solids with Microstructure, World Scientific, Singapore, 2003. S. Papargyri-Beskou, D. Polyzos, D. Beskos, Wave dispersion in gradient elastic solids and structures: a unified treatment, Int. J. Solids Struct. 46 (21) (2009) 3751–3759. J. Engelbrecht, A. Berezovski, F. Pastrone, M. Braun, Waves in microstructured materials and dispersion, Phil. Mag. 85 (33–35) (2005) 4127–4141. A. Berezovski, J. Engelbrecht, G.A. Maugin, Generalized thermomechanics with dual internal variables, Arch. Appl. Mech. 81 (2) (2010) 229–240. H. Askes, E.C. Aifantis, Gradient elasticity theories in statics and dynamics — a unification of approaches, Int. J. Fract. 139 (2) (2006) 297–304. A. Berezovski, J. Engelbrecht, M. Berezovski, Waves in microstructured solids: a unified viewpoint of modeling, Acta Mech. 220 (1–4) (2011) 349–363. C.I. Christov, G.A. Maugin, A.V. Porubov, On Boussinesq’s paradigm in nonlinear wave propagation, C. R. Mecanique 335 (9–10) (2007) 521–535. C.I. Christov, M. Velarde, Inelastic interaction of Boussinesq solitons, Internat. J. Bifur. Chaos Appl. Sci. Engrg. 4 (5) (1994) 1095–1112. C.I. Christov, G.A. Maugin, M. Velarde, Well-posed Boussinesq paradigm with purely spatial higher-order derivatives, Phys. Rev. E 54 (4) (1996) 3621–3638. K. Tamm, Wave propagation and interaction in mindlin-type microstructured solids: numerical simulation, Ph.D. Thesis, Tallinn University of Technology, 2011. J. Engelbrecht, A. Salupere, K. Tamm, Waves in microstructured solids and the Boussinesq paradigm, Wave Motion 48 (2011) 717–726. J. Janno, J. Engelbrecht, Solitary waves in nonlinear microstructured materials, J. Phys. A: Math. Gen. 38 (2005) 5159–5172. K. Tamm, A. Salupere, On the propagation of 1D solitary waves in Mindlin-type microstructured solids, Math. Comput. Simul. 82 (7) (2012) 1308–1320. A. Sergeeva, E. Pelinovsky, T. Talipova, Nonlinear random wave field in shallow water: variable Korteweg–de Vries framework, Nat. Hazards Earth Syst. Sci. 11 (2) (2011) 323–330. A.V. Porubov, G.A. Maugin, B.R. Andrievsky, Solitary wave interactions and reshaping in coupled systems, Wave Motion 48 (8) (2011) 773–781. A.V. Porubov, I.V. Andrianov, V.V. Danishevskyy, Nonlinear strain wave localization in periodic composites, Int. J. Solids Struct. 49 (2012) 3381–3387. C. Sun, J.D. Achenbach, G. Herrmann, Continuum theory for a laminated medium, J. Appl. Mech. 35 (3) (1968) 467. J. Engelbrecht, F. Pastrone, Waves in microstructured solids with nonlinearities in microscale, Proc. Est. Acad. Sci., Phys. Math. 52 (1) (2003) 12–21. T. Peets, M. Randrüüt, J. Engelbrecht, On modelling dispersion in microstructured solids, Wave Motion 45 (4) (2008) 471–480. G. Whitham, Linear and Nonlinear Waves, Wiley-Interscience, New York, 1974. B. Fornberg, A Practical Guide to Pseudospectral Methods, Cambridge University Press, Cambridge, 1998. A. Salupere, The pseudospectral method and discrete spectral analysis, in: E. Quak, T. Soomere (Eds.), Applied Wave Mathematics, Springer, Berlin, 2009, pp. 301–334.
A. Salupere, K. Tamm / Wave Motion 50 (2013) 1127–1139
1139
[28] L. Ilison, A. Salupere, Propagation of sech2 -type solitary waves in hierarchical KdV-type systems, Math. Comput. Simul. 79 (11) (2009) 3314–3327. [29] L. Ilison, A. Salupere, P. Peterson, On the propagation of localized perturbations in media with microstructure, Proc. Est. Acad. Sci., Phys. Math. 56 (2) (2007) 84–92. [30] E. Jones, T. Oliphant, P. Peterson, SciPy: open source scientific tools for Python, 2007. http://www.scipy.org. [31] M. Frigo, S. Johnson, The design and implementation of FFTW 3, Proc. IEEE 93 (2) (2005) 216–231. [32] P. Peterson, F2PY: Fortran to Python interface generator, 2005. http://cens.ioc.ee/projects/f2py2e/. [33] A. Hindmarsh, ODEPACK, A Systematized Collection of ODE Solvers, Vol. 1, North-Holland, Amsterdam, 1983. [34] L. Brillouin, Wave Propagation in Periodic Structures, Dover Publications Inc., New York, 1953. [35] A. Eringen, Microcontinuum Field Theories I: Foundations and Solids, Springer, New York, 1999. [36] K. Tamm, T. Peets, On the influence of internal degrees of freedom on dispersion in microstructured solids, Mech. Res. Commun. 47 (2013) 106–111. [37] A. Metrikine, On causality of the gradient elasticity models, J. Sound Vib. 297 (3–5) (2006) 727–742. [38] R. Camassa, D. Holm, An integrable shallow water equation with peaked solitons, Phys. Rev. Lett. 71 (11) (1993) 1661–1664. [39] H. Kalisch, J. Lenells, Numerical study of traveling-wave solutions for the Camassa–Holm equation, Chaos Solitons Fractals 25 (2) (2005) 287–298. [40] A. Salupere, On hidden solitons in KdV related systems, Math. Comput. Simulat. (2013) submitted for publication. [41] J. Janno, J. Engelbrecht, An inverse solitary wave problem related to microstructured materials, Inverse Problems 21 (2005) 2019–2034. [42] A. Salupere, J. Engelbrecht, P. Peterson, Long-time behaviour of soliton ensembles, part I—emergence of ensembles, Chaos Solitons Fractals 14 (2002) 1413–1424.