A wave finite element approach for modelling wave transmission through laminated plate junctions

We present a numerical method for computing reflection and transmission coefficients at joints connecting composite laminated plates. The method is based on modelling joints with finite elements with boundary conditions given by the solutions of the wave finite element method for the plates in the infinite half-spaces connected to the joint. There are no restrictions on the number of plates, inter-plate angles, and material parameters of individual layers forming the composite. An L-shaped laminated plate junction is discussed in more detail. Comparisons of numerically predicted scattering coefficients with semi-analytical solutions for the selected structures are presented. The results obtained are essential for statistical energy analysis and dynamical energy analysis based calculations of the wave energy distribution in full built-up structure.

A wave finite element approach for modelling wave transmission through laminated plate junctions Nurkanat Aimakov 1* , Gregor Tanner 1 & Dimitrios Chronopoulos 2 We present a numerical method for computing reflection and transmission coefficients at joints connecting composite laminated plates. The method is based on modelling joints with finite elements with boundary conditions given by the solutions of the wave finite element method for the plates in the infinite half-spaces connected to the joint. There are no restrictions on the number of plates, inter-plate angles, and material parameters of individual layers forming the composite. An L-shaped laminated plate junction is discussed in more detail. Comparisons of numerically predicted scattering coefficients with semi-analytical solutions for the selected structures are presented. The results obtained are essential for statistical energy analysis and dynamical energy analysis based calculations of the wave energy distribution in full built-up structure.
Composites are widely used within the transport sector, in particular in the aerospace, automotive and naval manufacturing industries 1,2 . In comparison to isotropic materials such as aluminium and stainless steel, composites provide similar stiffness and strength characteristics whilst being significantly lighter 3 . Furthermore, the mechanical properties of fibre-reinforced composites can be tailored to suit particular needs 1,4 . Over the past decades, these advantages of composites have led to a growing number of use-cases for composites in the construction of primary structural components in the aerospace and automotive industries.
However, despite their superior structural characteristics, composites exhibit reduced vibro-acoustic performance levels due to the large variety of propagating wave modes. Thus, modelling noise and vibration in composite structures plays an important role both at the design phase of a vehicle and at the post-built stage, when non-destructive testing techniques help monitor the structure's performance. Therefore, there is a need for numerical methods to evaluate the vibrational response of composite structures fast and accurately.
Vibrations of a complex structure are in general modelled using deterministic schemes such as finite element (FE) 5 , finite difference (FD) 6 or boundary element (BE) methods 7 . These methods are particularly useful in providing the full phase and amplitude information of the wave field in the low-frequency regime. However, at higher frequencies, these methods become inefficient and computationally expensive as the model-size increases drastically with the frequency. Moreover, mode shapes and eigenfrequencies which are essential in the modal approach become highly sensitive to geometrical and/or material uncertainties of meshes, thus producing inaccurate results [8][9][10] .
At higher frequencies, numerical approaches such as SEA 11,12 , the radiative transfer method [13][14][15] or DEA [16][17][18] are favoured. For all these methods, wave propagation characteristics such as dispersion relations and scattering coefficients at discontinuities in the structure are required and routinely used. For complex materials such as composites, the dispersion curves and associated mode types can be obtained numerically using the wave finite element (WFE) method. Reflection and transmission behaviour at joints can be solved using FE tools as has been done for isotropic materials, see next paragraphs, and will be presented for composites in "Wave finite element method for composite plates". Alternatively, these scattering coefficients can be estimated using semi-analytical methods based on force-balance equations at the interface, see 19 for isotropic materials and 20 for composite plates.
The WFE method is a technique to study wave motion in homogeneous or periodic structures. The vibroacoustic behaviour of the whole structure can then be described through the analysis of a single FE-cell or using one periodic segment of the structure, respectively 21 . Since only one period of the structure is used, the size of the WFE model does not depend on the dimensions of the waveguide, and the computational cost of the method is low. In addition, conventional FE matrices are used to discretise the periodic cell; thus, the full potential of existing conventional FE tools can be exploited. The WFE method was originally proposed by Mead 22 describing the harmonic wave propagation in one-dimensional periodic systems. An important contribution to the analysis of wave propagation in various periodic structures using FE models has been made by Abdel-Rahman 23 . The free wave propagation in one-dimensional isotropic and anisotropic (composite) waveguides was analysed by Mace et al. 24 and Duhamel et al. 25 . Dispersion relations of waves were computed using the one-dimensional WFE method for structures such as stiffened cylinders 26 , car tyres 27 , thin-walled structures 28 , inhomogeneous cylindrical 29 and fluid-filled pipes 30 , sandwich beams and panels 31,32 and laminated cylinders 33 . A detailed analysis of the numerical implementation and numerous applications of the one-dimensional WFE method can be found in the work by Waki 34,35 . It is worth mentioning the work of Mencik and Ichchou 36 suggesting a special substructuring technique to address numerical issues in the case of multi-layered structures.
The basics of the WFE method for two-dimensional periodic systems were presented in the work by Mead 22 . Later, Mead and Parthan 37 showed how the problem of defining the dispersion relations in the general direction over the plate's plane dimensions could be reduced to an array of one-dimensional WFE problems with varying lengths of the periodic segments. A rigorous mathematical framework for the WFE method for two-dimensional periodic isotropic and composite systems has been developed by Manconi and Mace 38,39 . Several representations of the eigenvalue problem leading to the computation of dispersion relations were postulated. Alimonti et al. 40 extended this work by presenting a contour integral method to compute the non-linear eigenvalue problem arising from the governing equations of motion upon fixing the frequency and the direction of propagation. Dispersion relations were computed for two-dimensional arbitrarily thick layered panels in 41,42 and periodic textile composites in 43 .
Mencik and Ichchou introduced the hybrid FE/WFE method for calculating reflection and transmission coefficients for one-dimensional waveguides coupled longitudinally 44 . In recent years, this method has been extended to other types of junctions [45][46][47] and to two-dimensional waveguides 48,49 . However, the structures considered in the references listed above were all isotropic.
Beyond the case of wave propagation in isotropic materials, Chronopoulos 50 computed scattering coefficients at a junction representing damage between two composite beams. Later, Apolowo and Chronopoulos 51 computed the scattering coefficients of two multi-layer composite plates coupled longitudinally to localise the structural damage in the context of structural health monitoring. An attempt to extend the work of Renno et al. 45 to composite plates has been made by Mitrou and Renno in 52 . The results were not reliable, however, as the energy scattering coefficients did not sum to unity as expected in lossless systems 45,49 .
Beyond the WFE method, Karunasena and Shah 53 studied reflection of guided waves in the region of a bonding material connecting two composite plates using the hybrid FE and semi-analytical FE method. Bosmans et al. 54 studied the scattering properties of orthotropic plate junctions with principal material axes aligned with the plate coordinates, that is, so-called specially orthotropic plates. Results were presented only for the particular case of bending wave transmission loss in right-angled plates, so-called L-junctions. Aimakov et al. 20 developed a semi-analytical approach for calculating scattering matrices of junctions of orthotropic plates with no restrictions on the angles of orientation and of the principal material axes. Lee et al. 55 presented the scattering coefficients of coupled composite plates with joint compliance and damping using the First-Order laminated plate theory 56 . However, as in the work of Bosmans et al., the principal material axes of laminates considered in this work are aligned with the plate coordinates, effectively reducing the complexity of the underlying governing equations. Furthermore, in 55 , a shear correction factor is introduced to correct transverse shear stiffness in the laminate, which must be defined for each laminate separately.
In this paper, we extend the hybrid FE/WFE method to composite laminated plates. The principal novelty of this work is a detailed derivation of reflection and transmission matrices for waves travelling in the structural junctions connecting composite laminated plates at arbitrary angles and with an arbitrary material orientation of the principal axis of the laminae. We give, in particular, the scattering coefficients containing the full angleof-incidence and frequency dependence. For ray-based methods such as the radiative transfer method and DEA, detailed information on the reflection/transmission behaviour of all propagating modes at complex junctions is needed. This includes information about the angle-of-incidence dependence of scattering and mode conversion coefficients.
The manuscript is organised as follows: in "Wave finite element method for composite plates", the WFE method for modelling composite plates is reviewed. An eigenvalue problem whose solutions yield wave numbers and mode shapes is set up. The classification of the wave numbers and the wave basis setting are described. Having established a wave basis representation of displacement and force vectors in individual plates, we combine these solutions with the equations of motion of the joint fulfilling continuity of displacements and force equilibrium at the joint boundaries. This then yields the desired scattering coefficients as described in "Hybrid FE/WFE method and scattering coefficients". In "Numerical case examples", we present numerical case studies for two coupled composite plates. In particular, the energy scattering coefficients for L-type junctions of regular cross-ply and angle-ply composite plates are computed. These results are compared with semi-analytical estimates of energy scattering coefficients based on the work of Aimakov et al. 20 . Finally, concluding remarks are put in "Conclusion".

