Multiscale Wavelet Finite Element Analysis in Structural Dynamics

2.1. Daubechies wavelet

Daubechies wavelets are compact supported orthonormal wavelets developed by Ingrid Daubechies and for order L, the scaling and wavelet functions are described by the ‘two-scale’ relation [8]:

The scaling and wavelet functions have the supports 0L−1and 1−L2L2respectively. The normalised wavelet function filter coefficients qLkand scaling function filter coefficients pLkhave the relation qLk=−1kpL1−k. The multiresolution scaling and wavelet basis functions corresponding to the subspaces Vj and Wj are defined as:

The scaling and wavelet functions defined in Eqs. (10)–(13) satisfy the following properties [8]:

Certain wavelet families have no explicit formulation, as is the case with the Daubechies wavelets. Therefore, Eq. (10) gives rise to a system of equations that require a normalising equation obtained from Eq. (14) to evaluate the scaling functions. The Daubechies wavelet of order L has L2−1vanishing moments from property (18) and consequently the scaling functions at scale j can represent a polynomial of order xmwhere 0≤m≤L2−1, i.e., [37]

The coefficients Mkj,mdenote the moments of the scaling function and it translates at Vj. The derivatives of the Daubechies wavelet scaling functions are evaluated by differentiating the refinement Eq. (10) m times, and are obtained as [12]:

A normalising condition is required to evaluate Eq. (20) which is obtained from the moments of the scaling functions.

2.2. Daubechies connection coefficients

As earlier mentioned, the Daubechies functions cannot be computed analytically and their derivatives are highly oscillatory, particularly at low wavelet orders and/or high order derivatives. Therefore, the integral of the products of the scaling functions and/or derivatives are computed as what is commonly known as connection coefficients [37]. There are two forms of connection coefficients that are of relevance to this study; the multiscale two-term connection coefficient Γa,bk,lj,d1,d2and multiscale connection coefficient Υkj,m. We define the two-term connection coefficient [30]

where a and b are the orders of the scaling function at multiresolution j, while the values d1and d2denote the order of the derivative of the scaling functions. X0,1x=10≤x≤10otherwiseis the characteristic function. The formulation presented is a modified algorithm of that described in [4] and allows for the evaluation of the connection coefficients for different values of a and b at different multiresolution scales j. From the ‘two-scale’ relation presented in Eq. (10),

Differentiating Eq. (23) m times

Substituting Eq. (24) into Eq. (22) and applying the ‘two-scale’ relation of the characteristic function, the two-term connection coefficient can be expressed as:

where 2−a≤k,r≤2j−1and 2−b≤l,s≤2j−1. Eq. (25) can be expressed in matrix form as:

where the square matrix Pa,bcontains the filter coefficients as expressed in Eq. (25) and Γja,bcontains the connection coefficients. To uniquely determine the connection coefficients, normalising conditions are required to generate a sufficient number of inhomogeneous equations via the multiscale moment condition from Eq. (19)

Defining the second form of the connection coefficient

Substituting Eq. (27) into (28)

However,

Thus

where ΓL,Lk,lj,0,0are the two-term connection coefficients with a=b=Land d1=d1=0and Mlj,mare the moments earlier described.

2.3. B-spline wavelets on the interval [0,1] (BSWI)

The BSWI are a family of wavelets that emanate from Basis splines functions (B-Splines) and the basic functions in subspace Vjof order m and scale j>0are expressed as [20]

with the knot sequence

tkjtk+1j…tk+mjt, is the m^th divided difference of the truncated power function t−x+m−1with respect to variable t. The general B-splines take the form

and have support Bm,kjxsupp=tkjtk+mj. The B-spline basis function has simple knots inside the unit interval and m-tuple knots at 0 and 1, as expressed in Eq. (33). The knots at 0 ad 1 coalesce and form multiple knots for BSWI while the internal knots are simple hence smoothness is unaffected. For the knot sequence on [0,1], tkjis given as [38]:

The number of inner scaling functions present in the formulation of BSWI is determined by the scale j. There must be at least one inner scaling function on the interval [0,1] and this gives rise to the minimum value of j necessary to ensure this condition is met and is defined as j0:

The basis Bm,kjxfrom the inner knots corresponds to the m^th cardinal B-splines, Nmx, at multiresolution j [38]:

where ϕm,kjxis the BSWI scaling function which can be differentiated m times. The corresponding B-wavelet with support suppψm,kjx=k2jk+2m−12jis expressed as:

B2m,kj+1,mxis the m^th derivative for the B-spline of order 2m and scale j+1and can be evaluated explicitly from Eq. (34). Given that the requirement j>j0ensures at least one inner B-wavelet is present, the scaling and wavelet function of the BSWI are obtained as [39]:

and the scaling function derivatives can be evaluated directly by differentiating Eq. (40).

3.1. Axial rod wavelet finite element

