Effective boundary condition for the reflection of shear waves at the periodic rough boundary of an elastic body

le this surface is independent of the choice of the origine in Fig. 2, both ym1 and (S − a/h) do depend on this choice Taking the limit ym1 → +∞ in the above equation, and using the matching condition (11), we finally get σ11 (0, x2, τ) = ( S − a h ) ∂σ01 ∂x1 (0, x2, τ) + C ∂σ 0 2 ∂x2 (0, x2, τ) , (20) where we have defined (with Ω = lim ym1→+∞ Ωm) C ≡ − ∫ Ω dy ∂W ∂y2 . (21) Effective boundary condition for the reflection of shear waves at the periodic rough boundary of an elastic body 311 2.4. The final homogenized problem The equations in the bulk (7) and the associated boundary conditions (13) and (20) could be used to solve the homogenized problem iteratively: first compute ( v0,σ0 ) for a flat stress-free boundary (compute also C in (21)) and use the results to get the right hand- side term in (20); then, compute ( v1,σ1 ) ; finally, recompose v0 + εv1 which approximate vε up to O ( ε2 ) . As discussed in [3], it is preferable to handle a unique problem and this is done by defining the fields ( vh, σh ) satisfying the following homogenized problem ∂σh ∂τ = ∇xvh, ∂v h ∂τ = divx σh, for x1 > 0, σh1 (0, x2, τ) = ε ( S − a h ) ∂σh1 ∂x1 (0, x2, τ) + εC ∂σ h 2 ∂x2 (0, x2, τ). (22) A quick analysis of (7), (13) and (20) shows that vh has the same expansion as v0 + v1, and thus as vε, up to O ( ε2 ) (the same for σh). Finally, coming back to the real space (through (2)), we get  ∂Σh ∂t = µ∇Vh, ρ∂V h ∂t = divΣh, in Dh, Σh1 (0, X2, t) = (eϕ− a) ∂Σh1 ∂X1 (0, X2, t) + hC ∂Σ h 2 ∂X2 (0, X2, t), the appropriate radiation condition at X1 → +∞. (23) The above homogenized problem depends on the choice of a because of the boundary condition (in an obvious manner; note that C does not depend on a) but also because of the domain Dh (Fig. 1). In the following section, we address the problem of the choice of a. 3. DEPENDANCE OF THE HOMOGENIZED PROBLEM ON THE α-VALUE 3.1. Energy conservation in the homogenized problem 3.1.1. The boundary energy in the equation of energy conservation The choice of the origine y1 = 0 affects the boundary condition in (23) and the posi- tion of the equivalent boundary (notably it is within the actual roughnesses for−e < a < e). In the original problem, the elastic energy is E = ∫ D dV [ ρ 2 V2 + 1 2µ |Σ|2 ] , and in the homogenized problem, it has to be of form Eh = ∫ Dh dV [ ρ 2 Vh2 + 1 2µ ∣∣∣Σh∣∣∣2]+ Eb, where Eb is the energy of the equivalent boundary. In the absence of losses, the equations of energy conservation read dE/dt = 0 = dEh/dt in the two problems. From (23), it is 312 Agne`s Maurel, Jean-Jacques Marigo, Kim Pham easy to see that the energy conservation in the homogenized problems reads d dt ∫ Dh dV [ ρ 2 Vh 2 + 1 2µ ∣∣∣Σh∣∣∣2]+ ∫ X1=0 dX2VhΣh1 + ∫ Lh dSΠ.n = 0, where Lh = ∂Dh\ (X1 = 0) is the boundary of Dh except the segment at X1 = 0. In view of the stability of numerical scheme, we want I ≡ ∫ X1=0 dX2VhΣh1 I ≡ ∫ X1=0 dX2VhΣh1 to be the time derivative of a positive energy Eb. Applying the boundary conditions of Eqs. (23), we get I = ∫ X1=0 dX2 [ (eϕ− a)∂Σ h 1 ∂X1 Vh + hC ∂Σ h 2 ∂X2 Vh ] , leading to I = dEb/dt (with (23) again), withEb = 1 2 ∫ X1=0 dX2 [ ρ(eϕ− a)Vh2 + hCaΣ 2 2 µ ] , hCa ≡ eϕ− a− hC. (24) Note that the integration by part of Vh∂X2Σ h 2 makes a boundary term to appear, [ Σh2V h ] X2 and something should be said at both extremities of the equivalent surface; this is disre- garded in the present paper. To ensure Eb ≥ 0, one has to choose a ≤ eϕ and to ensure Ca ≥ 0. In fact, this latter constraint is more restrictive since C ≥ 0. Indeed, coming back to the problem (17) satisfied by W, we have 0 = ∫ Ω dyW∆W = ∫ Ω dyW div (∇W + e2) = − ∫ Ω dy (∇W + e2)∇W + ∫ ∂Ω dl∇ (W + y2) .nW, (25) and the integral on ∂Ω vanishes because of the boundary conditions in (17). It results that C = ∫ Ω dy|∇W|2 ≥ 0. (26) In the following, we inspect the bound on a to ensure that hC ≤ eϕ− a which guaranties Eb ≥ 0 in (24). 3.1.2. Bounds for the a value ensuring a positive boundary energy Eb To find a bound for C (and thus on a), it is sufficient to remark that the homoge- nized problem (17) admits a variational formulation from which a principle of energy minimization can be associated, specifically E(S) ≤ E(S˜), with E(S˜) ≡ ∫ Ω ∣∣S˜− e2∣∣2 , (27) for any S˜ admissible field satisfying S˜.n|Γ = 0 and S˜→ e2 when y1 → +∞, and E(S) = ∫ Ω |∇W|2, with S = ∇ (W + y2) . (28) Effective boundary condition for the reflection of shear waves at the periodic rough boundary of an elastic body 313 Using (26), we must have C = E(S) ≤ E(S˜). (29) 9 Note that the integration by part of V h∂X2Σ h 2 makes a boundary term to appear, [Σ h 2V h]X2 and something should be said at both extremities of the equivalent surface; this is disregarded in the present paper. To ensure Eb ≥ 0, one has to choose a ≤ eϕ and to ensure Ca ≥ 0. In fact, this latter constraint is more restrictive since C ≥ 0. Indeed, coming back to the problem (17) satisfied by W , we have 0 = ∫ Ω dy W∆W = ∫ Ω dy Wdiv(∇W + e2) = − ∫ Ω dy (∇W + e2)∇W + ∫ ∂Ω dl∇(W + y2).n W (25) and the integral on ∂Ω vanishes because of the boundary conditions in (17). It results that C = ∫ Ω dy |∇W |2 ≥ 0. (26) In the following, we inspect the bound on a to ensure that hC ≤ eϕ− a which guaranties Eb ≥ 0 in (24). 2. Bounds for the a value ensuring a positive boundary energy Eb To find a bound for C (and thus on a), it is sufficient to remark that the homogenized problem (17) admits a variational formulation from which a principle of energy minimization can be associated, specifically E(S) ≤ E(S˜), with E(S˜) ≡ ∫ Ω |S˜− e2|2, (27) for any S˜ admissible field satisfying S˜.n|Γ = 0 and S˜→ e2 when y1 → +∞, and E(S) = ∫ Ω |∇W |2, ith S =∇(W + y2). (28) Using (26), we must have C = E(S) ≤ E(S˜). (29) S˜ S y1 0 1 FIG. 4. The chosen admissible field S˜; S˜ = 0 in the roughness and S˜ = e2 outside. We choose S˜ piecewise constant along y1: for the y1 values inside the roughness, S˜ = 0 and outside S˜ = e2 (Fig. 4). For this field, we have E(S˜) ≡ ∫ Ω |S˜− e2|2 = S, (30) from which C ≤ S imposes the final (sufficient but not necessary) condition Eb ≥ 0, if a ≤ 0, (31) and this holds for any shape of the roughnesses. Obviously, a less strict criterion can be found for particular a shape. Fig. 5 reports the variation of Ca varying a and ϕ for rectangular roughness shape as we shall consider in the forthcoming section. It is visible that increasing ϕ produce larger values of a allowing for Ca ≥ 0; from the estimate C ' ϕe/h− piϕ2/16 given in [5], we get Ca ' −a/h+ piϕ2/16, and the line a/h = piϕ2/16 reasonably corresponds to Ca = 0 (dotted black line in Fig. 5). Fig. 4. The chosen admissible field S˜; S˜ 0 in the roughness and S˜ e2 outside We choose S˜ piecewise constant along y1: for the y1 values inside the roughness, S˜ = 0 and outside S˜ = e2 (Fig. 4). For this field, we have S˜ = e2, (30) from which C ≤ S imposes the final (sufficient but not necessary) condition Eb ≥ 0, if a ≤ 0, (31) and this holds for any shape of the roughnesses. Obviously, a less strict criterion can be und for particular a shape. Fig. 5 reports the variation of Ca varyi g a and ϕ f rectangular roughness shape as we sh ll consider in the forthcoming ection. It is visible that increasing ϕ produce larger values of a allowing for Ca ≥ 0; from the estimate C ' ϕe/h− piϕ2/16 given in [5], we get Ca ' −a/h + piϕ2/16, and the line a/h = piϕ2/16 reasonably corresponds to Ca = 0 (dotted black line in Fig. 5). 3.2. Accuracy of the homogenized solution with respect to the actual solution In this section, we inspect the accuracy of the homogenization and to do so, we shall consider the particular scattering problem of the reflection of a wave hitting the rough- nesses at oblique incidence θ. We shall work with complex fields (owing to the physical fields are the real parts of the computed complex ones); in the harmonic regime, the com- plex fields (and we shall consider the displacement field U(X)) have a time dependance in e−ikct and it will be omitted in the following. For an incident wave of the form Uinc(X) = e−ik cos θX1+ik sin θX2 , we discriminate the right-going wave, with X1 dependence in eik cos θX1 (corresponding to increasing X1 values for increasing time) and left-going waves with X1-dependence in e−ik cos θX1 (the incident wave is a left-going wave coming from X1 → +∞). In the actual problem, this is done in the configuration of Fig. 6(a). The wavefield for X1 > 0 reads U(X) = [ e−ik cos θX1 + Reik cos θX1 ] eik sin θX2 +Uev(X). (32) 314 Agne`s Maurel, Jean-Jacques Marigo, Kim Pham 10 0.5 0.5 0 ' 0 1 e/h e/h a/h 0 FIG. 5. Ca in (24) (with (21)) in colorscale as a function of a/h and ϕ for rectangular shape roughness (and e/h = 0.5); the dotted black line shows a/h = piϕ2/16 where Ca ' 0. B. Accuracy of the homogenized solution with respect to the actual solution In this section, we inspect the accuracy of the homogenization and to do so, we shall consider the particular scattering problem of the reflection of a wave hitting the roughnesses at oblique incidence θ. We shall work with complex fields (owing to the physical fields are the real parts of the computed complex ones); in the harmonic regime, the complex fields (and we shall consider the displacement field U(X)) have a time dependance in e−ikct and it will be omitted in the following. For an incident wave of the form U inc(X) = e−ik cos θX1+ik sin θX2 , we discriminate the right-going wave, with X1 dependence in e ik cos θX1 (corresponding to increasing X1 values for increasing time) and left-going waves with X1-dependence in e −ik cos θX1 (the incident wave is a left-going wave coming from X1 → +∞). In the actual problem, this is done in the configuration of Fig. 6(a). The wavefield for X1 > 0 reads U(X) = [ e−ik cos θX1 +Reik cos θX1 ] eik sin θX2 + U ev(X), (32) The boundary condition on the roughnesses are Neumann boundary condition∇U.n = 0, with n the local normal to the roughnesses. U ev(X) is an evanescent field excited in the vicinity of the roughnesses and vanishing at X1 → +∞. It is worth noting that U ev(X) = 0 for ϕ = 0 and ϕ = 1, since the problem is reduced to a one-dimensional problem along X1. With n = e1 in these cases, the boundary conditions are ∂X1U(0, X2) = 0 (ϕ = 0) and ∂X1U(−e,X2) = 0 (ϕ = 1); it is easy to see that exact solutions are R = 1 and R = e2ike cos θ respectively (and U ev(X) = 0 in both cases). The reflection coefficient, with |R| = 1 by conservation of the energy, is characterized by its phase. As written in (32), the phase of the reflection is defined by R = U r|Σ U inc|Σ , (33) with U r = (U − U inc) and where we choose the plane Σ corresponding to boundary between the roughnesses and the plain elastic body, Fig. 6. The homogenized problem (23) can be translated in terms of a complex displacement field Uh. With V h → −iωUh, we get Σhi → µ∂XiUh, where the arrow means going toward complex fields. (23) can be written in terms of the complex field of displacement Uh only, namely ∆Uh + k2Uh = 0, for X1 > 0, ∂Uh ∂X1 (0, X2) = (eϕ− a) ∂ 2Uh ∂X21 (0, X2) + hC ∂ 2Uh ∂X22 (0, X2), lim X1→+∞ [ ∂Uh,r ∂X1 − ik cos θUh,r → 0 ] , with Uh,r ≡ Uh − U inc. (34) Fig. 5. Ca in (24) (with (21)) in colorscale as a function of a/h and ϕ for ect ngular s ape rough- ness (and e/h = 0.5); the dotted black line shows a/h = piϕ2/16 where Ca ' 0 The boundary condition on the roughnesses are Neumann boundary condition ∇U.n = 0, with n the local normal to the roughnesses. Uev(X) is an evanescent field ex- cited in the vicinity of the roughnesses and vanishing at X1 → +∞. It is worth noting that Uev(X) = 0 for ϕ = 0 and ϕ = 1, since the problem is reduced to a one-dimensional prob- lem along X1. With n = e1 in these cases, the boundary conditions are ∂X1U (0, X2) = 0 (ϕ = 0) and ∂X1U (−e, X2) = 0 (ϕ = 1); it is easy to see that exact solutions are R = 1 and R = e2ike cos θ respectively (and Uev(X) = 0 in both cases). The reflection coefficient, with —|R| = 1 by conservation of the energy, is character- ized by its phase. As written in (32), the phase of the reflection is defined by R = Ur|Σ U inc|Σ , (33) with Ur = ( U − inc) and where we choose the plane Σ corresponding to boundary between the roughnesses and the plain elastic body, Fig. 6. The homogenized problem (23) can be translated in terms of a complex displacement field Uh. With Vh → −iωUh, we get Σhi → µ∂Xi Uh, where the arrow means going toward complex fields. (23) can be writte in terms of the complex field of displacement Uh only, namely  ∆Uh + k2Uh = 0, for X1 > 0, ∂Uh ∂X1 (0, X2) = (eϕ− a)∂ 2Uh ∂X21 (0, X2) + hC ∂ 2Ub ∂X22 (0, X2) , lim X1→+∞ [ ∂Uh,r ∂X1 − ik cos θUh,r → 0 ] , with Uh,r ≡ Uh −Uinc. (34) Effective boundary condition for the reflection of shear waves at the periodic rough boundary of an elastic body 315 11 . . . .. . .. . X2 X1 ⌃ U r 0 e 'h U inc ✓ . . . .. . .. . ⌃ Uh,r X2 X1a 0 U inc ✓ (a) (b) FIG. 6. Reflection of a wave U inc at oblique incidence on the rough boundary of the elastic body; (a) in the real problem with U r the reflected wave and (b) in the homogenized problem with Uh,r the reflected wave (the system of coordinate has been shifted of a, here a < 0; light grey region show the actual roughnesses, but they do not exist in the homogenized problem). In the above system, the origine of the X1-axis has been shifted of a (Fig. 6); to keep the same definition of the reflection coefficient, we want Rh = Uh,r|Σ U inc|Σ , (35) and now Σ is the plane X1 = a. An exact solution of (23) can be found, accounting for (35), of the form Uh(X) = [ e−ik cos θ(X1−a) +Rheik cos θ(X1−a) ] eik sin θX2 , (36) and applying the boundary condition of (23) at X1 = 0, we have Rh = 1 + ikh cos θ za(θ) 1− ikh cos θ za(θ) e 2ika cos θ, za(θ) ≡ ( ϕ e h − a h ) + C tan2 θ. (37) To begin with, we report in Fig. 7 an example of the wavefield in the actual problem together with its homogenized counterparts obtained for a = −e and a = 0 (the configuration is described in the figure caption). The agreement |U − Uh|/|U | (in L2 norm) is of 20% for a = −e and of 5% for a = 0, which suggests that small a value is preferable. In fact, something is not very convenient in the form of (37). We have found that at leading order, the roughnesses behave as a flat free boundary, and as previously said, this will be exact in the two limiting cases ϕ = 0 (leading to R = 1) and ϕ = 1 (leading to R = e2ike cos θ). In both cases, the elementary problem (17) has a trivial solution, W = constant, since Γ is a boundary y1 = constant, whence C = 0. But neither e nor a being zero, Rh does not go to the expected values; strictly, this is possible since the homogenized solution is an approximation of the real one, but it is annoying to constat that going to higher order may degrade the prediction, even if this concerns limiting cases only. To avoid this, it is possible to renormalized the reflection coefficient (37) identifying (1 + iaε)/(1− iaε) to e2iaε which is true up to O(ε2). We get Rren = eiψ, ψ ≡ 2k cos θ (ϕe+ C tan2 θ) , (38) and this renormalized version goes to the expected values for ϕ = 0, 1. We report in Fig. 8 two series of results. The color scale panels show |Rh − R| as a function of ϕ for a varying between −e and e. The error decreases when a increases, which is a bad news since a has an upper bound (a ≤ 0) in order to ensure a positive energy supported by the effective boundary (Eb > 0), as discussed in the previous section. Next, the error increases when ϕ or θ increases and in both cases, this is because R departs from the value R = 1 obtained in the absence of roughnesses, see e.g. (38). Next, the upper panels reports |Rren − R| (black lines) as a function of ϕ together with the set of errors |Rh − R| corresponding to a < 0. As expected, the error vanishes for ϕ→ 0; this is less true for ϕ = 1 except for θ which does not question C (when C is concerned, only the singular value (a) 11 . . . .. . .. . X2 X1 ⌃ U r 0 e 'h U inc ✓ . . . .. . .. . ⌃ Uh,r X2 X1a 0 U inc ✓ (a) (b) FIG. 6. Reflection of a wave U inc at oblique incidence on the rough boundary of the elastic body; (a) in the real problem with U r the reflected wave and (b) in the homogenized problem with Uh,r the reflected wave (the system of coordinate has been shifted of a, here a < 0; light grey region show the actual roughnesses, but they do not exist in the homogenized problem). In the above system, the origine of the X1-axis has been shifted of a (Fig. 6); to keep the same definition of the reflection coefficient, we want Rh = U ,r|Σ U inc|Σ , (35) and now Σ is the plane X1 = a. An exact solution of (23) can be found, accounting for (35), of the form Uh(X) = [ e−ik cos θ(X1−a) +Rheik cos θ(X1−a) ] eik sin θX2 , (36) and applying the boundary condition of (23) at X1 = 0, we have Rh = 1 + ikh cos θ za(θ) 1− ikh cos θ za(θ) e 2ika cos θ, za(θ) ≡ ( ϕ e h − a h ) + C tan2 θ. (37) To begin with, we report in Fig. 7 an example of the wavefield in the actual problem together with its homogenized counterparts obtained for a = −e and a = 0 (the configuration is described in the figure caption). The agreement |U − Uh|/|U | (in L2 norm) is of 20% for a = −e and of 5% for a = 0, which suggests that small a value is preferable. In fact, something is not very convenient in the form of (37). We have found that at leading order, the roughnesses behave as a flat free boundary, and as previously said, this will be exact in the two limiting cases ϕ = 0 (leading to R = 1) and ϕ = 1 (leading to R = e2ike cos θ). In both cases, the elementary problem (17) has a trivial solution, W = constant, since Γ is a boundary y1 = constant, whence C = 0. But neither e nor a being zero, Rh does not go to the expected values; strictly, this is possible since the homogenized solution is an approximation of the real one, but it is annoying to constat that going to higher order may degrade the prediction, even if this concerns limiting cases only. To avoid this, it is possible to renormalized the reflection coefficient (37) identifying (1 + iaε)/(1− iaε) to e2iaε which is true up to O(ε2). We get Rren = eiψ, ψ ≡ 2k cos θ (ϕe+ C tan2 θ) , (38) and this r normalized version goes to the expected values for ϕ = 0, 1. We report in Fig. 8 two series of results. The color scale panels show |Rh − R| as a function of ϕ for a varying between −e and e. The error decreases when a increases, which is a bad news since a has an upper bound (a ≤ 0) in order to ensure a positive energy supported by the effective boundary (Eb > 0), as discussed in the previous section. Next, the error increases when ϕ or θ increases and in both cases, this is because R departs from the value R = 1 obtained in the absence of roughnesses, see e.g. (38). Next, the upper panels reports |Rren − R| (black lines) as a function of ϕ together with the set of errors |Rh − R| corresponding to a < 0. As expected, the error vanishes for ϕ→ 0; this is less true for ϕ = 1 except for θ which does not question C (when C is concerned, only the singular value (b) Fig. 6. Reflection of a w ve U inc at oblique incidence on the rough boundary of the elastic body; (a) in the real problem with Ur the reflected wave and (b) in the homogenized problem with Uh,r the reflected wave (the system of coordinate has been shifted of a, here a < 0; light grey region show the actual roughnesses, but they do not exist in the homogenized proble ) In the above system, the origine of the X1-axis has been shifted of a (Fig. 6); to keep the same defin tion of the reflection coefficient, we want Rh = Uh,r|Σ Uhic|Σ , (35) and now Σ is the plane X1 = a. An exact solution of (23) can be found, accounting for (35), of the form Uh(X) = [ e−ik cos θ(X1−a) Rheik cos θ(X1 a) eik sin θX2 , (36) and applying the b undary condition of (23) at X1 = 0, we have Rh = 1+ ikh cos θza(θ) 1− ikh cos θza(θ) e 2ika cos θ , za(θ) ≡ ( ϕ e h − a h ) + C tan2 θ. (37) To begin with, we report in Fig. 7 an example of the wavefield in the actual problem together with its homogenized counterparts obtained for a = −e and a = 0 (the configu- ration is desc ibed in the figure caption). The ag eement ∣∣∣U −Uh∣∣∣ /|U| (in L2 norm) is of 20% for a = −e and f 5% for = 0, which suggests that mall a value is prefer ble. In fact, somethi g is no very convenient in the form of (37). We have found that at leading order, the roughnesses behave as a flat free bou dary, and as previously said, this will be exa t in t e two limiting cases ϕ = 0 (leading to R = 1) nd ϕ = 1 (leading to R = e2ike cos θ). In both cases, the lementary problem (17) has a trivial solution, W = consta t, since Γ is a boundary y1 = constant, whence C = 0. But neither e nor a being zero, Rh does not go to the expected values; strictly, this is possible since the homogenized solution s an approximation of the real one, but it is ann ying to constat that going to higher order 316 Agne`s Maurel, Jean-Jacques Marigo, Kim Pham may degrade the prediction, even if this concerns limiting cases only. To avoid this, it is possible to renormalized the reflection coefficient (37) identifying (1 + iaε)/(1− iaε) to e2iaε which is true up to O ( ε2 ) . We get Rren = eiψ, ψ ≡ 2k cos θ (ϕe + C tan2 θ) , (38) and this renormalized version goes to the expected values for ϕ = 0, 1. 12 a = 0 ⌃ (b) Uh(X), a = e (c) Uh(X), a = 0(a) U(X) FIG. 7. Reflection of an incident wave on the free edge of an elastic body with step roughnesses (kh = 1, e/h = ϕ = 0.5, θ = 20◦). (a) the solution in the actual problem U(X), (b-c) the homogenized solutions Uh(X), (36)-(37), for a = −e and a = 0. ϕ = 1 would lead to Rren = R for C = 0; but until ϕ = 1, C increases). From (a) and (b) , we observe that Rren is a better estimation of R and we built it indeed to do a better job; unexpectedly, it happens that this is not the case for high ϕ values near the grazing angles and we do not have explanation for that. All the observations reported in Fig. 8 are recovered for all k and e values; increasing k or e only produce a global increase of the errors. e/h e/h a/h 0 ' 0 1 ' 0 1 ' 0 1 3 4 5 6 106 104 102 Eb > 0 |R R r e n | . (a) θ = 0 (b) θ = 40◦ (c) θ = 80◦ FIG. 8. Top panels report the error |R−Rren| as a function of ϕ (black lines) and the set of the errors |R−Rh| for 0 < a < e (grey zones). Bottom panels report the error |R−Rh| (in colorscale) in the plane (ϕ, a/h); the dotted black line indicates the limit value of a/h below which Eb > 0. IV. CONCLUDING REMARKS We have studied a homogenized problem which can replace the actual problem of the reflection of shear waves at the rough free boundary of an elastic body. Parameters characteristic of an equivalent flat boundary enter in a boundary condition which differs from the usual stress free condition. We have inspected different forms of the boundary conditions, being all equivalent up to the order of validity of the model. This has been done addressing (a) U(X) 12 a = 0 ⌃ (b) Uh(X), a = e (c) Uh(X), a = 0(a) U(X) FIG. 7. R flection of an incident wave on the fre edge of an elastic body with step roughnesses (kh = 1, e/h = ϕ = 0.5, θ = 20◦). (a) the solution in the actual problem U(X), (b-c) the homogenized solutions Uh(X), (36)-(37), for a = −e and a = 0. ϕ = 1 would lead to Rren = R for C = 0; but until ϕ = 1, C increases). From (a) and (b) , we observe that Rren is a better estimation of R and we built it indeed to do a better job; unexpectedly, it ha pens that this is not the case for high ϕ values near the grazing angles and we do not have explanation for that. All the observations reported in Fig. 8 are recovered for all k and e values; increasing k or e only produce a global increase of the errors. e/h e/h a/h 0 ' 0 1 ' 0 1 ' 0 1 3 4 5 6 106 104 102 Eb > 0 |R R r e n | . (a) θ = 0 (b) θ = 40◦ (c) θ = 80◦ FIG. 8. Top panels report the error |R−Rren| as a function of ϕ (black lines) and the set of the errors |R−Rh| for 0 < a < e (grey zones). Bottom panels report the error |R−Rh| (in colorscale) in the plane (ϕ, a/h); the dotted black line indicates the limit value of a/h below which Eb > 0. IV. CONCLUDING REMARKS We have studied a homogenized problem which can replace the actual problem of the reflection of shear waves at the rough free boundary of an elastic body. Parameters characteristic of an equivalent flat boundary enter in a boundary condition which differs from the usual stress free condition. We have inspected different forms of the boundary conditions, being all equivalent up to the order of validity of the model. This has been done addressing (b) Uh(X), a = −e 12 a = 0 ⌃ (b) Uh(X), a = e (c) Uh(X), a = 0(a) U(X) FIG. 7. Reflection of an incident wave on th free edge of n elastic body with step roughnesses (kh = 1, e/h ϕ = 0.5, θ = 20◦). (a) the s lution in the ctual problem U(X), (b-c) the homogenized s lutions Uh(X), 6 -(37), for a = −e and a = 0. ϕ = 1 would lead to Rren = R for C = 0; bu until ϕ = 1, C increases). From (a) and (b) , we observe that Rren is a better es imation of R and we built it indee to do a better job; unexpectedly, it happens that this is not the case for high ϕ values near the grazi angles and we do not have explanation for that. All the observations eported in Fig. 8 are r covered for all k and e values; increasing k r e only produce a global increase of the errors. e/h e/h a/h 0 ' 0 1 ' 0 1 ' 0 1 3 4 5 6 106 104 102 Eb > 0 |R R r e n | . (a) θ = 0 (b) θ = 40◦ (c) θ = 80◦ FIG. 8. Top panels report the error | −Rren| as a fu ction of ϕ (black lines) and the set of the errors | −Rh| for 0 a < e (grey zones). Bottom panels report the error | −Rh| (in colorscale) in the plane (ϕ, a/ ); the dotted black line indicates the limit value of a/h below which Eb > 0. IV. CONCLUDING REMARKS We have studied a homogenized problem which can replace the ctual problem of the reflection of she r waves at the rough free boundary of an elastic body. Parameters characteristic of an equivalent flat boundary enter in a boundary condition which differs from the usual stress free condition. We have inspected different forms of the boundary conditions, being all equivalent up to the order of validity of the model. This has bee done addressing (c) Uh(X), a = 0 Fig. 7. Reflection of an incident wave on the free edge of an elastic body with step roughnesses (kh = 1, e/h = ϕ = 0.5, θ = 20◦). (a) the solution in the actual problem U (X), (b-c) the homoge- nized solutions Uh(X), (36)–(37), for a = −e and a = 0 We report in Fig. 8 two series of results. Th color scale p nels show ∣∣∣Rh − R∣∣∣ as fu ction of ϕ for a varyi g between − and e. The rror cre ses whe a increases, which is a bad ew si ce a has an upp r bou d (a ≤ 0) in order to ensu e a po itive energy supported by the ffective boundary (Eb > 0), as discussed in th pr vi us sec- tion. Next, the error increases when ϕ or θ increases and in both cases, this is because R departs from the value R = 1 obtained in the absence of roughnesses, see e.g. (38). Next, the upper panels reports |Rren − R| (black lines) as a function of ϕ together with the set of errors ∣∣∣Rh − R∣∣∣ corresponding to a < 0. As expected, the error vanishes for ϕ → 0; this is less true for ϕ = 1 except for θ which does not question C (when C is concerned, only the singular value ϕ = 1 would lead to Rren = R for C = 0; but until ϕ = 1, C increases). From (a) and (b) , we observe that Rren is a better estimation of R and we built it indeed to do a better job; unexpectedly, it happens that this is not the case for high ϕ values near the grazing angles and we do not have explanation for that. All the observations reported in Fig. 8 are recovered for all k and e values; increasing k or e only produce a global increase of the errors. Effective boundary condition for the reflection of shear waves at the periodic rough boundary of an elastic body 317 12 a = 0 ⌃ (b) Uh(X), a = e (c) Uh(X), a = 0(a) U(X) FIG. 7. Reflection of an incident wave on the free edge of an elastic body with step roughnesses (kh = 1, e/h = ϕ = 0.5, θ = 20◦). (a) the solution in the actual problem U(X), (b-c) the homogenized solutions Uh(X), (36)-(37), for a = −e and a = 0. ϕ = 1 would lead to Rren = R for C = 0; but until ϕ = 1, C increases). From (a) and (b) , we observe that Rren is a better estimation of R and we built it indeed to do a better job; unexpectedly, it happens that this is not the case for high ϕ values near the grazing angles and we do not have explanation for that. All the observations reported in Fig. 8 are recovered for all k and e values; increasing k or e only produce a global increase of the errors. e/h e/h a/h 0 ' 0 1 ' 0 1 ' 0 1 3 4 5 6 106 104 102 Eb > 0 |R R r e n | . (a) θ = 0 (b) θ = 40◦ (c) θ = 80◦ FIG. 8. Top panels report the error |R−Rren| as a function of ϕ (black lines) and the set of the errors |R−Rh| for 0 < a < e (grey zones). Bottom panels report the error |R−Rh| (in colorscale) in the plane (ϕ, a/h); the dotted black line indicates the limit value of a/h below which Eb > 0. IV. CONCLUDING REMARKS We have studied a homogenized problem which can replace the actual problem of the reflection of shear waves at the rough free boundary of an elastic body. Parameters characteristic of an equivalent flat boundary enter in a boundary condition which differs from the usual stress free condition. We have inspected different forms of the boundary conditions, being all equivalent up to the order of validity of the model. This has been done addressing (a) θ = 0 12 a = 0 ⌃ (b) Uh(X), a = e (c) Uh(X), a = 0(a) U(X) FIG. 7. Reflection of a incident wave on the free edge of an elastic body with step roughnesses (kh = 1, e/h = ϕ = 0.5, θ = 20◦). (a) the solution in the actual problem U(X), (b-c) th homogenized solutions Uh(X), (36)-(37), for = −e and a = 0. ϕ = 1 would lead to ren = R for C = 0; but until ϕ = 1, C increases). From a) and (b) , w observe that Rren is a better estimatio of R and we built it in eed to do a better job; unexpectedly, it happens that this is not the case for high ϕ v lues near the grazing angles and we do not h ve explanation for that. All the observations reported in Fig. 8 are rec vered for all k and e values; increasing k or e only produce a global increase of the errors. e/h e/h a/h 0 ' 0 1 ' 0 1 ' 0 1 3 4 5 6 106 104 102 Eb > 0 |R R r e n | . (a) θ = 0 (b) θ = 40◦ (c) θ = 80◦ FIG. 8. Top panels r po t the error |R−Rren| as a function of ϕ (black lines) and the set of the errors |R−Rh| for 0 < a < e (grey z nes). Bottom panels r po t the error |R−Rh| (in colorscale) in the plane (ϕ, a/h); the dotted black lin indicates the limit value of a/h below which Eb > 0. IV. CONCLUDING EMARKS We have studied a homogenized problem which can replace the actual problem of the reflection of sh ar waves at t e ough free boundary of an elastic body. Parameters characteristic of an equivalent flat bou dary enter in a boundary condition which differs from the usual stress free condition. We have inspected different forms of the boundary conditions, being all equivalent up to the order of validity of the model. This has b en one addressing (b) θ = 40◦ 12 a = 0 ⌃ (b) Uh(X), a = e (c) Uh(X), a = 0(a) U(X) FIG. 7. Reflectio of an incident wav n the free edge of an elastic body with step roughnesses (kh = 1, e/h = ϕ = 0.5, θ = 20◦). (a) the solution in the actual problem U(X), (b-c) th homogenized solutions Uh(X), (36)-(37), for = −e and a = 0. ϕ = 1 would lead to Rren = R for C = 0; but until ϕ = 1, C increases). From (a) and (b) , we observe that Rren is a better estimatio of R and we built it in ed to do a b tt r job; unex ctedly, it happens at this is not the case for high ϕ v lues near the grazing angles and we do not h ve explanation for that. All th obs rvations reported in Fig. 8 are rec vered for all k and e values; increasing k or e only produce a global increase of the errors. e/h e/h a/h 0 ' 0 1 ' 0 1 ' 0 1 3 4 5 6 106 104 102 Eb > 0 |R R r e n | . (a) θ = 0 (b) θ = 40◦ (c) θ = 80◦ FIG. 8. Top anels report the error |R−Rren| as a function of ϕ (black lines) and the set of the errors |R−Rh| for 0 < a < e (grey z ne ). B t om anels repo t the err r |R−Rh| (in co orscale) in the plane (ϕ, a/h); the dotted black line indicates the limit value of a/h below which Eb > 0. IV. CONCLUDING EMARKS We ave studi a homogenized probl m whic c n replace the actual problem of the reflection of sh ar waves at t e rough free boundary of an elastic body. Parameters characteristic of n eq ivalent flat bou dary enter in b undary condition whic differs from the usual s ress fre condition. We hav i spected differ nt orms of the b undary condit ons, b ing al equivalen up to the order of validity of the mod l. This has b e one addressing (c) θ = 80◦ Fig. 8. Top panels report t e error |R− Rren| as a function of ϕ (black lines) and the set of the errors ∣∣∣R− Rh∣∣∣ for 0 < a < e (grey zones). Bottom panels report the error |R− |Rh| (in colorscale) in the plane (ϕ, a/h); the dotted black line indicates the limit value of a/h below which Eb > 0 4. CONCLUDING REMARKS We hav stu ied homogeniz pr blem wh ch can replace the actual problem of the reflection of shear waves at the rough free boundary of an elastic body. Parameters characteristic of an equivalent flat boundary enter in a boundary condition which differs from the usual stress free condition. We have inspected different forms of the boundary conditions, being all equivalent up to the order of validity of the model. This has been done addressing two aspects, (i) wether or not the homogenized problem is well suited for a numerical resolution in the time domain (which means free of numerical instabil- ities in the time computation) and (ii) which formulation gives the smallest error in the model when compared to the solution of the actual problem. The first aspect questions the equation of energy conservation; the homogenized boundary condition makes an additional term of energy to appear and this energy must be positive in order ensure a consistent computational method in the time domain (and as so, it is not optional). The second aspect is of less importance; firstly because no definitive answer can be given, the error being measured in a particular scattering problem with no guaranty that the optimal homogenized problem will be the same in another scattering problem. REFERENCES [1] D. Cioranescu and P. Donato. An introduction to homogenization, Vol. 4. The Clarendon Press Oxford University Press, New York, (1999). 318 Agne`s Maurel, Jean-Jacques Marigo, Kim Pham [2] J. J. Marigo and C. Pideri. The effective behavior of elastic bodies containing microcracks or microholes localized on a surface. International Journal of Damage Mechanics, 20, (8), (2011), pp. 1151–1177. https://doi.org/10.1177/1056789511406914. [3] M. David, J. J. Marigo, and C. Pideri. Homogenized interface model describing in- homogeneities located on a surface. Journal of Elasticity, 109, (2), (2012), pp. 153–187. https://doi.org/10.1007/s10659-012-9374-5. [4] A. Bonnet-Bendhia, D. Drissi, and N. Gmati. Simulation of muffler’s transmission losses by a homogenized finite element method. Journal of Computational Acoustics, 12, (03), (2004), pp. 447–474. https://doi.org/10.1142/s0218396x04002304. [5] J. J. Marigo and A. Maurel. Homogenization models for thin rigid structured surfaces and films. The Journal of the Acoustical Society of America, 140, (1), (2016), pp. 260–273. https://doi.org/10.1121/1.4954756. [6] B. Delourme, H. Haddar, and P. Joly. Approximate models for wave propagation across thin periodic interfaces. Journal de Mathe´matiques Pures et Applique´es, 98, (1), (2012), pp. 28–71. https://doi.org/10.1016/j.matpur.2012.01.003. [7] B. Delourme. High-order asymptotics for the electromagnetic scattering by thin peri- odic layers. Mathematical Methods in the Applied Sciences, 38, (5), (2015), pp. 811–833. https://doi.org/10.1002/mma.3110. [8] A. Maurel, J. J. Marigo, and A. Ourir. Homogenization of ultrathin metallo-dielectric struc- tures leading to transmission conditions at an equivalent interface. The Journal of the Optical Society of America B, 33, (5), (2016), pp. 947–956. https://doi.org/10.1364/josab.33.000947. [9] J. J. Marigo and A. Maurel. Two-scale homogenization to determine effective parameters of thin metallic-structured films. Proc. R. Soc. A, 472, (2016). [10] L. Rayleigh. On the dynamical theory of gratings. In Proceedings of the Royal Society of London, Series A, Vol. 79, (1907), pp. 399–416. [11] P. A. Martin and R. A. Dalrymple. Scattering of long waves by cylindrical obstacles and grat- ings using matched asymptotic expansions. Journal of Fluid Mechanics, 188, (1988), pp. 465– 490. https://doi.org/10.1017/s0022112088000801. APPENDIX 1. SIMPLE PROCEDURE TO GET C FOR RECTANGULAR INCLUSIONS Mode matching is a simple way to get C. We consider the solution W˜ = W + y2 satisfying ∆W˜ = 0,∇W˜.n|∂V = 0 and W˜ → y2 for y1 → ∞. The field W˜ can be written W˜(y) =  W−(y) = N− ∑ n=1 w−n cosh an (y1 + e) cosh ane W−n (y2) , 0 ≥ y1 ≥ −e W+(y) = y2 + N+ ∑ n=−N+,n 6=0 w+n e −|bn|y1W+n (y2) , y1 ≥ 0 (39) with an = npi/ϕ, bn = 2npi, and W−n (y2) = √ 2 ϕ cos ( any2 + npi 2 ) , W+n (y2) = e ibny2 , (40) the transverse functions (forming a basis) adapted for solutions being respectively peri- odic and with zero derivatives at y2 = ±ϕ/2 (Fig. 9). Effective boundary condition for the reflection of shear waves at the periodic rough boundary of an elastic body 319 13 two aspects, (i) wether or not the homogenized problem is well suited for a numerical resolution in the time domain (which means free of numerical instabilities in the time computation) and (ii) which formulation gives the smallest error in the model when compared to the solution of the actual problem. The first aspect questions the equation of energy conservation; the homogenized boundary condition makes an additional term of energy to appear and this energy must be positive in order ensure a consistent computational method in the time domain (and as so, it is not optional). The second aspect is of less importance; firstly because no definitive answer can be given, the error being measured in a particular scattering problem with no guaranty that the optimal homogenized problem will be the same in another scattering problem. [1] D. Cioranescu and P. Donato, The Clarendon Press Oxford University Press, New York 4, 118 (1999). [2] J.-J. Marigo and C. Pideri, International Journal of Damage Mechanics , 1056789511406914 (2011). [3] M. David, J.-J. Marigo, and C. Pideri, Journal of Elasticity 109, 153 (2012). [4] A. Bonnet-Bendhia, D. Drissi, and N. Gmati, Journal of computational acoustics 12, 447 (2004). [5] J.-J. Marigo and A. Maurel, The Journal of the Acoustical Society of America 140, 260 (2016). [6] B. Delourme, H. Haddar, and P. Joly, Journal de mathe´matiques pures et applique´es 98, 28 (2012). [7] B. Delourme, Mathematical Methods in the Applied Sciences 38, 811 (2015). [8] A. Maurel, J.-J. Marigo, and A. Ourir, The Journal of the Optical Society of America B 33, 947 (2016). [9] J.-J. Marigo and A. Maurel, in Proc. R. Soc. A, Vol. 472 (The Royal Society, 2016) p. 20160068. [10] L. Rayleigh, Proceedings of the Royal Society of London. Series A, 79, 399 (1907). [11] P. A. Martin and R. A. Dalrymple, Journal of Fluid Mechanics 188, 465 (1988). 1. Simple procedure to get C for rectangular inclusions Mode matching is a simple way to get C. We consider the solution W˜ = W + y2 satisfying ∆W˜ = 0, ∇W˜ .n|∂V = 0 and W˜ → y2 for y1 →∞. The field W˜ can be written W˜ (y) =  W−(y) = N−∑ n=1 w−n cosh an(y1 + e) cosh ane W−n (y2), 0 ≥ y1 ≥ −e W+(y) = y2 + N+∑ n=−N+,n6=0 w+n e −|bn|y1 W+n (y2), y1 ≥ 0, (39) with an = npi/ϕ, bn = 2npi, and W−n (y2) = √ 2 ϕ cos ( any2 + npi 2 ) , W+n (y2) = e ibny2 , (40) the transverse functions (forming a basis) adapted for solutions being respectively periodic and with zero derivatives at y2 = ±ϕ/2. y1 y2 0e/h '/2 '/2 1/2 1/2 W(y) W+(y) FIG. 9. Mode matching configuration. The solution W± is written for y1 > −e/h, (39), and the resolution involves only matching conditions at y1 = 0. Now, we will ask to W± to match (on average) their value and their first derivative at y1 = 0, and this latter matching on the derivative will be done accounting for the boundary conditions at y1 = 0 and |y2| > ϕ/2 (note that Fig. 9. Mode matching configuration. The solution W± is written for y1 > −e/h, (39), and the resolution involves only matching conditions at y1 = 0 Now, we will ask to W± to match (on average) their value and their first derivative at y1 = 0, and this latter matching on the derivative will be done accounting for the boundary conditions at y1 = 0 and |y2| > ϕ/2 (note that W− satisfies by construction the right boundary condition on Γ, both at the wall y1 = 0 and on the boundary of the roughness at y2 = ±ϕ/2). To that aim, we use the following relations ∫ ϕ/2 −ϕ/2 dy2W− (0, y2)W−m (y2) = ∫ ϕ/2 −ϕ/2 dy2W+ (0, y2)W−m (y2),∫ ϕ/2 −ϕ/2 dy2 ∂W− ∂y1 (0, y2)W+∗m (y2) = ∫ 1/2 −1/2 dy2 ∂W+ ∂y1 (0, y2)W+∗m (y2) , (41) with W+∗m the conjugate of W+m (W−m is real). The first relation is the matching of the values in the region where W−m is defined. The second relation has more information: we have used that the ∂y1W + = 0 for |y2| > ϕ/2, from which∫ 1/2 −1/2 dy2 ∂W+ ∂y1 (0, y2)W+∗n (y2) = ∫ ϕ/2 −ϕ/2 dy2 ∂W+ ∂y1 (0, y2)W+ ∗ m (y2) , (42) afterwards we ask, on average, ∂y1W + = ∂y1W − for |y2| < ϕ/2. We get for a matrix system for the two vectors w− = ( w−n ) n=0,...,N− and w + = ( w+n ) n=0,...,N+( I −tF∗ FA tanh(Ae) B )( w− w+ ) = ( S 0 ) , (43) with I the N− × N− identity matrix, A = diag (an) ,B = diag (|bn|), and Fmn =∫ ϕ/2 −ϕ/2 dy2W+∗m (y2)W−n (y2) and Sn = ∫ ϕ/2 −ϕ/2 dy2y2W−n (y2). The expressions of Fmn and Sn are given below Fmn = √ ϕ 2 [ sinc ((an − bm) ϕ/2) einpi/2 + sinc ((an + bm) ϕ/2) e−inpi/2 ] , Sn = −2 √ 2 ϕ 1 a2n . (44) 320 Agne`s Maurel, Jean-Jacques Marigo, Kim Pham The system of the form Mw = s with the matrix M being square (this is not always the case in systems written using mode matching). Next, M is invertible and the system can be solved to find w in the least squares sense (as done by the operation M\s in Matlab). Then, we want to determine C = − ∫ dy ∂W ∂y2 = ∫ 0 −e dy1 ∫ ϕ/2 −ϕ/2 [ 1− ∂W − ∂y2 ] . (45) where we have used that W (y1 ≥ 0, y2) = W+(y) − y2 is periodic, thus of vanishing contribution. It is now sufficient to write C = eϕ− w−n tanh ane/an [ W−n ]ϕ/2 −ϕ/2 to get C = eϕ+ 2 √ 2 ϕ∑n tanh ane an w−n . (46) The procedure of mode matching is longer to explain than to encode; below is a script working with Matlab. func t ion C=CoefC ( phi , e ,Nd,Np) nd = 1 : 2 :Nd; np=[−Np:−1 ,1 :Np ] ; Nd=length ( nd ) ; Np=length ( np ) ; an=nd∗pi/phi ; bn=2∗np∗pi ; f o r mm=1:Np, f o r nn=1:Nd, a=an ( nn ) ; b=bn (mm) ; n=nd ( nn ) ; EX=exp (1 i ∗n∗pi / 2 ) ; temp= s i n c ( ( a−b )∗ phi /(2∗ pi ) ) ∗EX+ s i n c ( ( a+b )∗ phi /(2∗ pi ) ) /EX ; F (mm, nn)= s q r t ( phi /2)∗ temp ; end end M=[ eye (Nd) ,−F ’ ; F∗diag ( an .∗ tanh ( an∗e ) ) , diag ( abs ( bn ) ) ] ; s= s q r t (2/ phi ) ∗ ( ( ( − 1 ) . ˆ nd−1)./an . ˆ 2 ) . ’ ; S=[ s ; zeros (Np, 1 ) ] ; w=M\S ;wm=q ( 1 :Nd) ; C=phi∗e+2∗ s q r t (2/ phi )∗sum(wm. ’ . ∗ tanh ( an∗e ) . / an ) ; end APPENDIX 2. THE SAME PROCEDURE TO SOLVE THE SCATTERING PROBLEM The procedure to solve the scattering problem for an incident plane wave is quasi the same than to solve the static elementary problem. Now, we look for U(X) being the Effective boundary condition for the reflection of shear waves at the periodic rough boundary of an elastic body 321 solution of (32) (and the whole story is to compute Uev(X) in X1 > 0, which requires to compute U (X1 < 0, X2) in the roughnesses). U(X) =  U−(X) = N− ∑ n=1 un cos an (X1 + e) cos ane U−n (X2) , 0 ≥ X1 ≥ −e U+(X) = e−iknX1U+0 (X2) + N+ ∑ n=−N+ RneiknX1U+n (y2) , X1 ≥ 0 (47) and with αn ≡ npi/ϕ as previously, the bases are defined as U−n (X2) = √ 2− δn0 ϕ cos ( αnX2 + npi 2 ) , a2n = k 2 − α2n, which ensures that U− satisfies the right Neumann boundary conditions for X1 < 0, |X2| = ϕ/2 and satisfies the Helmholtz equation ∆U− + k2U− = 0, and U+n (X2) = e iβnX2 , β2n ≡ k2 − ( k sin θ + 2npi h )2 . Again, U+n is chosen in order to ensure (i) that U + satisfies the Helmholtz equation ∆U+ + k2U+ = 0 and (ii) that U+ satisfies the so-called condition of pseudo-periodicity U+ (X1, h) = U+ (X1, 0) eiβ0h (this condition is imposed by the form of the incident wave). Now the same matching as (41) are used, leading to a matrix condition involving Gmn = ∫ ϕ/2 −ϕ/2 dX2U+∗m (X2)U−n (X2), Gmn =  √ ϕ sinc (βmϕ/2) , n = 0√ ϕ 2 [ sinc ((αn − βm) ϕ/2) einpi/2 + sinc ((αn + βm) ϕ/2) e−inpi/2 ] , n 6= 0 (48) We get a system of the form( −tG∗ I iK GA tan(Ae) )( R u ) = ( s1 s2 ) , (49) The source terms are now due to the incident wave and we get S1,n = G∗0n, S2,n = ik0δn0. (50) The scripts below allow to compute the reflection coefficients Rn (and R = R0). For a graphical representation of the solution, see below. funct ion [R , Rn , un]= Sca t te r ingPb ( k , theta , e , phi ,Np,Nn) np=−Np:Np; no=Np+1; nd=0:Nn; Np=length ( np ) ; Nd=length ( nd ) ; beta0=k∗ s in ( t h e t a ) ; betan=beta0 +2∗np∗pi ; 322 Agne`s Maurel, Jean-Jacques Marigo, Kim Pham k0= s q r t ( kˆ2−beta0 ˆ 2 ) ; kn= s q r t ( kˆ2−betan . ˆ 2 ) ; alphan=nd∗pi/phi ; an= s q r t ( kˆ2−alphan . ˆ 2 ) ; f o r mm=1:Np, b=betan (mm) ; G(mm, 1 ) = s q r t ( phi )∗ s i n c ( b∗phi /(2∗ pi ) ) ; f o r nn=2:Nd, a=alphan ( nn ) ; Ex=exp (1 i ∗nd ( nn )∗ pi / 2 ) ; temp= s i n c ( ( a−b )∗ phi /(2∗ pi ) ) ∗Ex+ s i n c ( ( a+b )∗ phi /(2∗ pi ) ) / Ex ; G(mm, nn)= s q r t ( phi /2)∗ temp ; end end M=[−G’ eye (Nd) ; diag (1 i ∗kn ) G∗diag ( an .∗ tan ( an∗e ) ) ] ; S1=G( no , : ) ’ ; S2=0∗kn ; S2 ( no)=1 i ∗k0 ; S=[ S1 ; S2 . ’ ] ; V=M\S ; Rn=V( 1 :Np ) ; un=V(Np+1: end ) ; R=Rn( no ) ; The graphical representation of the solution is done simply following (47). However, is N− is too large, this may conduce to numerical divergences, because an become com- plex, thus cos an (X1 + e) / cos ane diverge for X1 < 0. Obviously this limitation does not concerns the resolution of the system as presented above (since the mode matching is performed at X1 = 0 and does not question X1 ¡ 0). As given below, we simply use a trick to avoid divergence of the solution for X1 < 0) (truncating the solution for diverging cos ane. Alternatively, one has to consider the solution written as U−(X) = N− ∑ n=1 [ Aneian(X1+e) + Bne−ianX1 ] U−n (X2) , 0 ≥ X1 ≥ −e, (51) and to apply matching conditions at X1 = −e and X1 = 0, and at X1 = 0, it is∫ hϕ/2 −hϕ/2 dX2 ∂U− ∂X1 (−e, X2)U−m (X2) = 0, to get a final system on (An, Bn, Rn). func t ion Graphic ( k , theta , e , phi , Rn , un , L , Nrough , dx ) %truncat ioni fneeded , alphan = ( 0 : length ( un ) ) ∗ pi/phi ; an= s q r t ( kˆ2−alphan . ˆ 2 ) ; [max ,Nd]=min ( abs ( cos ( an∗e)−1e10 ) ) ; %% Np=( length (Rn)−1)/2; np=−Np:Np; nd=0:Nd−1; Np=length ( np ) ; Nd=length ( nd ) ; Effective boundary condition for the reflection of shear waves at the periodic rough boundary of an elastic body 323 beta0=k∗ s in ( t h e t a ) ; betan=beta0 +2∗np∗pi ; k0= s q r t ( kˆ2−beta0 ˆ 2 ) ; kn= s q r t ( kˆ2−betan . ˆ 2 ) ; N1= f l o o r ( e/dx ) + 1 0 ;N2= f l o o r ( phi/dx ) + 1 0 ; x1m= l i n s p a c e (−e , 0 ,N1) ’∗ ones ( 1 ,N2 ) ; y2m=ones (N1, 1 ) ∗ l i n s p a c e (−phi /2 , phi /2 ,N2 ) ; Unn=1/ s q r t ( phi ) ;Um=un ( 1 ) / cos ( an ( 1 )∗ e )∗ cos ( an ( 1 ) ∗ ( x1m+e ) ) . ∗Unn ; f o r i i =2 : length ( nd ) , Unn=1/ s q r t ( phi /2)∗ cos ( alphan ( i i )∗y2m+nd ( i i )∗ pi / 2 ) ; Um=Um+ un ( i i )/ cos ( an ( i i )∗ e )∗ cos ( an ( i i ) ∗ ( x1m+e ) ) . ∗Unn ; end N1= f l o o r ( L/dx ) + 1 0 ;N2= f l o o r ( . 5 / dx ) + 1 0 ; x1p= l i n s p a c e ( 0 , L , N1) ’∗ ones ( 1 ,N2 ) ; x2p=ones (N1, 1 ) ∗ l i n s p a c e (−1/2 ,1/2 ,N2 ) ; Upn=exp (1 i ∗beta0 ∗x2p ) ; Up=exp(−1 i ∗k0∗x1p ) . ∗Upn ; f o r i i =1 : length ( np ) , Upn=exp (1 i ∗betan ( i i )∗ x2p ) ; Up=Up+Rn( i i )∗ exp (1 i ∗kn ( i i )∗ x1p ) . ∗Upn ; end f i g u r e f o r i i =1 :Nrough , phase=exp (1 i ∗ ( i i −1)∗beta0 ) ; pcolor (x1m , y2m+( i i −1) , r e a l (Um. ∗ phase ) ) , hold on , pcolor ( x1p , x2p +( i i −1) , r e a l (Up. ∗ phase ) ) , shading f l a t , a x i s auto end a x i s equal