Wave finite element method for composite plates
Governing equations of motion. Consider a unit cell of a periodic or homogeneous composite plate with arbitrary lay-up through the thickness direction and plane dimensions d x and d y . (Note that for a homogeneous structure, the length scales of the unit cell are somewhat arbitrary and typically represented by a single finite element in the plate directions.) It can be modelled using three-dimensional solid elements (such as SOLID185 in the FE software ANSYS) stacked up one on top of the other, representing different composite layers. Figure 1 represents a unit cell of a five-layer plate meshed with SOLID185 elements and a nodal displacements vector labelling convention. The nodal displacements vector q is organised as The nodal forces vector f is arranged in the same manner. The number of degrees of freedom must be the same for each pair of edges on opposite faces. The number of mesh cells in the x, y and z direction are labelled by n x , n y and n z . The number of degrees of freedom per edge is labelled as m; for plates modelled with SOLID185 elements m = 3(n z + 1) . Consequently, the sizes of nodal displacement sub-vectors can be represented as Assuming that the structure undergoes harmonic vibration with angular frequency ω and no external forces are applied, we can write the governing equation of motion of the unit cell as where M and K are the mass and stiffness matrices, respectively. The parameter η denotes a uniform structural damping coefficient. The dimension of Eq. (2) is m(n x + 1)(n y + 1).
We consider a plane wave travelling across the plate to be of the form exp(−ik x x − ik y y + iωt) , where k x and k y are the x and y components of the wave vector k . Periodic structure theory 57,58 requires that where x = e −ik x d x is the propagation factor in the x direction. Now, we denote as q j X , X = L, R, I displacement sub-vectors of nodes placed at x j = d x j/n x , j = 1, 2, . . . , n x − 1 . Consequently, as in Eq. (3), one can relate internal and edge degrees of freedom to bottom ones as Following Eqs. (3) and (4), we can express the nodal displacements vector q in terms of displacement sub-vectors of nodes � . Figure 1. The set up of N plates joined together and a periodic segment of a plate with two different alternating plies modelled with three-dimensional finite elements. A ± 1,N are the amplitudes of incoming and outgoing waves travelling from infinity towards the junction and from the junction to infinity, respectively. The degrees of freedom are grouped into internal q I , edge q L , q R , q B , q T and corner q LB , q RB , q LT , q RT degrees of freedom. where Finally, by applying the periodic structure theory and force equilibrium in the y direction, which can be written as one can get from Eq. (8) the following eigenvalue problem for the propagation factor y