Assume each WFE is divided into equal segments, ns, connected by r=ns+1elemental nodes, as shown in Figure 1, with axial deformation ui. The total number of degrees of freedom (DOFs) within each WFE is denoted by n=rfor n,r∈N. Vector ue=u1u2⋯ur−1urTcontains all the axial DOFs in physical space, as illustrated in Figure 2(a), where ui=uxirepresents the elemental node axial deformation DOF at node i corresponding to coordinate position xi. The nodal natural coordinates is ξi=xi−x1Le(0 ≤ ξi ≤ 1, 1 ≤ i ≤ r). The Daubechies and BSWI scaling functions ϕz,kjxare used as the interpolating functions and for a family of order z at multiresolution scale j, the axial deformation

contains the unknown wavelet coefficients az,kj. This gives rise to the vector uecontaining the axial deformations at all elemental nodes in physical space.

The matrix Rrw=Φzjξ1Φzjξ2⋯Φzjξr−1ΦzjξrTcontains the scaling function vectors Φzjξiapproximating the axial deformation at the corresponding elemental nodes and ae=az,hjaz,h+1j⋯az,2j−2jaz,2j−1jT. The axial deformation at any point along the rod element can be generalised as:

The matrix Trw=Rrw−1is the axial rod wavelet transformation matrix with the scripts r and w denoting rod and wavelet respectively. The wavelet based axial rod shape functions can be evaluated as Nr,eξ=ΦzjξTrwwithin each element.

Suppose the axial rod is subjected to nodal point loads fxiand distributed loadingfdx, then the potential energy within the axial rod Πacan be generally expressed as [40]:

where E is the Young’s modulus, A is the cross-sectional area and l is the length of the rod. Therefore, given the relation highlighted in Eq. (44), the axial stain energy Ueawithin each WFE of length Leis expressed in natural coordinates as:

The stiffness matrix of the rod element in wavelet space, kr,ewis computed using the first derivative of the scaling functions and is symmetric.

In order for one to obtain the stiffness matrix in physical space, the element properties and transformation matrix Trware applied to the wavelet space stiffness matrix in Eq. (47).

The load vector containing the axial point loads of the WFE in physical space is obtained as:

and the equivalent nodal load vector for the distributed load fdxin physical space is

When applying the Daubechies wavelet family, the WFE has a total of n=2j+L−2DOFs. The wavelet space stiffness matrix is evaluated from the multiscale two-term connection coefficients Γa,bk,lj,d1,d2a=b=Land d1=d1=1and is given as:

where 22jis the normalising factor and the matrix Γj,1,1has the entries ΓL,Lk,lj,1,1for the limits 2−L≤k,l≤2j−1. Similarly, the distributed forces acting on the element require the form Υkj,mfor limits 2−L≤k,l≤2j−1of connection coefficients and the value of m depends on the order of the function fdxof the forces. In the case of the BSWI formulations, the total DOFs is n=2j+m−1and the condition j≥j0must be satisfied. Therefore, the wavelet space stiffness matrices of the BSWI axial rod are computed as:

3.2. Euler Bernoulli beam wavelet finite element

According to Euler Bernoulli beam theory, it is assumed that the shear deformation effects are neglected because before and after bending occurs, the plane cross-sections remain plane and perpendicular to the axial centroidal axis of the beam. The beam WFE of length Le, is divided into ns equally spaced elemental segments connected by r elemental nodes at coordinate values xi∈x1xrand i∈Nas illustrated in Figure 1. The WFE has the transverse displacement v and rotation θ taken into account, with corresponding transverse forces fyand moments m´respectively. The transverse displacement and rotation DOFs must be present at each elemental end node to ensure inter-element compatibility [4, 5, 6]. However, the DOFs at the internal elemental nodes can be tailored according to the desired requirements and this in turn will affect the total number of elemental segments and nodes present in each element. In this case the internal WFE nodes only have the transverse displacement present and the total number of DOFs within each beam element is n as illustrated in Figure 2(b). Therefore, there are n−2displacement DOFs and 2 rotation DOFs in total for each WFE and consequently r=n−2elemental nodes and ns=n−3elemental segments. Let the vector {ve}=v1θ1v2v3⋯vr−2 vr−1 vrθrTdenote all the physical DOFs within the beam element. The displacement and rotation DOFs corresponding to coordinate position xi∈x1xri∈Nand1≤i≤rin local coordinates are denoted as vi=vxiand θi=θxi. The nodal natural coordinate ξi=xi−x1Le(0 ≤ ξi ≤ 1, 1 ≤ i ≤ r). The deflection and rotation at any point of the wavelet based beam finite element can be approximated by applying the wavelet scaling functions ϕz,kjxof order z at multiresolution scale j as interpolating functions.

Therefore, the DOFs present within the entire beam element can be represented as

Rbw=Φzjξ11LeΦz′jξ1Φzjξ2⋯Φzjξr−1Φzjξr1LeΦz′jξrTand vector becontains the unknown wavelet coefficients bz,kjrepresenting the beam wavelet space DOFs.

From Eq. (54), the transverse displacement and rotation at any point of the beam element can be expressed as:

where Tbw=Rbw−1is the beam wavelet transformation matrix which is used to obtain the wavelet based shape functions for the beam Nb,eξ=ΦzjξTbw. The potential energy Πbwithin a Euler Bernoulli beam subjected to concentrated forces fyi, distributed force fdxand bending moments m´ican be generally expressed as [40]:

