Abstract
Free vibration analysis of beams with single delamination undergoing bending-torsion coupling is made, using traditional finite element technique. The Galerkin weighted residual method is applied to convert the coupled differential equations of motion into to a discrete problem, where, in addition to the conventional mass and stiffness matrices, a delamination stiffness matrix, representing the extra stiffening effects at the delamination tips, is introduced. The linear eigenvalue problem resulting from the discretization along the length of the beam is solved to determine the frequencies and modes of free vibration. Both “free mode” and “constrained mode” delamination models are considered in formulation, and it is shown that the continuity (both kinematic and force) conditions at the beam span-wise locations corresponding to the extremities of the delaminated region, in particular, play a great role in “free mode” model formulation. Current trends in the literature are examined, and insight into different types of modeling techniques and constraint types are introduced. In addition, the data previously available in the literature and those obtained from a finite element-based commercial software are utilized to validate the presented modeling scheme and to verify the correctness of natural frequencies of the systems analyzed here. The paper ends with general discussions and conclusions on the presented theories and modeling approaches.
1. Introduction
Materials or components built up in layers are used in various manufacturing sectors ranging from shipbuilding to aerospace. The ever-increasing application of such layered structural elements is primarily due to their many attractive features such as high specific stiffness, high specific strength, good buckling resistance, and formability into complex shapes, to name a few. Delamination is the most common failure mode in layered structural elements and may be introduced during processing or subsequent service conditions, leading to the reduction or loss of structural strength and stiffness [1], which, in turn, will modify the vibration characteristics of the structure, i.e., natural frequency and mode shape. Natural frequencies will decrease caused by the stiffness reduction, which may subsequently lead to resonance, should the reduced frequency match the working frequency. Therefore, it is imperative to understand the effect of the delamination and its influence on the vibration characteristics of structures. Depending on the number of delamination, size, and location, the defective sublaminate generally displays different vibration frequencies and modes.
One of the earliest models for vibration analysis of composite beams with delaminations was proposed by Kulkarni and Fredericks [2] where they studied a one-dimensional (1D) problem with a single delamination. They investigated a defective, circumferentially symmetric (circular) cylindrical shell, with a small length delamination located at the middle surface. In their analysis, it was assumed that the flexural rigidity of the delamination region can be taken as the sum of the flexural rigidities of the all individual delaminated layers. This analysis resulted in system natural frequencies considerably lower than those found experimentally. Ramkumar et al. [3] investigated the vibration characteristics of defective composite beams, where a single through-the-width delamination was modeled as four Timoshenko beams connected at delamination edges. The predicted natural frequencies were consistently found to be lower than those measured experimentally. This discrepancy was attributed to the contact between the delaminated free surfaces during vibrations, and the authors suggested that the inclusion of the contact effect might improve the analytical prediction. However, it was later discovered that underpredicted natural frequencies for off midplane delaminations, resulting from the free mode model, are in fact due to unrestricted penetration of the beams into each other. Considering free vibration characteristics of isotropic defective beams with through-the-width split, Wang et al. [4] developed an improved theoretical model to investigate the effect of delamination on the system’s frequencies. Assuming each beam segment, modeled as an Euler beam, vibrates freely without touching each other, they included the coupling between flexural and axial vibrations of the delaminated sublaminates in their model. Although this improved analytical model resulted in frequency data closer to experimental results, it predicted a few physically inadmissible opening modes, addressed by Mujumdar and Suryanarayan [5]. The “constrained mode” model [5] assumes that the delaminated beam segments undergo identical transverse displacements and predicts vibration behaviour more accurately for off midplane delamination configurations. However, opening delamination modes, commonly seen in experimental analysis [6, 7], cannot be captured using the constrained model. Tracy and Pardoen [8] also used the constrained model to assess the effect of delamination on natural frequency of symmetric laminated Euler beams containing midplane delamination. Based on both the classical and high-order shear deformation beam theories, and using the FEM, Nagesh Babu and Hanagud [9] studied the effects of delamination on the system’s natural frequencies. Performing experiments, Shen and Grady [7] observed the presence of opening modes, undetected by the constrained model. By exploiting the Galerkin method, they also investigated the effects of delamination on the natural frequency and mode shape of composite laminated beams. Introducing coupling between longitudinal and bending motion, Stamos et al. [10] presented a delaminated composite beam model and furthermore suggested an inverse method to determine position and size of delamination based on the degradation of the first two natural frequencies. Lee et al. [11] analyzed the free vibration of multi-delaminated composite beam-columns subjected to axial compression load analytically and experimentally. In their work, the characteristic equation of the defective system is obtained by dividing the multi-delaminated beam-columns into segments and by imposing recurrence relation from the continuity conditions on each subbeam-column. To verify their results, experimental results are obtained for isotropic single delaminated beam-columns and confirmed that the sizes, location, and number of delaminations have significant effect on the system’s vibration and stability. Chen et al. [12] developed an analytical model for free vibration of a delaminated layered composite beam in the prebuckled state. They also presented a new constrained model, including both compressive axial force and bending-extension coupling effects. Saravanos and Hopkins [13] investigated the damped-free vibration of unstressed defective laminated composite beams, both analytically and experimentally, where they considered the kinematic perturbations induced by a delamination crack as additional degrees of freedom. In a later work, Chrysochoidis and Saravanos [14] assessed the effects of delamination on the dynamic response of composite beam specimens with surface attached piezoelectric actuators and sensors using experiments and same analytical method developed earlier in [13]. Later, they developed an analytical model using coupled linear layerwise laminate theory and a FE model to investigate both static and dynamic response of flexible defective composite beams, including delaminations and active piezoelectric sensors [15]. Della and Shu [16, 17] reported an analytical solution method, and more recently, Erdelyi and Hashemi used Dynamic Stiffness Matrix (DSM) [18] and presented a FEM-based novel assembly technique [19] to investigate the behaviour of a delaminated slender beam and compared their results with those reported by Wang et al. [4] and Lee [20]. In a recent work, Szekrényes [21] performed stability analysis on delaminated composite beams undergoing coupled flexural-longitudinal vibration. Starting with Timoshenko beam theory and then reducing the model to Euler–Bernoulli solution, they presented a formulation, considering both free and constrained models. In addition, the equality of axial forces in the top and bottom beams was derived and shown in an exact way. Also, the continuity of the effective bending moments was related to the equilibrium equations, and it was also concluded and reported that if the normal force is compressive in one of the half-periods of the vibration and reaches a critical value, then it can lead to delamination buckling [21–23].
In two recent publications, the authors exploited the classical [24] and frequency-dependent (dynamic) FEM [25] formulations for the free vibration of intact prestressed beams. The free vibration/buckling of intact and delaminated axially-loaded beams has also been investigated before (e.g., [17–23]). However, to the authors’ best knowledge, the vibration of prestressed delaminated beams, subjected to combined axial load and end moment, has not been reported in the open literature.
In the recent years, layered, sandwich and composite elements have found their place and have been integrated in some commercial software, such as ANSYS®, used to perform static and dynamic structural analysis. However, it is worth noting that modeling of delaminated configurations in such environments in not a straightforward procedure can involve cumbersome, complex, time-consuming, and error-prone processes. It requires manual model creation, involving the use of, for example, Multipoint Constraint Rigid Link (ANSYS element MPC184) to enforce the displacement and slope continuity at the edges of delamination region [26]. Therefore, in the present study, the free vibration behaviour of a delaminated prestressed beam subjected to combined axial force and end moment is investigated using the finite element method (FEM). By implementing the Galerkin weighted residual method, the governing differential equations of motion are discretized, leading to element matrices. As it will be explained later in this paper, in addition to the mass and stiffness matrices commonly resulting from the conventional FEM formulation, the presented procedure also leads to three extra stiffness matrices: a geometric stiffness matrix caused by the axial load, a delamination stiffness matrix, resulting from the relationships between the two delamination region’s extremities, and a coupling geometric stiffness matrix caused by the applied end moment. The assembly of the element matrices and application of relevant boundary conditions then leads to the system’s linear eigenvalue problem. The resulting eigenproblem is then solved to extract the system’s natural frequencies and modes. The correctness of the proposed FEM method is first validated, where the free vibration analysis of a simple delaminated beam (in the absence of prestress) is conducted, and the frequency results are compared with those presented in the literature. The application of the presented model and formulation is then further demonstrated through the illustrative examples of unstressed delaminated [09/0]2s composite beams, as well as bending-torsion coupled beams, subjected to various prestress load combinations (i.e., axial forces and end moments). The paper concludes with a brief discussion on the presented method and numerical results. While the FEM model presented in this paper assumes isotropic materials, further research is underway to extend it to sandwich and fibre-reinforced laminated composite beams, characterized by an extensional response coupled with flexural/torsional and coupled bending-torsion vibration.
2. Mathematical Model
Figure 1 shows a prestressed two-layered beam of length and thickness having a single through-the-width delamination with a length , starting from . This delamination divides the beam into four segments, which are assumed to have a slenderness ratio of greater than 10 and, as a result will be analyzed as four interconnected Euler–Bernoulli beams. Beam segments 2 and 3 have and thicknesses, respectively. The structure is subjected to an axial force of and constant end moment of , applied at its ends, and .
Assuming that each uniform beam segment is made of linearly elastic, homogeneous, isotropic material with constant mechanical, material, and geometric properties, it can be shown that the partial differential equations of motion, governing free, linear, bending-torsion coupled vibrations of the Euler–Bernoulli beam segment subjected to axial force and end moment are as follows [24]:where is the lateral displacement along direction, stands for the torsional twist about x-axis, is Young’s modulus, represents the shear modulus, and are, respectively, the torsion constant and the polar moment of inertia, and subscript represents the beam segment’s number, where . and represent axial force and end moment for beam, respectively. A prime, (′), denotes spatial derivative with respect to , and the over dot, (•), stands for time derivative. The beam’s torsional rigidity () is assumed to be very large compared with its warping rigidity (), and ends are free to warp, i.e., state of uniform torsion.
As can be observed from Equations (1) and (2), the system is coupled by the end moment, . Exploiting the simple harmonic motion assumption to separate variables in space and time, and introducing the (sinusoidal) functions and for lateral and torsional displacements, respectively, in expressions (1) and (2) leads to following form of the governing differential equations, in frequency domain:where denotes the frequency and and are the amplitudes of flexural and torsional displacements, respectively, of the beam segment .
2.1. Free Mode Delamination Model
In contrast to the “constraint mode” model, the “free mode” delamination scheme considers no constraint or interaction between the upper and lowest layers in the defective zone, i.e., segments 2 and 3. The “free mode” model predicts the behaviour of the beam more accurately but it also exhibits some physically impossible modes of vibration. To solve the four pairs of coupled equations as one system numerically, the global end conditions at both extremities of the beam, as well as continuity requirements at the delamination tips, are to be enforced. The continuity conditions for lateral displacement, slope, torsional displacement, and axial force at the left side of the defective region () are expressed as
Similarly, continuity conditions for the right delamination tip are written as
Additionally, the requirement for the delamination tip faces to remain planar after deformation, at the left delamination tip, results inwhere is the axial displacement of beam section . Combining expression (10) with the similar expression from the right delamination tip leads towhere
It was originally assumed by Mujumdar and Suryanarayan [5] and was later shown in an exact way by Szekrényes [21] that for small deformations of a uniform beam, the axial displacement will behave according to the following expression (13):where represents the local position along the length and is the axial force induced in the segment . Compressive force during vibration, , consists of applied static axial force, , and perturbed axial force, . The latter, preventing interlaminar slip, is induced in the delaminated layers during vibration. Knowing that the change in axial displacement for each segment is caused by induced perturbed component of axial force rather than its static component, substituting Equation (13) into (11) yields
Taking the perturbation of Equation (8) for continuity of axial force on the left delamination tip, one obtains
Since and is constant, which means: . Substituting into Equation (14) results inwhere the coefficient is defined as
At the left and right delamination tips, respectively, continuity of bending moments leads towhere . The internal bending moment, , shear force, , and torque, , in each of the individual beam segments are written as
Using Equations (20) and (18), the following continuity of bending moments is obtainedwhere the parameter is defined as
Likewise, to satisfy the continuity of shear forces at the left delamination tip, one should haveand finally, the continuity of twisting moment at left delamination tip yields
Using Galerkin weighted residual formulation, the integral forms of the governing differential equations (3) and (4), also representing the virtual flexural () and torsional () works are written aswhere and are weighting functions corresponding to flexural and torsional displacements, respectively, for beam segment . Then, performing integrations by parts twice on Equation (27) and once on Equation (28) leads to
Considering expressions (20)–(22), it can be shown that for all the system’s global classic end conditions (at and ), the corresponding relationships in boundary terms () and () will vanish. For example, for clamped-free boundary conditions, one has zero displacements and zero virtual displacements at the clamped end () and null resultant shear force, , bending moment, , and twisting moment, at the free end (). This leads to
The remaining terms in (), corresponding to delamination edges (i.e., and ), can be resolved by applying the continuity conditions (23) and (25), with the following as a result:
The terms () and () go to zero directly because of shear force continuity conditions. The following remaining terms, however, do not vanish:
Likewise, the remaining terms in () can be resolved by applying the continuity conditions (7, 11):
The terms () and () vanish as a result of twisting moment continuity conditions for left delamination tip (Equation (26)) and similar continuity condition for right tip. With the end conditions and continuity requirements satisfied and knowing that the expressions (27) and (28) must satisfy the virtual work principle, i.e., , the system discretization leads towhere stands for the number of elements in beam segment . Following the traditional Euler–Bernoulli FE development, cubic Hermite shape functions for flexural displacement and linear shape functions for torsion displacement are used to relate the displacement within each element () to their corresponding nodal displacements:where, and are row vectors of cubic shape functions and linear shape functions, respectively. , , , and are vectors containing the nodal real and virtual displacements for each element (), defined as
Substituting the above approximate displacements in (34),where and are real and virtual nodal displacement vectors for the two nodes at the two ends of the delamination, while and are the elemental nodal vectors of real and virtual displacements (36) and (37), rewritten as
Assembly of terms in (38) and enforcing the system’s global boundary conditions leads to the following linear eigenvalue problem:where and and are the global bending and torsion stiffness matrices, respectively. and are global bending-torsion and torsion-bending stiffness matrices, respectively, resulting from coupled terms in bending and torsion equations. is delamination stiffness matrix generated from the term outside the integral in Equation (38), and and represent the overall (global) mass and stiffness matrices, respectively. Finally the system’s natural frequencies and corresponding modes are the system’s eigenvalues and eigenvectors evaluated from eigenproblem (40), i.e., for arbitrary , . The above procedure is achieved through a code developed in Matlab®.
2.2. Constrained Mode Model
In the simplified “constrained mode” model, the delaminated layers are “constrained” to have the same transverse deformations. The delaminated beam configuration is analyzed as assembly of three beam segments I, II, and III, as shown in Figure 1, where 0 ≤ x ≤ x1 represents segment I, x1 ≤ x ≤ x2 constitutes segment II, and x2 ≤ x ≤ L is segment III.
For beam segment I, the governing equations are
For segment II,and for segment III,
The continuity conditions for lateral displacement, slope, shear force, bending moment, and torsion torque at the left edge of delamination, , are
The continuity conditions at the right delamination tip are also written in a similar way to expressions (44)–(49). The general solution steps, namely, the assembly of element equations, application of boundary conditions, forming the system’s linear eigenproblem, and extraction of the natural frequencies and modes for the “constrained mode” are identical to those stated above for the “free mode” model.
3. Numerical Results
In this section, numerical tests performed to confirm the predictability and accuracy of the proposed FE models for delaminated beams, are presented. A Matlab® code was developed to create the model, assembly, and the application of the system boundary conditions, and the resulting linear eigenproblem was solved using the “eig” function.
In the first two examples, the natural frequencies of defective steel and graphite/epoxy beam configurations, in the absence of axial load and end moment (i.e., ), are investigated to validate and verify the applicability of the formulation against existing data. The systems are assumed to have a central split, about the midsection (L1 = L4) and of various lengths up to 60% of the span (0 ≤ L2/L ≤ 0.6). The split is occurring symmetrically along the midplane of the beam and surrounded by intact beam segments, as has also been presented and studied in [4, 18]. The defective FEM models are then created and used to evaluate the natural frequencies and mode shapes of various prestressed delaminated beam configurations. As the benchmarks for comparison and validation purposes, the results from references [4, 5, 7, 13, 18] were used.
3.1. Validation of Presented Formulation
3.1.1. Steel Beam
Consider a delaminated steel beam with elastic modulus of E = 200 GPa, mass density of ρ = 7800 kg/m3, length of 8 m, and a rectangular cross-sectional area of width of 0.4 m and depth of 0.2 m. The defect is assumed to be a midplane (H2 = 0.5H1) delamination of variable length of 0 ≤ L2/L ≤ 0.6 (refer to Figure 1). A sample of convergence study is presented in Figure 2, where the analytical results presented by Wang et al. [4] are used as the benchmark to evaluate the performance of the presented formulation. The un-prestressed defective beam, in this case, is assumed to be clamped- clamped and isotropic, with P0 = 0, Mzz = 0, H2/H = 0.3, and L2/L = 0.4. To validate the present FE solution method, the first two nondimensional natural frequencies, λ2, of a defective clamped-clamped beam with a through-the-width delamination occurring symmetrically about the midsection (L1 = L4) on the midplane (H2 = H3) for various delamination lengths are compared with the analytical results reported by Wang et al. [4], as well as the DSM [18] and FEM data [19] by Erdelyi and Hashemi, where
In the formulations presented by Wang et al. [4] and Erdelyi and Hashemi [18, 19], no prestress effect and torsional vibration are investigated. Therefore, for comparison purposes, only bending equation is solved and applied axial force and end moment are set to zero; P0 = 0 and Mzz = 0. Erdelyi and Hashemi [19] presented FEM-based defective beam models and the frequency data calculated from six- and ten-element meshes of 2-node FEM elements. Tables 1 and 2 summarize the system’s first two natural frequencies obtained using the presented method, in comparison with those reported in the literature [4, 6, 17–19]. Referring to Tables 1 and 2 of reference [18], they showed exact match between the 1st frequency values obtained from both meshes. The 2nd frequency values obtained from the 6- and 10-element meshes, however, were found to be slightly different. Therefore, in what follows, frequency data obtained from the latter mesh are used for comparison.
As can be observed from Tables 1 and 2, the natural frequencies calculated from the present method are in excellent agreement with the analytical results reported by Wang et al. [4], Della and Shu [16, 17], FEM [19], and layerwise FEM [6] data, with discrepancies less than 0.1%.
3.1.2. Laminated ([0/90]2s) T300/934 Graphite/Epoxy Cantilever Beam
For the second validation case, consider an unstressed, defective, cantilevered, laminated beam ([0/90]2), made with T300/934 graphite/epoxy, containing a single, variable-length central delamination, located at the midplane, as also reported by Saravanos and Hopkins [13]. The dimensions of the beam are the same as the one investigated in the experimental study by Shen and Grady [7], i.e., 10 × 0.5 × 0.0435 inches. Since the presented model in this paper is based on isotropic materials, the classic lamination theory (CLT) [27] is used to find the properties of equivalent single layer beam for each of the four segments in Figure 1. In addition, since the beam is unstressed, the axial load and end moment are set equal to zero. As can be observed from Figure 3, the FEM results obtained from the present theory show excellent agreement with both analytical and experimental data reported by Saravanos and Hopkins [13] and Shen and Grady [7]. It is to be noted that the analytical and experimental frequency data are extracted from Figures 4 and 5, reported in references [7, 13], respectively.
3.2. Vibration Analysis of Prestressed Delaminated Beams
To further investigate the validity and practical applicability of the proposed formulation and the effect of axial load on the frequency response of the defective beam structures, following illustrative examples of prestressed symmetrically delaminated beams (Figure 1), characterized by various delamination lengths, applied axial force and end moments are investigated. In each case, the change of normalized fundamental frequency, , versus the normalized end moment, , compressive axial load, , for “free mode” is investigated, wherewhere and , respectively, stand for the critical buckling load and buckling moment of an intact (with no delamination) beam. Figure 4 shows the system’s fundamental natural frequency, , of a homogeneous beam with clamped-clamped boundary condition, , and central delamination located in midplane , versus normalized axial force in “free mode” model for different delamination sizes, . As can be observed, when the applied axial load, , approaches the system’s buckling load, the natural frequency, , turns into zero, i.e., the structure has buckled. Moreover, as expected, the larger the delamination size, the smaller the and the fundamental frequency. This reconfirms the fact that structural stiffness decreases when the delamination size increased.
Figure 6 demonstrates the change in the first and second natural frequencies with respect to normalized buckling load P0/Pcr for Mzz/Mb = 0.4 and H2/H = 0.3. Results are shown for different delamination lengths of L2/L = 0.2, 0.3, 0.4, and 0.5. The figure reveals similar trends to those observed in Figure 2. Fundamental frequency results of “constrained mode” are compared to “free mode” data in Figure 7. As can be observed, the difference between predicted frequencies from these two modes increases as the applied axial load increases. But for both natural modes, the larger the compressive axial load, the smaller the frequency. It can be observed that as P0 increases, the natural frequency λ2 first decreases slowly, but by getting close to buckling load it decreases drastically and eventually goes to zero. As can also be seen from Figure 7, the “constrained mode” comparing to “free mode,” overpredicts the natural frequencies especially for compressive load values close to buckling load. Additionally, from the same trends observed in Figures 2–4, it can be inferred that compressive axial load causes a reduction in the structural stiffness leading, in turn, to lower natural frequencies. Similarly, increasing the delamination length results in lower stiffness and natural frequencies.
Figure 5 presents the normalized natural frequency, λ2, versus the normalized end moment, Mzz/Mb, for both constrained and free modes with H2/H = 0.5, P0/Pcr = 0.4, and different values of L2/L = 0.3, 0.4, and 0.5. Referring to Figure 5, the same trends as in the previous graphs (for axial force) is observed for variation of natural frequency versus buckling moment. The difference between predicted frequencies from these two modes increases as the applied end moment increases, and the larger the end moment, the smaller the frequency. It can also be observed that as Mzz increases, the natural frequency λ2 first decreases slowly, but by getting close to buckling load it decreases drastically and eventually goes to zero. In addition, the “constrained mode” comparing to “free mode,” overpredicts the natural frequencies. Once again, it can also be seen that as the end moment increases, the stiffness decreases, leading to lower natural frequencies.
Figure 8 depicts the mode shape associated with fundamental natural frequency of an intact beam in comparison with the first opening mode of vibration for a defective beam, characterized by a central delamination on the midplane, obtained using “free mode” model. This mode, however, cannot be detected using the constrained model. As also reported in the literature [19, 22], not all the mode shapes obtained from “free mode” model are admissible since unlike “constrained mode” no constraints are imposed between the two delaminated layers, which results in some physically impossible mode shapes. It is also worth noting that the experimentally observed delamination openings have been reported to be much fewer than those suggested by the free model, which would rather justify the use of the constrained model for further analysis [21].
Finally, it has been recently shown and experimentally verified by Szekrényes [22] that the opening is amplitude dependent. In other words, the delamination opening takes place only at certain critical amplitudes. The mode shapes are asymmetric and can be approximated by the superposition of the global shape of the entire beam and the local buckling eigenshape of the delaminated part based on a dynamic stability analysis [21]. Further investigation on this topic, however, is beyond the scope of this paper.
4. Discussion and Concluding Remarks
A systematic FEM-based formulation and numerical solutions for the free vibration modeling and analysis of delaminated beams under axial compressive load and end moment were presented. The delaminated beam configurations are modeled and analyzed as the assemblage of four Euler–Bernoulli beams attached at the defect edges. Both the “free mode” and “constrained mode” models were created to carry out buckling and vibration analyses, and the results were validated against the analytical data and other numerical results available in the open literature. Parametric studies were also performed to investigate the effects of end moment and axial load on the natural frequencies and mode shapes of delaminated beams. The results exhibited a monotonic relation between the system’s natural frequency and both loads (i.e., axial compressive force and end moment). Both compressive axial load and end moment reduce the structural stiffness and result in lower predictions for natural frequencies. Similarly, increasing the length of delamination results in lower stiffness and natural frequencies. Based on the “free mode” model, the first and second mode buckling loads of the delaminated beams were also studied. In general, the “constrained mode,” when compared to “free mode,” was shown to overpredict the natural frequencies and specifically when compressive load approached to values close to buckling load. Finally, the difference between predicted frequencies from “free mode” and “constrained mode” delamination models increases as the applied axial load or end moment increases.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Disclosure
This paper presents the results of a recent research conducted by the first author (Mir Tahmaseb Kashani), under the supervision of the second author (Seyed M. Hashemi).
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
The supports provided by Ryerson University, the Natural Sciences and Engineering Research Council of Canada (NSERC), and Ontario Graduate Scholarship (OGS) are acknowledged.