Scientific Reports
The dimension of the matrices D and Z is 2m. Therefore, the solution of the eigenvalue problem (11) consists of 2m propagation factors y,i and the correspondent eigenvectors φ q,i φ f ,i T provided that the circular frequency ω and the wave number component k x are fixed. The wave number components k y,i can be computed as Consequently, by varying the wave number component k x , one can extract wave vector curves (k x , k y ) for a fixed value of ω . It is worth noting that the so obtained wave number components k y,i can be real, imaginary or complex, making the correspondent plane waves propagating, evanescent or attenuating. www.nature.com/scientificreports/ Computation of wave numbers and eigenvectors. As mentioned in 35 , the solution of the eigenvalue problem (11) might be prone to ill-conditioning of the eigenvalue matrix Z and eigenvectors [q j L ,f j L ] T . An equivalent form of the eigenvalue problem can be formulated by eliminating f L and f R in Eq. (8) as suggested in 59 : with I being the m-by-m identity matrix. The matrices N and L consist of block parts of the reduced dynamic stiffness matrix D with no matrix inversion as in Eq. (11), thus reducing numerical ill-conditioning of the method to some extent. However, there might be a difference of several orders of magnitude between I and −D RL(RR) and D LL(LR) ; therefore, the condition numbers of the matrices N and L can still be large, causing numerical errors in the evaluation of eigenvalues and eigenvectors.
We use the following form of the eigenvalue problem (13), proposed by Fan et al. 47,60 where σ = �D RR � 2 m 2 , 2 representing the largest singular value of a matrix. This formulation is an improved version of the eigenvalue problem (13), and the factor σ is introduced here to reduce the condition number of the matrices N and L . All the formulations presented are equivalent, and the eigenvalue solutions are the same. In fact, writing the characteristic equation of the eigenvalue problem (14) yields Note that eigenvectors in Eq. (14) consist of left and right nodal displacements sub-vectors therefore to get an eigenvector in the form as in Eq. (11) one can apply the following transformation Incoming and outgoing waves and the wave basis. The eigensolutions of Eq. (14) can be separated into m pairs of roots, ± y,i ; these correspond to negative waves (with superscripts "−") and positive waves (with superscripts " +"). In the context of scattering properties at junctions, the negative and positive waves will also refer to incoming and outgoing waves moving towards or away from the junction, respectively, see the left-handside of Fig. 1. Furthermore, the waves can be categorised as propagating, evanescent or attenuating.
In the absence of damping, the standard wave classification 24,49,61 consists of checking whether | y,i | = 1 or not to identify propagating or evanescent waves, respectively. If the damping parameter is not zero, we establish an algorithm for the ith wave as shown in Table 2, where Re(iωφ * q,i φ f ,i ) is the energy flux of the ith wave in the positive y direction. The real parameter c in Table 2 is chosen empirically to separate the least attenuating waves from evanescent or strongly attenuating waves.
In the case of isotropic materials and some special types of composite plates, e.g. cross-ply laminates, which consist of layers with ply direction angles 0 • or 90 • , the transfer matrix Z in Eq. (11) is symplectic. That means that the eigenvalues of Eqs. (11) and (14) appear in pairs ( y,i , 1/ y,i ) , and k + y,i = −k − y,i . However, for general composites, the transfer matrix Z is not symplectic, and k + y,i � = −k − y,i . Figure 2 shows an example of a bending wave vector curve, where one can see the inequality between incoming and outgoing wave number components k y represented by red squares and blue dots, respectively. Furthermore, we note the difference between phase and group angles in composite plates shown with θ ± 1,2 and α ± 1,2 , respectively. These features were discussed for a similar wave vector curve in 20,62 .
Recall that in the absence of damping, we expect that waves are either propagating or evanescent. While the wave number components k ± y,i for propagating waves are indeed purely real, those of non-propagating waves can still have both a real and imaginary part for composites, thus seemingly producing an energy flux in the y direction 20 . However, it turns out that these waves are actually purely evanescent and decay or increase exponentially towards the line y = 0 , say, however, with the decaying/increasing direction axis not necessarily aligned with the y axis, see 63,64 . Figure 3 presents the contour plot of the real part of the out-of-plane displacement field created by the outgoing bending propagating (a) and evanescent (b) waves in a 45 • / − 45 • /45 • / − 45 • /45 • composite plate. In the case of the evanescent wave, the displacement at y = 0 along the x axis is oscillating, as expected. In contrast, at y > 0 , it decays away along the inclined null-lines, i.e., lines at which displacement is zero (black straight lines in Fig. 3b). This leads to the oscillating displacement shape projection along the y axis, which is why the energy flux along the y axis appears to be non-zero.
The angle of decay/increase of the evanescent waves depends on the angle of orientation of the principal material axes of plies of the plate, i.e. on the ply direction angle. It can be computed as arctan(Re(k ± y,i )/k x ) . Figure 4 presents the decay direction angles of longitudinal, shear and bending evanescent waves with respect to the y axis as a function of the ply direction angle of an orthotropic plate with material parameters from Table 1. It can be noted that the decay of evanescent waves can be inclined with respect to the y axis with angles up to 60 • ; for example, the line corresponding to the evanescent shear wave in Fig. 4. We can see that the decay angles are zero for ply direction of angles 0 • and ±90 • , that is, the evanescent waves decay/increase along the y axis. In such cases, the plate is specially orthotropic if it consists of only one layer 20 or balanced if it consists of multiple layers 1,56 . www.nature.com/scientificreports/ Identifying correctly whether a wave is incoming/outgoing and propagating/attenuating is not enough to further proceed with calculating scattering coefficients. In fact, for a range of frequencies and wave number components k x , we obtain a set of unclassified branches of propagating, evanescent and attenuating waves. To identify the eigensolutions of Eq. (14) corresponding to the same wave type, we apply the so-called MAC criterion 65 . For an eigenvector solution � i = φ ± q,i φ ± f ,i T defined at frequency ω for a fixed wave number component k x , we find an eigenvector solution � j = φ ± q,j φ ± f ,j T at the frequency ω + dω with sufficiently small d ω such that is maximised. The same criterion can be utilised to classify wave types for a range of wave number components k x and a fixed frequency with corresponding step size d k x being sufficiently small.   Table 1. Table 2. Properties of eigenvalues and associated waves. The expression Re(iωφ * q,i φ f ,i ) denotes the power flow of the ith wave.