where E is the Young’s modulus, I is the moment of inertia and l is the length of the beam. The strain energy Uebwithin each beam element of length Lecan expressed in terms of the approximation of the transverse displacement via scaling functions as highlighted in Eq. (55).

This gives rise to the beam WFE stiffness matrix in wavelet space

The vector Φ′′zjξ=ϕz,h″jξϕz,h+1″jξ⋯ϕz,2j−2″jξϕz,2j−1″jξcontains the second derivative of the scaling functions. Taking into account the material properties of the beam, the wavelet space stiffness matrix is transformed into physical space via the transformation matrix Tbw.

The transverse kinetic energy of the beam element is expressed as

where v̇ξ=∂vξ∂t, ρis the density and A is the cross-sectional area of the beam. Applying the scaling functions to approximate the displacements within the beam, the kinetic energy becomes

The mass matrix in physical space of the Euler Bernoulli beam element, mb,ep, can be evaluated as:

The vectors containing the element concentrated point loads, bending moments and equivalent distributed loads in physical space respectively are subsequently evaluated as:

In various engineering problems, the loading conditions analysed vary in location and/or magnitude with respect to time, e.g., a train travelling over a track, and this is generally referred to as moving load problems. Assume a moving load of magnitude P travels across a beam element, as illustrated in Figure 3, from the left at a constant speed of c m·s⁻¹ and is represented by the function xt=Pδx−x0[41]. δxis the Dirac Delta function and x0is the distance travelled by the moving load at time t. The potential work of the load at this instant at position ξ0=x0Lein natural coordinates is [1, 30]:

Therefore, the element load vector in physical space is evaluated as

Assuming the moving load transverses to a new position ξ0within the same WFE, the numerical values of the shape functions, and consequently load vector, will change accordingly. All other WFEs representing the system with no loading present have zero entries within the load vectors at that particular time t. When the moving load is acting on a new WFE, the scaling functions corresponding to the WFE subjected to the moving load are used to obtain the load vector for that particular element.

When applying the Daubechies wavelet family of order L at multiresolution j, the total DOFs within a single element is n=2j+L−2and for this specific layout, the total number of elemental nodes is r=2j+L−4and corresponding elemental segments ns=2j+L−5. The Daubechies wavelet space stiffness and mass matrices of the Euler Bernoulli beam WFE are obtained from the connection coefficients and are expressed as:

Correspondingly, the connection coefficients of the form Υkj,mfor 2−L≤k≤2j−1are used to evaluated the distributed loads and the value of m is based on the load function fdx. For the BSWI family of order m and at scale j, there are n=2j+m−1total DOFs, r=2j+m−3elemental nodes and ns=2j+m−4elemental segments within the each WFE for this layout. The stiffness and mass matrices in wavelet space can be evaluated directly and are obtained as:

3.3. Transversely varying functionally graded Euler Bernoulli beam wavelet finite element

Functionally graded materials are a recent evolution of composite materials where the material constituents, hence properties, vary continuously in the desired spatial directions. The need for such revolutionary materials arose to overcome limitations of conventional composite materials, for instance, desirable properties would diminished when applied to highly intense thermal environments or material debonding due to increased stress concentration at material interfaces [42]. In the formulation of the wavelet based functionally grade beam as presented in Figure 4(a), of height h, length l and width b, the material distribution is modelled based on the power law of transverse gradation [43]

As illustrated in Figure 4(b), the transverse variation of the effective material properties P(y) (Young’s modulus) can be infinitely altered via the non-negative volume fraction power law exponent, n. Pratiois the ratio of the upper and lower surface material properties Puand Plorespectively.

The beam WFE has axial deformation uiand transverse deflection viDOFs at all elemental nodes and rotation θiDOFs only present at elemental end nodes with corresponding axial forces fxi, transverse forces fyiand bending moments θias illustrated in Figure 4(c). The wavelet scaling functions are implemented as interpolating functions and the axial deformation, deflection and rotation at any point of the beam element are described by Eqs. (42) and (53) respectively. However, in order to ensure that the defined DOFs are positioned correctly, the layout of the element determines the order of scaling functions selected. In this case, the order of the scaling functions selected to approximate the axial displacement is z−2if the scaling function order approximating the bending DOFs is z. The vector containing the total number of DOFs, s, present in the functionally graded beam element is he=u1v1θ1u2v2u3v3⋯ur−1vr−1urvrθrTand subsequently

where the vector cecontains the unknown wavelet space element DOFs and

Therefore, the DOFs present within the entire beam element can be represented as

and consequently

The wavelet transformation matrix Tpw=Rpw−1. The strain energy of the functionally graded beam element, Ue, is defined as

where Leis the length of the element and Eythe effective Young’s modulus obtained from Eq. (70). Let

EeA, EeBand EeCdenote axial, axial-bending coupling and bending stiffness of the WFE respectively. The wavelet based physical space elemental stiffness matrix of the beam, kew, is

The kinetic energy of the functionally graded beam element, Λe, is defined as