Incoming Outgoing
Propagating ⇔ Re(k y,i ) > c Im(k y,i )

Hybrid FE/WFE method and scattering coefficients
In this section, we derive the energy scattering coefficients using continuity and equilibrium conditions at the boundaries of the joint.
Governing equations of the joint. We consider N composite plates that are connected together via a joint element, see Fig. 5a. The local y axes of all plates are directed away from the joint, whereas the x axes are aligned with the x J axis of the joint. The joint is assumed to be periodic in the x J axis, and the width of the periodic segment is assumed to be equal to d y . According to Fig. 5b, the nodal displacements vector of the joint is organised in the following way The nodal displacement vector q O consists of degrees of freedom of the internal nodes. It is required that the node arrangement on the face containing q LB,k , q L,k and q LT,k , k = 1, . . . , N is coherent with one on the left face of the kth plate. Therefore, we enforce that |q LB,k | = |q L,k | = |q LT,k | = |q L,k | = m k and n J,k = n x,k , where n J,k is the number of mesh cells in the joint element along the x direction of the face containing q LB,k , q L,k and q LT,k .
The nodal forces vector f J is arranged in the same way.
. . . Reflection and transmission at interfaces of the joint. Equation (23) considers only interface nodes, that is, nodes shared between the plates and the joint. The force equilibrium and displacement continuity conditions must be satisfied at these nodes. Specifically, q LB,k and f LB,k , k = 1, . . . , N denote nodal displacements and forces at nodes shared between the kth plate and the joint. Therefore, one can represent q LB,k and f LB,k in the wave mode shapes basis of the kth plate using Eq. (18) as follows where R k is the matrix that transforms the local coordinate system of the kth plate to the global coordinate system. Since the local x k axis is aligned with the global x J axis, the transformation matrix R k can be written as where ψ k denotes the angle of rotation between y k and y J . Now, we can concatenate individual expressions (26) for q LB,k and f LB,k to express q J,red and f J,red in Eq. (23) as with Inserting Eq. (28) in Eq. (23) yields The expression for S defines a scattering matrix relating the amplitudes A − and A + of incoming and outgoing waves, respectively. The dimension of the square matrix S is N k=1 m k , and the scattering coefficients have the form s nm ij (ω, k x ) , relating an incoming wave i in the plate n and a reflected or transmitted wave j in the plate m at angular frequency ω and wave number component k x . For the associated energy fluxes, we obtain the energy scattering coefficients as where (25) q J,red = T J q E , q J,red = � q LB,1 · · · q LB,N � T www.nature.com/scientificreports/ In the absence of damping, total energy must be conserved, hence the sum of the energy scattering coefficients over the outgoing modes equals one, that is,

Numerical case examples
In this section, we present two numerical case studies to show the applicability of the method discussed above. In the first example, we consider an L junction of composite symmetric cross-ply laminated plates, i.e. laminates that consist of an odd number of orthotropic layers with principal material directions alternating between 0 • and 90 • for cross-ply laminates. In the second example, an L-type junction of composite angle-ply laminated plates is considered. In this case, laminates consist of an odd number of orthotropic layers with principal material directions alternating between α • and −α • , α ∈ (0 • , 90 • ).
In all cases, we compute the energy scattering coefficients with respect to the wave number component k x for a fixed frequency f. The results for cross-and angle-ply laminated plates are compared with those obtained using the semi-analytical approach based on the solution of wave equations in the line-junction approximation 20 . We assume that the system is undamped; however, to facilitate the wave tracking process described in "Incoming and outgoing waves and the wave basis", a small damping coefficient η = 0.00001 is applied.

Cross-ply laminated plates.
A five-layer symmetric cross-ply laminated plate of the total thickness of h = 0.005 m is considered. The material characteristics of all layers are the same and given in Table 1. The lamination scheme is A periodic cell of length d x = 0.001 m and width d y = 0.001 m is modelled in ANSYS with 3 SOLID185 elements per ply, i.e. 15 finite solid elements in total. The usage of only one element in the cross-section is justified since the laminates considered are homogeneous in their plane dimensions. However, there must be at least 6-10 FE elements per wavelength to obtain accurate results. In other words, the wave numbers k ≤ 2π 10 max(d x ,d y ) can be computed accurately. One can use more elements in the cross-section to alleviate round-off errors due to truncation of inertia terms in the dynamic stiffness matrix if needed 35 .
As presented in "Computation of wave numbers and eigenvectors", solving (11) or (13) can yield propagating wave vector pairs k x , k y for a fixed frequency ω 0 or dispersion curves k y = k y (ω, k 0 x ) for a fixed wave number component k 0 x . Figure 6 presents bending, shear and longitudinal wave vector curves for a fixed frequency f 0 = 3000 Hz on the left side and the correspondent dispersion curves for a fixed wave number component k 0 x = 5 m −1 on the right side. These numerical dispersion relations can be used to calculate the group velocity angles and, therefore, propagation angles of transmitted waves via modified Snell's law.
Next, we consider two identical cross-ply composite plates connected through an L-joint. The corresponding FE model is similar to the model shown in Fig. 5, but now only connecting N = 2 plates. The FE model of the joint consists of 55 SOLID185 elements; thus, the dimensions of the joint stiffness and mass matrices are 164-by-164.
N m=1 j t nm ij = 1.  20 for various incident modes as a function of wave number component k x at frequency f = 3000 Hz. We note that the semi-analytical approach is based on further approximations, such as treating the junction as a one-dimensional line and is considered here only as a reference. Deviations between the WFE approach and the semi-analytical approach are mainly due to inaccuracies in the semi-analytical treatment. All solid and dashed lines represent numerical results, whereas circles and squares show semi-analytical results. Excellent agreement between numerical and semi-analytical energy scattering coefficients can be noted for incoming longitudinal and shear waves in Fig. 7b,c , respectively. Hence, modelling coupled thin composite plates using the semi-analytical method is sufficient for accurately estimating in-plane wave energy scattering coefficients.
For the case of an incoming bending wave, discrepancies between the WFE and the semi-analytical results appear. For instance, in Fig. 7a, it can be seen that semi-analytical results underestimate energy reflection and hence, overestimate energy transmission of an L-junction of plates. In fact, the maximum difference observed between bending using the WFE and semi-analytical results is ∼ 20 to 22%. This can be referred to the fact that the semi-analytical model results are based on the assumption that the joint can be represented as a shared onedimensional line between plates 20 . This assumption breaks down at higher frequencies since the influence of the internal joint geometry becomes significant. Notably, the shear strain becomes more critical in the dynamic response of the joint and this effect is not considered in the semi-analytical model 49,66 .  Figure 8 presents bending, shear and longitudinal wave vector curves of an angle-ply laminated plate for a fixed frequency f = 3000 Hz on the left side and the correspondent dispersion curves for a fixed wave number component k x = 5 m −1 on the right side. Note that there are two shear wave dispersion curves present on both sides of Fig. 8. The second shear wave (plotted in green) exhibit a negative group velocity phenomenon-they will be denoted as S 2 from now on. Details about these features in an orthotropic plate are discussed in 20 .
A comparison of energy scattering coefficients obtained from the current approach and from the semianalytical method is given in Fig. 9 for various incident modes as a function of wave number component k x at frequency f = 3000 Hz. Again, numerical and semi-analytical energy scattering coefficients agree well for incoming longitudinal and shear modes, see Fig. 9b,c .
Deviations in the shape of energy reflection and transmission coefficients can be seen for an incoming bending wave in the case of SOLID185-based numerical results, see Fig. 9a. This phenomenon can be explained by the fact that the semi-analytical approach produces effective one-layer plates joined along the shared edge, thus losing the complexity of the connection between individual layers of plates at the junction.
As in the case of cross-ply laminated plates, curves describing the coupling between shear and bending waves are symmetric for a range of wave number components |k x | < 7.6 m −1 -the correspondent scattering coefficients obey t

Conclusion
A hybrid FE/WFE model has been developed predicting the scattering properties for different junctions of twodimensional anisotropic composite plates. The influence of the angle of incidence on the distribution of the power flow of incident bending, shear and longitudinal type waves has been investigated. Numerical results presented were compared with semi-analytical evaluations of scattering coefficients. The method gives for the first time a detailed recipe for computing scattering coefficients for the generic case of an arbitrary number of composite plates connected at a junction without restrictions on the angles at which the plate meet or the orientation of the principal axes of individual plates. The results of this paper can be used for the computation of wave energy distributions in SEA by providing coupling loss factors 67 and for angle-of-incidence dependent scattering coefficients entering the radiative transfer and DEA method.