Multiscale Wavelet Finite Element Analysis in Structural Dynamics

Abstract

Over the recent past, various numerical analysis techniques have been formulated and used to obtain approximate solutions for numerous engineering problems to aid predict the behaviour of systems accurately and efficiently. One such approach is the Wavelet Finite Element Method (WFEM) which involves combining the classical Finite Element Method (FEM) with wavelet analysis. The key desirable properties exhibited by some wavelet families, such as compact support, multiresolution analysis (MRA), smoothness, vanishing moments and the ‘two-scale’ relations, make the use of wavelets in WFEM advantageous, particularly in the analysis of problems with strong nonlinearities, singularities and material property variations present. The wavelet based finite elements (WFEs) of a rod and beam are formulated using the Daubechies and B-spline wavelet on the interval (BSWI) wavelet scaling functions as interpolating functions due to their desirable properties, thus making it possible to alter the local scale of the WFE without changing the initial model mesh. Specific benchmark cases are presented to exhibit and compare the performance of the WFEM with FEM in static, dynamic, eigenvalue and moving load transient response analysis for homogenous systems and functionally graded materials, where the material properties continuously vary spatially with respect to the constituent materials.

1. Introduction

In the analysis of complex structural problems, it is often challenging to formulate and apply exact closed-form solutions, as the realistic nature of such engineering systems exhibits varying complexities, high gradients and strong irregularities, e.g., suddenly varying loading conditions, contrasting material composition or geometric variations. Based on the existing mathematical tools available, such systems may require certain assumptions and generalisations to be implemented in order to simplify the model, which may lead to inability to correctly describe the properties and behaviour of the system under described conditions. However, the preferred approach is to find an approximate numerical solution, whilst retaining these complexities as accurately as possible, to better describe and predict the behaviour of such systems. This has given rise to numerical methods such as the classical Finite Element Method which employs polynomial interpolating functions to obtain approximate solutions for various engineering problems. Although this numerical analysis technique has grown in popularity, its use to tackle problems with regions of the solution domain where the gradient of the field variables are expected to vary suddenly or fast, bring on difficulties in the analysis of a complex system [1]. In order to improve on the accuracy and better represent the system’s behaviour, higher order polynomial interpolating functions or finer meshes may be employed and this in turn significantly increases the computational costs; which is undesirable. Moreover, the resolution of the elements can only be analysed to a specific scale once the orders of the governing polynomial functions have been selected. Subsequently, overcoming these challenges has been the driving force in the formulation of other numerical approximation techniques such as the Wavelet Finite Element Method [1, 2, 3, 4, 5, 6].

The initial development of wavelet analysis came from separate efforts that led to the foundation of modern wavelet theory. Grossman and Morlet [7] used wavelet analysis as a tool for signal analysis of seismic data and are credited with the introduction of the term and methodology of wavelets as it is known today. Ingrid Daubechies is recognised for her major breakthrough and contribution by constructing a family of orthonormal wavelet with compact support known as the Daubechies wavelets [8]. Wavelet analysis was used mainly by mathematicians as a decomposition tool for data functions and operators and its application has vastly grown in various disciplines at an exponential rate e.g., medicine [9], finance [10] and astronomy [11]. Likewise, the range of wavelet families and bases available for selection has also increased and this is credited to the properties of wavelets that allow it to be tailored to suite numerous avenues for design manipulation to meet the necessary and specific requirements for its application. The properties of different wavelet families vary, and therefore the decision on which family is the ‘most adequate’, is paramount to its application. Nevertheless, the more general aspects of wavelets formulations make it an important and convenient tool for mathematical manipulation allowing for the decomposition of a function into a set of coefficients that are dependent on scale and location. The ‘two-scale’ relation gives rise to one of the most key features of wavelet theory, multiresolution analysis (MRA), which allows for the convenient transformation of wavelet basis functions between different resolution scales [8]. Furthermore, the compact support property of wavelets ensures that the wavelet basis functions are finitely bound (non-zero over a finite range). The vanishing moments of wavelets allow the basic functions of wavelets to represent polynomials and other complex functions.

These desirable properties of wavelets have led to the use of wavelet basis functions as interpolating functions, in contrast to conventional polynomial functions as used in classical FEM, in the formulation of the wavelet based finite element method. For example, MRA permits for specific WFEs to be selected and analysed locally at finer scales without altering the initial system model, thus improving the accuracy of the solution, particularly in areas with high gradients or singularities present. Furthermore, rapid convergence of the method and compact support lead to a reduction in computational costs since fewer elements are required to achieve acceptable levels of accuracy [4, 5]. Due to the adaptability of wavelets, different wavelet families are being developed and customised for specific problems. However, it must be noted that when selecting a particular wavelet basis function for WFEM, key requirements, such as compatibility, completeness and convergence, must be satisfied and should allow for the easy implementation and treatment of boundary conditions.

The Daubechies wavelet based finite element was first introduced to solve a 1D and 2D second order Neumann problem via the formulation of a tensor product finite element [2]. The Daubechies wavelet Galerkin finite element was then used to analyse the bending of plates and beams [12] giving rise to the formulation of a wavelet based beam finite element [6] and two dimensional Daubechies wavelet plate finite element [13] for static analysis. The Daubechies wavelet base finite element stiffness matrices and load vectors were presented by Chen et al. at multiresolution scale j = 0 [14] and different multiresolution scales [4]. The Daubechies plate finite element was developed by Diaz et al. for the static analysis of plates based on Mindlin-Reissner plate theory [15], where shear deformation is taken into consideration through the thickness of the plate, and compared it with Kirchhoff plate theory formulations [16]. This wavelet family has also been used in the analysis of many other structural problems, including formulation of the Rayleigh-Euler and Rayleigh-Timoshenko beam elements [17], the wavelet based spectral finite element to study elastic wave propagation in 1-D connected waveguides [18] and also to investigate the thermal stress distribution along the vertical direction of the tank wall [19]. Overall, the wavelet family performed decently in providing accurate solutions for the various structural analysis problems tackled. However, the Daubechies wavelet lacks an explicit expression for the wavelet and scaling functions and possesses unusual smoothness characteristics, particularly for lower orders, making it challenging to evaluate the numerical integrals necessary for the formulation of the element matrices and load vectors. The evaluation of the connection coefficients is therefore necessary for the formulation of these element matrices and vectors.

In a bid to overcome the limitations presented by the Daubechies wavelet, further research has been carried out to identify other potential wavelet families that can be implemented in WFEM. Basic spline functions were initially used as interpolating functions for the free vibration analysis of frame structures [3]. Chui and Quak [20] constructed the semi-orthogonal B-spline Wavelet on the Interval, which has the desirable properties of multiresolution, compact support, explicit expressions, smoothness and symmetry. The BSWI was employed to construct the wavelet based C⁰ type plane elastomechanics element and Mindlin plate element [21] as well as truncated conical shell wavelet finite elements [22]. Xiang et al. [5] significantly contributed to the use of BSWI in WFEM by constructing the axial rod, beam (Timoshenko and Euler Bernoulli) and spatial bar WFEs with a multiresolution lifting scheme. Furthermore, this research was extended to the static and dynamic analysis of plates based on Kirchhoff plate theory using BSWI based wavelet finite elements [23, 24]. Xiang et al. [25] were able to illustrate that the shear-locking phenomenon of a rotating Rayleigh-Timoshenko shaft was significantly eliminated when the BSWI based WFEs were employed. Majority of the problems examined by this point were of static analysis and this led Musuva and Mares [26] to develop and implement the Daubechies and BSWI homogenous beam WFEs for the analysis of dynamic response and moving load problems. The vibration and dynamic response analysis was carried out for frame structures using the two wavelet families [27] and the WFEM was compared with an analytical wavelet approach using coiflets for the analysis of vehicle-bridge interaction for fast moving loads [28]. Furthermore, the Daubechies and BSWI wavelets were used to construct a functionally graded beam wavelet finite element under various moving load conditions [29, 30].

Other different wavelet families have been selected and employed in the formulation of the WFEM to solve a wide variety of structural analysis problems and research in this field is still ongoing. The trigonometric Hermite wavelet, which can be explicitly expressed, was used to construct beam [31] and thin plate WFEs [32] for static and free vibration analysis. The Hermite Cubic Spline Wavelet on the Interval (HCSWI), polynomial wavelets [33] and the second generation wavelets [34] are other wavelet based approaches that have been introduced and researched on. A more comprehensive synthesis and summary of wavelet based numerical methods for various engineering problems is presented in [35].

A generalised Wavelet based Finite Element Method framework is presented based on the BSWI and Daubechies wavelet families to derive rod and beam WFEs for homogenous and functionally graded materials for static and dynamic structural problems. A brief introduction of wavelet analysis is described in Section 2, with emphasis given to the Daubechies wavelets, BSWI, multiresolution and connection coefficients formulations. In Section 3, the wavelet based finite elements for a rod, Euler Bernoulli homogeneous beam and transversely varying functionally graded beam are presented. The evaluation of the element matrices and various load vectors, including the WFEM moving load formulation, are presented. A comparison on the performance of the Daubechies and BSWI WFEMs are highlighted via numerical examples for a variety of static and dynamic structural problems in Section 4 followed by conclusions.

2. Wavelet and multiresolution analysis

Wavelets are a class of basic functions that represent functions locally, both in space and time, and allow for the analysis of functions to be carried out at different resolutions (scales) [36]. The wavelet basis emanates from a set of wavelet coefficients associated with a particular location in time and different multiresolution scales. The scaling and wavelet functions stem from multiresolution analysis (MRA), which is a key and desirable property of wavelets, and refers to the simultaneous appearance of multiple scales in function decompositions in the Hilbert space $L^2\left(\mathbb{R}\right)$ using a sequence of closed subspaces V_j, which is represented mathematically as [36]:

$$\cdots V_{-2}\subset V_{-1}\subset V_0\subset V_1\subset V_2\subset \dots$$E1

Therefore in principle, in order for multiresolution to occur, the closed subspaces $V_j$ satisfy the following properties:

$$\overline{\bigcup _{j\in \mathbb{Z}}{V}_j}=L^2\left(\mathbb{R}\right)$$E2

$$\bigcap _{j\in \mathbb{Z}}{V}_j=\left\{0\right\}$$E3

$$\begin{array}{c}{f}_2\left(x\right)=f\left(2x\right)\forall x\\ f\in V_j\iff f_2\in V_{j+1}j\in \mathbb{Z}\end{array}$$E4

$$\begin{array}{c}{f}_n\left(x\right)=f\left(x-n\right)\\ f\in V_0\iff f_n\in V_{0}n\in \mathbb{Z}\end{array}$$E5

The orthogonal complement subspace $W_j$ of $V_j$ contains the additional ‘detail’ for subspace $V_{j+1}$ i.e., $V_{j+1}=V_0\oplus W_0\oplus W_1\oplus W_2\cdots \oplus W_j$. The union of the subspaces $V_j$ leads to the space $L^2\left(\mathbb{R}\right)$ from the condition in Eq. (2) [36]. The scaling $\varphi \left(x\right)\in L^2\left(\mathbb{R}\right)$ and wavelet $\psi \left(x\right)\in L^2\left(\mathbb{R}\right)$ functions correspond to the subspaces $V_j$ and $W_j$ respectively. The difference between current subspace V_j and subsequent subspace $V_{j+1}$ is represented by the wavelet space $W_j$ which becomes automatically orthogonal to all other $W_j\text{for}k<j$ due to the inclusion in and orthogonality to $V_j$. For the fundamental space $V_0$, the scaling function $\varphi \left(x\right)$ and its translates $\varphi \left(x-k\right)$ produce an orthonormal basis for V₀. The orthonormal basis for the next space $V_1$ is the rescaled function $\sqrt{2}\varphi \left(2x-k\right)$. Thus, the orthonormal basis of $V_j$ is defined as:

$${\varphi }_k^j\left(x\right)=2^{\frac{j}{2}}\varphi \left(2^{j}x-k\right)k\in \mathbb{Z}$$E6

Provided Eq. (6) and the above mentioned properties are satisfied, the wavelet orthonormal basis for subspace $W_j$ at scale j is

$${\psi }_k^j\left(x\right)=2^{\frac{j}{2}}\psi \left(2^{j}x-k\right)k\in \mathbb{Z}$$E7

The orthogonal subspaces $W_j$ result from the decomposition of $L^2\left(\mathbb{R}\right)$ and subsequently the functions within these subspaces inherit the scale and shift invariance properties from the scaling function subspaces $V_j$ and are orthonormal [8]. The projections of a function $f\in L^2\left(\mathbb{R}\right)$ at scale j in the subspaces $V_j$ and $W_j$, defined as $P_{j}f$ and $Q_{j}f$ respectively, are expressed as:

$$\begin{array}{c}{P}_{j}f=\sum _k{a}_k^j{\varphi }_k^j\left(x\right)\\ Q_{j}f=\sum _k{b}_k^j{\psi }_k^j\left(x\right)\end{array}$$E8

where $a_k^j$ and $b_k^j$ are coefficients in the subspaces $V_j$ and $W_j$ respectively. Thus, if all the conditions described above are met, then the scaling and wavelet functions satisfy [8]

$$\begin{array}{c}{\int }_{-\infty }^{\infty }\varphi \left(x\right)\mathit{dx}\ne 0\\ {\int }_{-\infty }^{\infty }\psi \left(x\right)\mathit{dx}=0\end{array}$$E9

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]:

$${\varphi }_L\left(x\right)=\sum _{k=0}^{L-1}{p}_L\left(k\right){\varphi }_L\left(2x-k\right)$$E10

$${\psi }_L\left(x\right)=\sum _{k=0}^{L-1}{q}_L\left(k\right){\varphi }_L\left(2x-k\right)$$E11

The scaling and wavelet functions have the supports $\left[0,L-1\right]$ and $\left[1-\frac{L}{2},\frac{L}{2}\right]$ respectively. The normalised wavelet function filter coefficients $q_L\left(k\right)$ and scaling function filter coefficients $p_L\left(k\right)$ have the relation $q_L\left(k\right)={\left(-1\right)}^k{p}_L\left(1-k\right)$. The multiresolution scaling and wavelet basis functions corresponding to the subspaces V_j and W_j are defined as:

$${\varphi }_{L,k}^j\left(x\right)=2^{\frac{j}{2}}{\varphi }_L\left(2^{j}x-k\right)$$E12

$${\psi }_{L,k}^j\left(x\right)=2^{\frac{j}{2}}{\psi }_L\left(2^{j}x-k\right)$$E13

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

$${\int }_{-\infty }^{\infty }{\varphi }_L\left(x\right)\mathit{dx}=1$$E14

$${\int }_{-\infty }^{\infty }{\varphi }_{L,k}^j\left(x\right){\varphi }_{L,l}^j\left(x\right)\mathit{dx}={\delta }_{k,l}$$E15

$${\int }_{-\infty }^{\infty }{\psi }_{L,k}^j\left(x\right){\psi }_{L,k}^j\left(x\right)\mathit{dx}={\delta }_{k,l}$$E16

$${\int }_{-\infty }^{\infty }{\varphi }_{L,k}^j\left(x\right){\psi }_{L,l}^j\left(x\right)\mathit{dx}=0$$E17

$${\int }_{-\infty }^{\infty }{x}^m{\psi }_L\left(x\right)\mathit{dx}=0m=0,1,\dots ,\frac{L}{2}-1$$E18

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 $\frac{L}{2}-1$ vanishing moments from property (18) and consequently the scaling functions at scale j can represent a polynomial of order $x^m$ where $0\le m\le \frac{L}{2}-1$, i.e., [37]

$$x^m=\sum _k{M}_k^{j,m}{\varphi }_{L,k}^j\left(x\right)$$E19

The coefficients $M_k^{j,m}$ denote the moments of the scaling function and it translates at V_j. The derivatives of the Daubechies wavelet scaling functions are evaluated by differentiating the refinement Eq. (10) m times, and are obtained as [12]:

$${\varphi }_L^{\left(m\right)}\left(x\right)=2^m\sum _{k=0}^{L-1}p\left(k\right){\varphi }_L^{\left(m\right)}\left(2x-k\right)$$E20

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

$$\sum _{k=-\infty }^{k=\infty }{k}^m{\varphi }^{\left(m\right)}\left(x-k\right)=m!$$E21

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,b}\mathrm{\Gamma }}_{k,l}^{j,d_1,d_2}$ and multiscale connection coefficient ${\mathrm{{\rm Y}}}_k^{j,m}$. We define the two-term connection coefficient [30]

$${{}_{a,b}\mathrm{\Gamma }}_{k,l}^{j,d_1,d_2}=2^j{\int }_{-\infty }^{\infty }{\mathcal{X}}_{\left[0,1\right]}\left(\xi \right){\varphi }_a^{\left(d_1\right)}\left(2^j\xi -k\right){\varphi }_b^{\left(d_2\right)}\left(2^j\xi -l\right)\mathit{d\xi }$$E22

where a and b are the orders of the scaling function at multiresolution j, while the values $d_1$ and $d_2$ denote the order of the derivative of the scaling functions. ${\mathcal{X}}_{\left[0,1\right]}\left(x\right)=\left\{\begin{array}{l}10\le x\le 1\hfill \\ 0\mathit{\text{otherwise}}\end{array}\right.$ is 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),

$${\varphi }_L\left(2^j\xi -k\right)=\sum _{r}p\left(r\right){\varphi }_L\left(2^{j+1}\xi -2k-r\right)$$E23

Differentiating Eq. (23) m times

$$2^{\mathit{jm}}{\varphi }_L^{\left(m\right)}\left(2^j\xi -k\right)=2^{\left(j+1\right)m}\sum _{r}p\left(r\right){\varphi }_L^{\left(m\right)}\left(2^{j+1}\xi -2k-r\right)$$E24

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:

$${{}_{a,b}\mathrm{\Gamma }}_{k,l}^{j,d_1,d_2}=2^{d_1+d_2-1}\sum _{r,s}\left[p_a\left(r-2k\right)p_b\left(s-2l\right)+p_a\left(r-2k+2^j\right)p_b\left(s-2l+2^j\right)\right]{\mathrm{\Gamma }}_{r,s}^{j,d_1,d_2}$$E25

where $2-a\le k,r\le 2^j-1$ and $2-b\le l,s\le 2^j-1$. Eq. (25) can be expressed in matrix form as:

$$\begin{array}{l}{}_{\left(\left(a+2^j-2\right)\left(b+2^j-2\right)\times 1\right)}\left\{{}_{a,b}{\mathbf{\Gamma }}^j\right\}=\\ 2^{d_1+d_2-1}{}_{\left(\left(a+2^j-2\right)\left(b+2^j-2\right)\times \left(a+2^j-2\right)\left(b+2^j-2\right)\right)}\left[{}_{a,b}\mathit{P}\right]{}_{\left(\left(a+2^j-2\right)\left(b+2^j-2\right)\times 1\right)}\left\{{}_{a,b}{\mathbf{\Gamma }}^j\right\}\end{array}$$E26

where the square matrix $\left[{}_{\mathit{a},\mathit{b}}\mathit{P}\right]$ contains the filter coefficients as expressed in Eq. (25) and $\left\{{}_{a,b}{\mathrm{\Gamma }}^j\right\}$ contains 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)

$${\xi }^m=2^{\frac{j}{2}}\sum _k{{}_{L}M}_k^{j,m}{\varphi }_L\left(2^j\xi -k\right)$$E27

Defining the second form of the connection coefficient

$${\mathrm{{\rm Y}}}_k^{j,m}=2^{\frac{j}{2}}{\int }_0^1{x}^m{\varphi }_L\left(2^j\xi -k\right)\mathit{d\xi }=2^{\frac{j}{2}}{\int }_{-\infty }^{\infty }{\mathcal{X}}_{\left[0,1\right]}\left(\xi \right){\xi }^m{\varphi }_L\left(2^j\xi -k\right)\mathit{d\xi }$$E28

Substituting Eq. (27) into (28)

$${\mathrm{{\rm Y}}}_k^{j,m}=2^j\sum _l{M}_l^{j,m}{\int }_{-\infty }^{\infty }{\mathcal{X}}_{\left[0,1\right]}\left(x\right){\varphi }_L\left(2^{j}x-l\right){\varphi }_L\left(2^{j}x-k\right)\mathit{dx}$$E29

However,

$${{}_{L,L}\mathrm{\Gamma }}_{k,l}^{j,0,0}=2^j{\int }_{-\infty }^{\infty }{\mathcal{X}}_{\left[0,1\right]}\left(x\right){\varphi }_L\left(2^{j}x-l\right){\varphi }_L\left(2^{j}x-k\right)\mathit{dx}$$E30

Thus

$${\mathrm{{\rm Y}}}_k^{j,m}=\sum _l{M}_l^{j,m}{{}_{L,L}\mathrm{\Gamma }}_{k,l}^{j,0,0}$$E31

where ${{}_{L,L}\mathrm{\Gamma }}_{k,l}^{j,0,0}$ are the two-term connection coefficients with $a=b=L$ and $d_1=d_1=0$ and $M_l^{j,m}$ are 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 $V_j$ of order m and scale $j>0$ are expressed as [20]

$$B_{m,k}^j\left(x\right)=\left(t_{k+m}^j-t_k^j\right){\left[t_k^j,\dots ,t_{k+m}^j\right]}_f{\left(t-x\right)}_{+}^{m-1}$$E32

with the knot sequence

$$\begin{array}{c}{\left\{t_k^j\right\}}_{k=-m+1}^{2^j+m-1}\\ t_k^j\le t_{k+1}^j\end{array}$$E33

${\left[t_k^j,t_{k+1}^j,\dots ,t_{k+m}^j\right]}_t$, is the m^th divided difference of the truncated power function ${\left(t-x\right)}_{+}^{m-1}$ with respect to variable t. The general B-splines take the form

$$\begin{array}{cc}& B_{m,k}^j\left(x\right)=\frac{x-t_k^j}{t_{k+m-1}^j-t_k^j}{B}_{m-1,k}^j\left(x\right)+\frac{t_{k+m}^j-x}{t_{k+m}^j-t_{k+1}^j}{B}_{m-1,k+1}^j\left(x\right)\\ & B_{1,k}^j\left(x\right)=& \left\{\begin{array}{c}\begin{array}{cc}1& k\le x\le k+1\end{array}\\ \begin{array}{cc}0& \mathit{\text{otherwise}}\end{array}\end{array}\right)\end{array}$$E34

and have support ${}_{\mathit{supp}}{B}_{m,k}^j\left(x\right)=\left[t_k^j,t_{k+m}^j\right]$. 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], $t_k^j$ is given as [38]:

$$t_k^j=\left\{\begin{array}{cc}& \begin{array}{cc}0& -m+1\le k<1\end{array}\\ & \begin{array}{cc}{2}^{-j}k& 1\le k<2^j\end{array}\\ & \begin{array}{cc}1& 2^j\le k\le 2^j+m-1\end{array}\end{array}\right.$$E35

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 $j_0$:

$$2^{j_0}\ge 2m-1$$E36

The basis $B_{m,k}^j\left(x\right)$ from the inner knots corresponds to the m^th cardinal B-splines, $N_m\left(x\right)$, at multiresolution j [38]:

$$N_m\left(x\right)=m\left[0,1,\dots ,m\right]{\left(t-x\right)}_{+}^{m-1}$$E37

$${\varphi }_{m,k}^j\left(x\right)=B_{m,k}^j\left(x\right)=N_m\left(2^{j}x-k\right)0\le k<2^j-m+1$$E38

where ${\varphi }_{m,k}^j\left(x\right)$ is the BSWI scaling function which can be differentiated m times. The corresponding B-wavelet with support ${}_{\mathit{supp}}{\psi }_{m,k}^j\left(x\right)=\left[\frac{k}{2^j},\frac{k+2m-1}{2^j}\right]$ is expressed as:

$${\psi }_{m,k}^j\left(x\right)=\frac{1}{2^{m-1}}\sum _{l=0}^{2m-2}{\left(-1\right)}^l{N}_{2m}\left(l+1\right)B_{2m,2i+l}^{j+1,\left(m\right)}\left(x\right)$$E39

$B_{2m,k}^{j+1,\left(m\right)}\left(x\right)$ is the m^th derivative for the B-spline of order 2m and scale $j+1$ and can be evaluated explicitly from Eq. (34). Given that the requirement $j>j_0$ ensures at least one inner B-wavelet is present, the scaling and wavelet function of the BSWI are obtained as [39]:

$${\varphi }_{m,k}^j\left(x\right)=\left\{\begin{array}{cc}\begin{array}{cc}& B_{m,k}^{j_0}\left(2^{j-j_0}x\right)\\ & B_{m,0}^{j_0}\left(2^{j-j_0}x-2^{-j_0}k\right)\\ & B_{m,2^j-k-m}^{j_0}\left(1-2^{j-j_0}x\right)\end{array}& \begin{array}{c}-m+1\le k\le -1\\ 0\le i\le 2^j-m\\ 2^j\le i\le 2^j+m-1\end{array}\end{array}\right.$$E40

$${\psi }_{m,k}^j\left(x\right)=\left\{\begin{array}{ll}\begin{array}{cc}& {\psi }_{m,k}^{j_0}\left(2^{j-j_0}x\right)\\ & {\psi }_{m,0}^{j_0}\left(2^{j-j_0}x-2^{-j_0}k\right)\\ & {\psi }_{m,2^j-k-2m+1}^{j_0}\left(1-2^{j-j_0}x\right)\end{array}\hfill & \begin{array}{c}-m+1\le k\le -1\\ 0\le i\le 2^j-m\\ 2^j\le i\le 2^j+m-1\end{array}\end{array}\right.$$E41

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

3. The wavelet finite element method

3.1. Axial rod wavelet finite element

Assume each WFE is divided into equal segments, n_s, connected by $r=n_s+1$ elemental nodes, as shown in Figure 1, with axial deformation $u_i$. The total number of degrees of freedom (DOFs) within each WFE is denoted by $n=r$ for $n,r\in \mathbb{N}$. Vector $\left\{{\mathit{u}}_{\mathit{e}}\right\}={\left\{u_1{u}_2\cdots u_{r-1}{u}_r\right\}}^T$ contains all the axial DOFs in physical space, as illustrated in Figure 2(a), where $u_i=u\left(x_i\right)$ represents the elemental node axial deformation DOF at node i corresponding to coordinate position x_i. The nodal natural coordinates is ${\xi }_i=\frac{x_i-x_1}{L_e}$ (0 ≤ ξ_i ≤ 1, 1 ≤ i ≤ r). The Daubechies and BSWI scaling functions ${\varphi }_{z,k}^j\left(x\right)$ are used as the interpolating functions and for a family of order z at multiresolution scale j, the axial deformation

$$u\left(\xi \right)=\sum _{k=h}^{2^j-1}{a}_{z,k}^j{\varphi }_{z,k}^j\left(\xi \right)$$E42

contains the unknown wavelet coefficients $a_{z,k}^j$. This gives rise to the vector $\left\{{\mathit{u}}_{\mathit{e}}\right\}$ containing the axial deformations at all elemental nodes in physical space.

$${}_{\left(n\times 1\right)}\left\{{\mathit{u}}_{\mathit{e}}\right\}={}_{\left(n\times n\right)}\left[{\mathit{R}}_{\mathit{r}}^{\mathit{w}}\right]{}_{\left(n\times 1\right)}\left\{{\mathit{a}}_{\mathit{e}}\right\}$$E43

Figure 1.

Wavelet finite element layout.

Figure 2.

(a) Axial rod and (b) Euler Bernoulli beam wavelet finite element layout.

The matrix $\left[{\mathit{R}}_{\mathit{r}}^{\mathit{w}}\right]={\left[\begin{array}{ccc}\left\{{\mathbf{\Phi }}_{\mathit{z}}^{\mathit{j}}\left({\xi }_1\right)\right\}& \left\{{\mathbf{\Phi }}_{\mathit{z}}^{\mathit{j}}\left({\xi }_2\right)\right\}& \begin{array}{ccc}\cdots & \left\{{\mathbf{\Phi }}_{\mathit{z}}^{\mathit{j}}\left({\xi }_{r-1}\right)\right\}& \left\{{\mathbf{\Phi }}_{\mathit{z}}^{\mathit{j}}\left({\xi }_r\right)\right\}\end{array}\end{array}\right]}^{\mathrm{T}}$ contains the scaling function vectors $\left\{{\mathbf{\Phi }}_{\mathit{z}}^{\mathit{j}}\left({\xi }_i\right)\right\}$ approximating the axial deformation at the corresponding elemental nodes and $\left\{{\mathit{a}}_{\mathit{e}}\right\}={\left[\begin{array}{ccc}{a}_{z,h}^j& a_{z,h+1}^j& \begin{array}{ccc}\cdots & a_{z,2^j-2}^j& a_{z,2^j-1}^j\end{array}\end{array}\right]}^{\mathrm{T}}$. The axial deformation at any point along the rod element can be generalised as:

$$u\left(\xi \right)={}_{\left(1\times n\right)}\left\{{\mathbf{\Phi }}_{\mathit{z}}^{\mathit{j}}\left(\mathit{\xi }\right)\right\}{}_{\left(n\times n\right)}\left[{\mathit{T}}_{\mathit{r}}^{\mathit{w}}\right]{}_{\left(n\times 1\right)}\left\{{\mathit{u}}_{\mathit{e}}\right\}$$E44

The matrix $\left[{\mathit{T}}_{\mathit{r}}^{\mathit{w}}\right]={\left[{\mathit{R}}_{\mathit{r}}^{\mathit{w}}\right]}^{-1}$ is 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 $\left\{{\mathit{N}}_{\mathit{r},\mathit{e}}\left(\xi \right)\right\}=\left\{{\mathbf{\Phi }}_{\mathit{z}}^{\mathit{j}}\left(\xi \right)\right\}\left[{\mathit{T}}_{\mathit{r}}^{\mathit{w}}\right]$ within each element.

Suppose the axial rod is subjected to nodal point loads $f_{\mathit{xi}}$ and distributed loading$f_d\left(x\right)$, then the potential energy within the axial rod ${\Pi }^a$ can be generally expressed as [40]:

$${\Pi }^a={\int }_0^l\frac{\mathit{EA}}{2}{\left(\frac{\mathit{du}\left(x\right)}{\mathit{dx}}\right)}^2\mathit{dx}-\sum _{i}u\left(x_i\right)f_{\mathit{xi}}-{\int }_0^l{f}_d\left(x\right)u\left(x\right)\mathit{dx}$$E45

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 $U_e^a$ within each WFE of length $L_e$ is expressed in natural coordinates as:

$$U_e^a=\frac{1}{2}\frac{\mathit{EA}}{L_e}{\left\{{\mathit{u}}_{\mathit{e}}\right\}}^T{\int }_0^1{\left[{\mathit{T}}_{\mathit{r}}^{\mathit{w}}\right]}^T{\left\{\frac{d{\mathbf{\Phi }}_{\mathit{z}}^{\mathit{j}}\left(\xi \right)}{\mathit{d\xi }}\right\}}^T\left\{\frac{d{\mathbf{\Phi }}_{\mathit{z}}^{\mathit{j}}\left(\xi \right)}{\mathit{d\xi }}\right\}\left[{\mathit{T}}_{\mathit{r}}^{\mathit{w}}\right]\mathit{d\xi }\left\{{\mathit{u}}_{\mathit{e}}\right\}$$E46

The stiffness matrix of the rod element in wavelet space, $\left[{\mathit{k}}_{\mathit{r},\mathit{e}}^{\mathit{w}}\right]$ is computed using the first derivative of the scaling functions and is symmetric.

$${}_{\left(n\times n\right)}\left[{\mathit{k}}_{\mathit{r},\mathit{e}}^{\mathit{w}}\right]={\int }_0^1{\left\{{\mathbf{\Phi }}_{\mathit{z}}^{\prime \mathit{j}}\left(\xi \right)\right\}}^T\left\{{\mathbf{\Phi }}_{\mathit{z}}^{\prime \mathit{j}}\left(\xi \right)\right\}\mathit{d\xi }$$E47

In order for one to obtain the stiffness matrix in physical space, the element properties and transformation matrix $\left[{\mathit{T}}_{\mathit{r}}^{\mathit{w}}\right]$ are applied to the wavelet space stiffness matrix in Eq. (47).

$${}_{\left(n\times n\right)}\left[{\mathit{k}}_{\mathit{r},\mathit{e}}^{\mathit{p}}\right]=\frac{\mathit{EA}}{L_e}{{}_{\left(n\times n\right)}\left[{\mathit{T}}_{\mathit{r}}^{\mathit{w}}\right]}^T{}_{\left(n\times n\right)}\left[{\mathit{k}}_{\mathit{r},\mathit{e}}^{\mathit{w}}\right]{}_{\left(n\times n\right)}\left[{\mathit{T}}_{\mathit{r}}^{\mathit{w}}\right]$$E48

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

$${}_{\left(n\times 1\right)}\left\{{\mathit{f}}_{\mathit{r},\mathit{e}}^{\mathit{n},\mathit{p}}\right\}=\sum _i{\left[{\mathit{T}}_{\mathit{r}}^{\mathit{w}}\right]}^T{\left\{{\mathbf{\Phi }}_{\mathit{z}}^{\mathit{j}}\left({\xi }_i\right)\right\}}^T{f}_{\mathit{xi}}$$E49

and the equivalent nodal load vector for the distributed load $f_d\left(x\right)$ in physical space is

$${}_{\left(n\times 1\right)}\left\{{\mathit{f}}_{\mathit{r},\mathit{e}}^{\mathit{d},\mathit{p}}\right\}=L_e{\int }_0^1{f}_d\left(\xi \right){\left[{\mathit{T}}_{\mathit{r}}^{\mathit{w}}\right]}^T{\left\{{\mathbf{\Phi }}_{\mathit{z}}^{\mathit{j}}\left(\xi \right)\right\}}^T\mathit{d\xi }$$E50

When applying the Daubechies wavelet family, the WFE has a total of $n=2^j+L-2$ DOFs. The wavelet space stiffness matrix is evaluated from the multiscale two-term connection coefficients ${{}_{a,b}\mathrm{\Gamma }}_{k,l}^{j,d_1,d_2}$ $a=b=L$ and $d_1=d_1=1$ and is given as:

$${}_{\left(\left(2^j+L-2\right)x\left(2^j+L-2\right)\right)}{}^{\mathrm{ D}}\left[{\mathit{k}}_{\mathit{r},\mathit{e}}^{\mathit{w}}\right]=2^{2j}\left[{\mathbf{\Gamma }}^{j,1,1}\right]$$E51

where $\left(2^{2j}\right)$ is the normalising factor and the matrix $\left[{\mathbf{\Gamma }}^{j,1,1}\right]$ has the entries ${{}_{L,L}\mathrm{\Gamma }}_{k,l}^{j,1,1}$ for the limits $2-L\le k,l\le 2^j-1$. Similarly, the distributed forces acting on the element require the form ${{\rm Y}}_k^{j,m}$ for limits $2-L\le k,l\le 2^j-1$ of connection coefficients and the value of m depends on the order of the function $f_d\left(x\right)$ of the forces. In the case of the BSWI formulations, the total DOFs is $n=2^j+m-1$ and the condition $j\ge j_0$ must be satisfied. Therefore, the wavelet space stiffness matrices of the BSWI axial rod are computed as:

$${}_{\left(\left(2^j+m-1\right)\times \left(2^j+m-1\right)\right)}{}^{\mathit{ BS}}\left[{\mathit{k}}_{\mathit{r},\mathit{e}}^{\mathit{w}}\right]={\int }_0^1{\left\{{\mathbf{\Phi }}_{\mathit{m}}^{\prime \mathit{j}}\left(\xi \right)\right\}}^T\left\{{\mathbf{\Phi }}_{\mathit{m}}^{\prime \mathit{j}}\left(\xi \right)\right\}\mathit{d\xi }$$E52

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 L_e, is divided into n_s equally spaced elemental segments connected by r elemental nodes at coordinate values $x_i\in \left[x_1,x_r\right]$ and $i\in \mathbb{N}$ as illustrated in Figure 1. The WFE has the transverse displacement v and rotation θ taken into account, with corresponding transverse forces $f_y$ and moments ${\displaystyle \acute{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-2$ displacement DOFs and 2 rotation DOFs in total for each WFE and consequently $r=n-2$ elemental nodes and $n_s=n-3$ elemental segments. Let the vector {${\mathit{v}}_{\mathit{e}}\}={\left\{\begin{array}{ccccccccc}{v}_1& {\theta }_1& v_2& v_3& \cdots & v_{r-2}& v_{r-1}& v_r& {\theta }_r\end{array}\right\}}^T$ denote all the physical DOFs within the beam element. The displacement and rotation DOFs corresponding to coordinate position $x_i\in \left[x_1,x_r\right]i\in \mathbb{N}\text{and}\left(1\le i\le r\right)$ in local coordinates are denoted as $v_i=v\left(x_i\right)$ and ${\theta }_i=\theta \left(x_i\right)$. The nodal natural coordinate ${\xi }_i=\frac{x_i-x_1}{L_e}$ (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 ${\varphi }_{z,k}^j\left(x\right)$ of order z at multiresolution scale j as interpolating functions.

$$\begin{array}{c}v\left(\xi \right)=\sum _{k=h}^{2^j-1}{b}_{z,k}^j{\varphi }_{z,k}^j\left(\xi \right)\\ \theta \left(\xi \right)=\frac{\partial v\left(\xi \right)}{\partial x}=\frac{1}{L_e}\sum _{k=h}^{2^j-1}{b}_{z,k}^j\frac{\partial {\varphi }_{z,k}^j\left(\xi \right)}{\partial \xi }\end{array}$$E53

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

$${}_{\left(n\times 1\right)}\left\{{\mathit{v}}_{\mathit{e}}\right\}={}_{\left(n\times n\right)}\left[{\mathit{R}}_{\mathit{b}}^{\mathit{w}}\right]{}_{\left(n\times 1\right)}\left\{{\mathit{b}}_{\mathit{e}}\right\}$$E54

$\left[{\mathit{R}}_{\mathit{b}}^{\mathit{w}}\right]={\left[\begin{array}{ccc}\left\{{\mathbf{\Phi }}_{\mathit{z}}^{\mathit{j}}\left({\xi }_1\right)\right\}& \frac{1}{L_e}\left\{{\mathbf{\Phi }}_{\mathit{z}}^{\prime \mathit{j}}\left({\xi }_1\right)\right\}& \begin{array}{ccc}\left\{{\mathbf{\Phi }}_{\mathit{z}}^{\mathit{j}}\left({\xi }_2\right)\right\}& \cdots & \begin{array}{ccc}\left\{{\mathbf{\Phi }}_{\mathit{z}}^{\mathit{j}}\left({\xi }_{r-1}\right)\right\}& \left\{{\mathbf{\Phi }}_{\mathit{z}}^{\mathit{j}}\left({\xi }_r\right)\right\}& \frac{1}{L_e}\left\{{\mathbf{\Phi }}_{\mathit{z}}^{\prime \mathit{j}}\left({\xi }_r\right)\right\}\end{array}\end{array}\end{array}\right]}^{\mathrm{T}}$ and vector $\left\{{\mathit{b}}_{\mathit{e}}\right\}$ contains the unknown wavelet coefficients $b_{z,k}^j$ representing the beam wavelet space DOFs.

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

$$v\left(\xi \right)={}_{\left(1\times n\right)}\left\{{\mathbf{\Phi }}_{\mathit{z}}^{\mathit{j}}\left(\xi \right)\right\}{}_{\left(n\times n\right)}\left[{\mathit{T}}_{\mathit{b}}^{\mathit{w}}\right]{}_{\left(n\times 1\right)}\left\{{\mathit{v}}_{\mathit{e}}\right\}\theta \left(\xi \right)=\frac{1}{L_e}{}_{\left(1\times n\right)}\left\{{\mathbf{\Phi }}_{\mathit{z}}^{\prime \mathit{j}}\left(\xi \right)\right\}{}_{\left(n\times n\right)}\left[{\mathit{T}}_{\mathit{b}}^{\mathit{w}}\right]{}_{\left(n\times 1\right)}\left\{{\mathit{v}}_{\mathit{e}}\right\}$$E55

where $\left[{\mathit{T}}_{\mathit{b}}^{\mathit{w}}\right]={\left[{\mathit{R}}_{\mathit{b}}^{\mathit{w}}\right]}^{-1}$ is the beam wavelet transformation matrix which is used to obtain the wavelet based shape functions for the beam $\left\{{\mathit{N}}_{\mathit{b},\mathit{e}}\left(\xi \right)\right\}=\left\{{\mathbf{\Phi }}_{\mathit{z}}^{\mathit{j}}\left(\xi \right)\right\}\left[{\mathit{T}}_{\mathit{b}}^{\mathit{w}}\right]$. The potential energy ${\Pi }^b$ within a Euler Bernoulli beam subjected to concentrated forces $f_{\mathit{yi}}$, distributed force $f_d\left(x\right)$ and bending moments ${{\displaystyle \acute{m}}}_i$ can be generally expressed as [40]:

$${\Pi }^b={\int }_0^l\frac{\mathit{EI}}{2}{\left(\frac{d^{2}v}{d{x}^2}\right)}^2\mathit{dx}-\sum _i{f}_{\mathit{yi}}v\left(x_i\right)-{\int }_0^l{f}_d\left(x\right)\mathit{vdx}-\sum _k{\acute{m}}_k\frac{\mathit{dv}\left(x_k\right)}{\mathit{dx}}$$E56

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

$$U_e^b=\frac{1}{2}\frac{EI}{{L_e}^3}{\left\{{\mathit{v}}_{\mathit{e}}\right\}}^T{\int }_0^1{\left[{\mathit{T}}_{\mathit{b}}^{\mathit{w}}\right]}^T{\left\{\frac{d^2{\mathbf{\Phi }}_{\mathit{z}}^{\mathit{j}}\left(\xi \right)}{d{\xi }^2}\right\}}^T\left\{\frac{d^2{\mathbf{\Phi }}_{\mathit{z}}^{\mathit{j}}\left(\xi \right)}{d{\xi }^2}\right\}\left[{\mathit{T}}_{\mathit{b}}^{\mathit{w}}\right]\mathit{d\xi }\left\{{\mathit{v}}_{\mathit{e}}\right\}$$E57

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

$${}_{\left(n\times n\right)}\left[{\mathit{k}}_{\mathit{b},\mathit{e}}^{\mathit{w}}\right]={\int }_0^1{\left\{{\mathbf{\Phi }}_{\mathit{z}}^{″\mathit{j}}\left(\xi \right)\right\}}^T\left\{{\mathbf{\Phi }}_{\mathit{z}}^{″\mathit{j}}\left(\xi \right)\right\}\mathit{d\xi }$$E58

The vector $\left\{{{\mathbf{\Phi }}^{\prime \prime }}_{\mathit{z}}^{\mathit{j}}\left(\xi \right)\right\}=\left\{\begin{array}{ccccc}{\varphi }_{z,h}^{″j}\left(\xi \right)& {\varphi }_{z,h+1}^{″j}\left(\xi \right)& \cdots & {\varphi }_{z,2^j-2}^{″j}\left(\xi \right)& {\varphi }_{z,2^j-1}^{″j}\left(\xi \right)\end{array}\right\}$ 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 $\left[{\mathit{T}}_{\mathit{b}}^{\mathit{w}}\right]$.

$${}_{\left(n\times n\right)}\left[{\mathit{k}}_{\mathit{b},\mathit{e}}^{\mathit{p}}\right]=\frac{EI}{{L_e}^3}{{}_{\left(n\times n\right)}\left[{\mathit{T}}_{\mathit{b}}^{\mathit{w}}\right]}^T{}_{\left(n\times n\right)}\left[{\mathit{k}}_{\mathit{b},\mathit{e}}^{\mathit{w}}\right]{}_{\left(n\times n\right)}\left[{\mathit{T}}_{\mathit{b}}^{\mathit{w}}\right]$$E59

The transverse kinetic energy of the beam element is expressed as

$${\Lambda }_e^b=\frac{1}{2}\mathit{\rho A}{L}_e{\int }_0^1\dot{v}{\left(\xi \right)}^T\dot{v}\left(\xi \right)\mathit{d\xi }$$E60

where $\dot{v}\left(\xi \right)=\frac{\partial v\left(\xi \right)}{\partial t}$, $\rho$ 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

$${\Lambda }_e^b={\left\{{\dot{\mathit{v}}}_{\mathit{e}}\right\}}^T\frac{1}{2}\mathit{\rho A}{L}_e{\int }_0^1{\left[{\mathit{T}}_{\mathit{b}}^{\mathit{w}}\right]}^T{\left\{{\mathbf{\Phi }}_{\mathit{z}}^{\mathit{j}}\left(\xi \right)\right\}}^T\left\{{\mathbf{\Phi }}_{\mathit{z}}^{\mathit{j}}\left(\xi \right)\right\}\left[{\mathit{T}}_{\mathit{b}}^{\mathit{w}}\right]\mathit{d\xi }\left\{{\dot{\mathit{v}}}_{\mathit{e}}\right\}$$E61

The mass matrix in physical space of the Euler Bernoulli beam element, $\left[{\mathit{m}}_{\mathit{b},\mathit{e}}^{\mathit{p}}\right]$, can be evaluated as:

$$\left[{\mathit{m}}_{\mathit{b},\mathit{e}}^{\mathit{p}}\right]=\mathit{\rho A}{L}_e{\left[{\mathit{T}}_{\mathit{b}}^{\mathit{w}}\right]}^T{\int }_0^1{\left\{{\mathbf{\Phi }}_{\mathit{z}}^{\mathit{j}}\left(\xi \right)\right\}}^T\left\{{\mathbf{\Phi }}_{\mathit{z}}^{\mathit{j}}\left(\xi \right)\right\}\mathit{d\xi }\left[{\mathit{T}}_{\mathit{b}}^{\mathit{w}}\right]$$E62

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

$$\begin{array}{c}{}_{\left(n\times 1\right)}{}^{}\left\{f_{b,e}^{n,p}\right\}={\displaystyle \sum }_{i=1}^r{}_{\left(n\times n\right)}{}^{}\left[T_b^w\right]{}^T{}_{\left(n\times 1\right)}{}^{}\left\{{\Phi }_z^j\left({\xi }_i\right)\right\}{}^T{f}_{yi}\\ {}_{\left(n\times 1\right)}{}^{}\left\{f_{b,e}^{m,p}\right\}={\displaystyle \sum }_k{}_{\left(n\times n\right)}{}^{}\left[T_b^w\right]{}^T{}_{\left(n\times 1\right)}{}^{}\left\{{\Phi }_z^{\prime j}\left({\xi }_k\right)\right\}{}^T{\acute{m}}_k\\ {}_{\left(n\times 1\right)}{}^{}\left\{f_{b,e}^{d,p}\right\}=L_e\underset{0}{\overset{1}{{\displaystyle \int }}}{f}_d\left(\xi \right){}_{\left(n\times n\right)}{}^{}\left[T_b^w\right]{}^T{\left\{{\Phi }_z^j\left(\xi \right)\right\}}^{T}d\xi \end{array}$$E63

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 $\left(x,t\right)=\mathit{P\delta }\left(x-x_0\right)$ [41]. $\delta \left(x\right)$ is the Dirac Delta function and $x_0$ is the distance travelled by the moving load at time t. The potential work of the load at this instant at position ${\xi }_0=\frac{x_0}{L_e}$ in natural coordinates is [1, 30]:

$${\mathrm{\Omega }}_e^b\left({\xi }_0\right)={\int }_0^1\mathit{P\delta }\left(\xi -{\xi }_0\right)v\left(\xi \right)\mathit{d\xi }=P{\left\{{\mathit{v}}_{\mathit{e}}\right\}}^T{\left[{\mathit{T}}_{\mathit{b}}^{\mathit{w}}\right]}^T{\left\{{\mathbf{\Phi }}_{\mathit{z}}^{\mathit{j}}\left({\xi }_0\right)\right\}}^T$$E64

Figure 3.

Layout of a beam WFE subjected to a moving point load.

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

$$\left\{{\mathit{f}}_{\mathit{b},\mathit{e}}^{\mathit{p},\mathit{p}}\left(t\right)\right\}=P{\left[{\mathit{T}}_{\mathit{b}}^{\mathit{w}}\right]}^T {}^t{\left\{{\mathbf{\Phi }}_{\mathit{z}}^{\mathit{j}}\left({\xi }_0\right)\right\}}^T$$E65

Assuming the moving load transverses to a new position ${\xi }_0$ within 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=2^j+L-2$ and for this specific layout, the total number of elemental nodes is $r=2^j+L-4$ and corresponding elemental segments $n_s=2^j+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:

$${}_{\left(\left(2^j+L-2\right)\times \left(2^j+L-2\right)\right)}{}^{\mathrm{ D}}\left[{\mathit{k}}_{\mathit{b},\mathit{e}}^{\mathit{w}}\right]=2^{4j}\left[{\mathbf{\Gamma }}^{j,2,2}\right]$$E66

$${}_{\left(\left(2^j+L-2\right)\times \left(2^j+L-2\right)\right)}{}^{\mathrm{ D}}\left[{\mathit{m}}_{\mathit{b},\mathit{e}}^{\mathit{w}}\right]=\left[{\mathbf{\Gamma }}^{j,0,0}\right]$$E67

Correspondingly, the connection coefficients of the form ${\mathrm{{\rm Y}}}_k^{j,m}$ for $2-L\le k\le 2^j-1$ are used to evaluated the distributed loads and the value of m is based on the load function $f_d\left(x\right)$. For the BSWI family of order m and at scale j, there are $n=2^j+m-1$ total DOFs, $r=2^j+m-3$ elemental nodes and $n_s=2^j+m-4$ elemental segments within the each WFE for this layout. The stiffness and mass matrices in wavelet space can be evaluated directly and are obtained as:

$${}_{\left(\left(2^j+m-1\right)\times \left(2^j+m-1\right)\right)}{}^{\mathit{ BS}}\left[{\mathit{k}}_{\mathit{b},\mathit{e}}^{\mathit{w}}\right]={\int }_0^1{\left\{{\mathbf{\Phi }}_{\mathit{m}}^{″\mathit{j}}\left(\mathit{\xi }\right)\right\}}^T\left\{{\mathbf{\Phi }}_{\mathit{m}}^{″\mathit{j}}\left(\mathit{\xi }\right)\right\}\mathit{d\xi }$$E68

$${}_{\left(\left(2^j+m-1\right)\times \left(2^j+m-1\right)\right)}{}^{\mathit{ BS}}\left[{\mathit{m}}_{\mathit{b},\mathit{e}}^{\mathit{w}}\right]={\int }_0^1{\left\{{\mathbf{\Phi }}_{\mathit{m}}^{\mathit{j}}\left(\mathit{\xi }\right)\right\}}^T\left\{{\mathbf{\Phi }}_{\mathit{m}}^{\mathit{j}}\left(\mathit{\xi }\right)\right\}\mathit{d\xi }$$E69

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]

$$P\left(y\right)=P_{\mathit{lo}}\left({\left[P_{\mathit{\text{ratio}}}-1\right]\left(\frac{y}{h}+\frac{1}{2}\right)}^n+1\right)$$E70

Figure 4.

(a) Cross-section of transversely varying functionally graded beam. (b) Effective Young’s modulus variation of steel-alumina functionally graded beam for different n. (c) Functionally graded beam layout.

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. $P_{\mathit{\text{ratio}}}$ is the ratio of the upper and lower surface material properties $P_u$ and $P_{\mathit{lo}}$ respectively.

The beam WFE has axial deformation $u_i$ and transverse deflection $v_i$ DOFs at all elemental nodes and rotation ${\theta }_i$ DOFs only present at elemental end nodes with corresponding axial forces $f_{\mathit{xi}}$, transverse forces $f_{\mathit{yi}}$ and bending moments ${\theta }_i$ as 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-2$ if 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 $\left\{{\mathit{h}}_{\mathit{e}}\right\}={\left\{\begin{array}{ccccccccccccc}{u}_1& v_1& {\theta }_1& u_2& v_2& u_3& v_3& \cdots & u_{r-1}& v_{r-1}& u_r& v_r& {\theta }_r\end{array}\right\}}^T$ and subsequently

$$\begin{array}{c}u\left(\xi \right)={}_{\left(1\times s\right)}{}^a\left\{{\Phi }_{z-2}^j\left(\xi \right)\right\}{}_{\left(s\times 1\right)}{}^{}\left\{c_e\right\}\\ v\left(\xi \right)={}_{\left(1\times s\right)}{}^t\left\{{\Phi }_z^j\left(\xi \right)\right\}{}_{\left(s\times 1\right)}{}^{}\left\{c_e\right\}\\ \theta \left(\xi \right)=\frac{\partial v\left(\xi \right)}{\partial x}=\frac{1}{L_e}\frac{\partial v\left(\xi \right)}{\partial \xi }=\frac{1}{L_e}{}_{\left(1\times s\right)}{}^t\left\{{\Phi }_z^{\prime j}\left(\xi \right)\right\}{}_{\left(s\times 1\right)}{}^{}\left\{c_e\right\}\end{array}$$E71

where the vector $\left\{{\mathit{c}}^{\mathit{e}}\right\}$ contains the unknown wavelet space element DOFs and

$$\begin{array}{c}{}_{\left(1\times s\right)}{}^{\mathrm{ a}}\left\{{\mathbf{\Phi }}_{\mathit{z}-2}^{\mathit{j}}\left(\xi \right)\right\}=\left\{\begin{array}{cccccccccc}{\varphi }_{z-2,h}^j\left(\xi \right)& 0& 0& {\varphi }_{z-2,h+1}^j\left(\xi \right)& 0& \cdots & 0& {\varphi }_{z-2,2^j-1}^j\left(\xi \right)& 0& 0\end{array}\right\}\\ {}_{\left(1\times s\right)}{}^{\mathrm{ t}}\left\{{\mathbf{\Phi }}_{\mathit{z}}^{\mathit{j}}\left(\xi \right)\right\}=\left\{\begin{array}{cccccccc}0& {\varphi }_{z,i}^j\left(\xi \right)& {\varphi }_{z,i+1}^j\left(\xi \right)& 0& \cdots & 0& {\varphi }_{z,2^j-2}^j\left(\xi \right)& {\varphi }_{z,2^j-1}^j\left(\xi \right)\end{array}\right\}\\ {}_{\left(1\times s\right)}{}^{\mathrm{ t}}\left\{{\mathbf{\Phi }}_{\mathit{z}}^{\prime \mathit{j}}\left(\xi \right)\right\}=\left\{\begin{array}{cccccccc}0& {\varphi }_{z,i}^{\prime j}\left(\xi \right)& {\varphi }_{z,i+1}^{\prime j}\left(\xi \right)& 0& \cdots & 0& {\varphi }_{z,2^j-2}^{\prime j}\left(\xi \right)& {\varphi }_{z,2^j-1}^{\prime j}\left(\xi \right)\end{array}\right\}\end{array}$$E72

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

$${}_{\left(s\times 1\right)}\left\{{\mathit{h}}_{\mathit{e}}\right\}={}_{\left(s\times s\right)}\left[{\mathit{R}}_{\mathit{p}}^{\mathit{w}}\right]{}_{\left(s\times 1\right)}\left\{{\mathit{c}}_{\mathit{e}}\right\}$$E73

and consequently

$$\begin{array}{c}u\left(\xi \right)={}_{\left(1\times s\right)}{}^a\left\{{\Phi }_z^j\left(\xi \right)\right\}{}_{\left(s\times s\right)}{}^{}\left[T_p^w\right]{}_{\left(s\times 1\right)}{}^{}\left\{h_e\right\}\\ v\left(\xi \right)={}_{\left(1\times s\right)}{}^t\left\{{\Phi }_z^j\left(\xi \right)\right\}{}_{\left(s\times s\right)}{}^{}\left[T_p^w\right]{}_{\left(s\times 1\right)}{}^{}\left\{h_e\right\}\\ \theta \left(\xi \right)=\frac{1}{L_e}{}_{\left(1\times s\right)}{}^t\left\{{\Phi }_z^{\prime j}\left(\xi \right)\right\}{}_{\left(s\times s\right)}{}^{}\left[T_p^w\right]{}_{\left(s\times 1\right)}{}^{}\left\{h_e\right\}\end{array}$$E74

The wavelet transformation matrix $\left[{\mathit{T}}_{\mathit{p}}^{\mathit{w}}\right]={\left[{\mathit{R}}_{\mathit{p}}^{\mathit{w}}\right]}^{-1}$. The strain energy of the functionally graded beam element, $U_e$, is defined as

$$\begin{array}{l}{U}_e=\frac{b}{2}\underset{-\frac{h}{2}}{\overset{\frac{h}{2}}{{\displaystyle \int }}}\underset{0}{\overset{1}{{\displaystyle \int }}}E\left(y\right)[\frac{1}{L_e}{\left(\frac{\partial u\left(\xi \right)}{\partial \xi }\right)}^T\left(\frac{\partial u\left(\xi \right)}{\partial \xi }\right)-\frac{y}{L_e{}^2}{\left(\frac{{\partial }^{2}v\left(\xi \right)}{\partial {\xi }^2}\right)}^T\left(\frac{\partial u\left(\xi \right)}{\partial \xi }\right)\\ -\frac{y}{L_e{}^2}{\left(\frac{\partial u\left(\xi \right)}{\partial \xi }\right)}^T\left(\frac{{\partial }^{2}v\left(\xi \right)}{\partial {\xi }^2}\right)+\frac{y^2}{L_e{}^3}{\left(\frac{{\partial }^{2}v\left(\xi \right)}{\partial {\xi }^2}\right)}^T\left(\frac{{\partial }^{2}v\left(\xi \right)}{\partial {\xi }^2}\right)]d\xi dy\end{array}$$E75

where $L_e$ is the length of the element and $E\left(y\right)$ the effective Young’s modulus obtained from Eq. (70). Let

$$\begin{array}{c}{}^A{E}_e={\int }_{-\frac{h}{2}}^{\frac{h}{2}}E\left(y\right)\mathit{dy}={\int }_{-\frac{h}{2}}^{\frac{h}{2}}{\left[E_u-E_l\right]\left(\frac{y}{h}+\frac{1}{2}\right)}^n+E_l\mathit{dy}\\ {}^B{E}_e={\int }_{-\frac{h}{2}}^{\frac{h}{2}}yE\left(y\right)\mathit{dy}={\int }_{-\frac{h}{2}}^{\frac{h}{2}}y\left({\left[E_u-E_l\right]\left(\frac{y}{h}+\frac{1}{2}\right)}^n+E_l\right)\mathit{dy}\\ {}^C{E}_e={\int }_{-\frac{h}{2}}^{\frac{h}{2}}{y}^{2}E\left(y\right)\mathit{dy}={\int }_{-\frac{h}{2}}^{\frac{h}{2}}\left(y^2\left[E_u-E_l\right]{\left(\frac{y}{h}+\frac{1}{2}\right)}^n+E_l\right)\mathit{dy}\end{array}$$E76

${}^A{E}_e$, ${}^B{E}_e$ and ${}^C{E}_e$ denote axial, axial-bending coupling and bending stiffness of the WFE respectively. The wavelet based physical space elemental stiffness matrix of the beam, $\left[{\mathit{k}}_{\mathit{e}}^{\mathit{w}}\right]$, is

$$\begin{array}{c}{}_{\left(s\times s\right)}{}^{\mathrm{ A}}\left[{\mathit{k}}_{\mathit{e}}^{\mathit{w}}\right]={\int }_0^1\frac{b{}^A{E}_e}{L_e}{\left[{\mathit{T}}_{\mathit{p}}^{\mathit{w}}\right]}^T{}^a{\left\{\frac{\partial {\mathbf{\Phi }}_{\mathit{z}-2}^{\mathit{j}}\left(\xi \right)}{\partial \xi }\right\}}^T{}^a\left\{\frac{\partial {\mathbf{\Phi }}_{\mathit{z}-2}^{\mathit{j}}\left(\xi \right)}{\partial \xi }\right\}\left[{\mathit{T}}_{\mathit{p}}^{\mathit{w}}\right]\mathit{d\xi }\hfill \\ {}_{\left(s\times s\right)}{}^{\mathrm{ B}}\left[{\mathit{k}}_{\mathit{e}}^{\mathit{w}}\right]={\int }_0^1\frac{b{}^B{E}_e}{{L_e}^2}{\left[{\mathit{T}}_{\mathit{p}}^{\mathit{w}}\right]}^T{}^t{\left\{\frac{{\partial }^2{\mathbf{\Phi }}_{\mathit{z}}^{\mathit{j}}\left(\xi \right)}{\partial {\xi }^2}\right\}}^T{}^a\left\{\frac{\partial {\mathbf{\Phi }}_{\mathit{z}-2}^{\mathit{j}}\left(\xi \right)}{\partial \xi }\right\}\left[{\mathit{T}}_{\mathit{p}}^{\mathit{w}}\right]\mathit{d\xi }\hfill \\ {}_{\left(s\times s\right)}{}^{\mathrm{ C}}\left[{\mathit{k}}_{\mathit{e}}^{\mathit{w}}\right]={\int }_0^1\frac{b{}^B{E}_e}{{L_e}^2}{\left[{\mathit{T}}_{\mathit{p}}^{\mathit{w}}\right]}^T{}^a{\left\{\frac{\partial {\mathbf{\Phi }}_{\mathit{z}-2}^{\mathit{j}}\left(\xi \right)}{\partial \xi }\right\}}^T{}^t\left\{\frac{{\partial }^2{\mathbf{\Phi }}_{\mathit{z}}^{\mathit{j}}\left(\xi \right)}{\partial {\xi }^2}\right\}\left[{\mathit{T}}_{\mathit{p}}^{\mathit{w}}\right]\mathit{d\xi }\hfill \\ {}_{\left(s\times s\right)}{}^{\mathrm{ D}}\left[{\mathit{k}}_{\mathit{e}}^{\mathit{w}}\right]={\int }_0^1\frac{b{}^C{E}_e}{{L_e}^3}{\left[{\mathit{T}}_{\mathit{p}}^{\mathit{w}}\right]}^T{}^t{\left\{\frac{{\partial }^2{\mathbf{\Phi }}_{\mathit{z}}^{\mathit{j}}\left(\xi \right)}{\partial {\xi }^2}\right\}}^T{}^t\left\{\frac{{\partial }^2{\mathbf{\Phi }}_{\mathit{z}}^{\mathit{j}}\left(\xi \right)}{\partial {\xi }^2}\right\}\left[{\mathit{T}}_{\mathit{p}}^{\mathit{w}}\right]\mathit{d\xi }\hfill \\ {}_{\left(s\times s\right)}\left[{\mathit{k}}_{\mathit{e}}^{\mathit{p}}\right]={}^A\left[{\mathit{k}}_{\mathit{e}}^{\mathit{p}}\right]-{}^B\left[{\mathit{k}}_{\mathit{e}}^{\mathit{p}}\right]-{}^C\left[{\mathit{k}}_{\mathit{e}}^{\mathit{p}}\right]+{}^D\left[{\mathit{k}}_{\mathit{e}}^{\mathit{p}}\right]\hfill \end{array}$$E77

The kinetic energy of the functionally graded beam element, ${\Lambda }_e$, is defined as

$$\begin{array}{c}{\Lambda }_e=\frac{1}{2}\underset{o}{\overset{b}{{\displaystyle \int }}}dz\underset{-\frac{h}{2}}{\overset{\frac{h}{2}}{{\displaystyle \int }}}\underset{0}{\overset{1}{{\displaystyle \int }}}\rho \left(y\right)(L_e\left(\dot{u}\left(\xi \text{,}t\right)\dot{u}\left(\xi \text{,}t\right)\right)-y\left(\dot{u}\left(\xi \text{,}t\right)\frac{\partial \dot{v}\left(\xi \text{,}t\right)}{\partial x}\right)-y\left(\frac{\partial \dot{v}\left(\xi \text{,}t\right)}{\partial \xi }\dot{u}\left(\xi \text{,}t\right)\right)\\ +\frac{y^2}{L_e}\left(\frac{\partial \dot{v}\left(\xi \text{,}t\right)}{\partial x}\frac{\partial \dot{v}\left(\xi \text{,}t\right)}{\partial x}\right)+L_e\left(\dot{v}\left(\xi \text{,}t\right)\dot{v}(\xi \text{,}t)\right))d\xi dy\end{array}$$E78

$\rho \left(y\right)$ is the effective density also obtained from Eq. (70). Let the inertial coefficients be denoted as:

$$\begin{array}{c}{}^A{\rho }_e={\int }_{-\frac{h}{2}}^{\frac{h}{2}}\rho \left(y\right)\mathit{dy}={\int }_{-\frac{h}{2}}^{\frac{h}{2}}{\left[{\rho }_u-{\rho }_l\right]\left(\frac{y}{h}+\frac{1}{2}\right)}^n+{\rho }_l\mathit{dy}\\ {}^B{\rho }_e={\int }_{-\frac{h}{2}}^{\frac{h}{2}}\mathit{y\rho }\left(y\right)\mathit{dy}={\int }_{-\frac{h}{2}}^{\frac{h}{2}}y\left({\left[{\rho }_u-{\rho }_l\right]\left(\frac{y}{h}+\frac{1}{2}\right)}^n+{\rho }_l\right)\mathit{dy}\\ {}^C{\rho }_e={\int }_{-\frac{h}{2}}^{\frac{h}{2}}{y}^2\rho \left(y\right)\mathit{dy}={\int }_{-\frac{h}{2}}^{\frac{h}{2}}\left(y^2\left[{\rho }_u-{\rho }_l\right]{\left(\frac{y}{h}+\frac{1}{2}\right)}^n+{\rho }_l\right)\mathit{dy}\end{array}$$E79

The wavelet based physical space elemental mass matrix of the beam, $\left[{\mathit{m}}_{\mathit{e}}^{\mathit{p}}\right]$, is

$$\begin{array}{c}{}_{\left(s\times s\right)}{}^{\mathrm{ A}}\left[{\mathit{m}}_{\mathit{e}}^{\mathit{w}}\right]={\int }_0^{1}b{}^A{\rho }_e{L}_e{\left[{\mathit{T}}_{\mathit{p}}^{\mathit{w}}\right]}^T{}^a{\left\{{\mathbf{\Phi }}_{\mathit{z}-2}^{\mathit{j}}\left(\xi \right)\right\}}^T{}^a\left\{{\mathbf{\Phi }}_{\mathit{z}-2}^{\mathit{j}}\left(\xi \right)\right\}{\left[{\mathit{T}}_{\mathit{p}}^{\mathit{w}}\right]}^T\mathit{d\xi }\hfill \\ {}_{\left(s\times s\right)}{}^{\mathrm{ B}}\left[{\mathit{m}}_{\mathit{e}}^{\mathit{w}}\right]={\int }_0^{1}b{}^B{\rho }_e{\left[{\mathit{T}}_{\mathit{p}}^{\mathit{w}}\right]}^T{}^a{\left\{{\mathbf{\Phi }}_{\mathit{z}-2}^{\mathit{j}}\left(\xi \right)\right\}}^T{}^t\left\{\frac{\partial {\mathbf{\Phi }}_{\mathit{z}}^{\mathit{j}}\left(\xi \right)}{\partial \xi }\right\}{\left[{\mathit{T}}_{\mathit{p}}^{\mathit{w}}\right]}^T\mathit{d\xi }\hfill \\ {}_{\left(s\times s\right)}{}^{\mathrm{ C}}\left[{\mathit{m}}_{\mathit{e}}^{\mathit{w}}\right]={\int }_0^{1}b{}^B{\rho }_e{\left[{\mathit{T}}_{\mathit{p}}^{\mathit{w}}\right]}^T{}^t{\left\{\frac{\partial {\mathbf{\Phi }}_{\mathit{z}}^{\mathit{j}}\left(\xi \right)}{\partial \xi }\right\}}^T{}^a\left\{{\mathbf{\Phi }}_{\mathit{z}-2}^{\mathit{j}}\left(\xi \right)\right\}{\left[{\mathit{T}}_{\mathit{p}}^{\mathit{w}}\right]}^T\mathit{d\xi }\hfill \\ {}_{\left(s\times s\right)}{}^{\mathrm{ D}}\left[{\mathit{m}}_{\mathit{e}}^{\mathit{w}}\right]={\int }_0^{1}b{}^C{\rho }_e{\left[{\mathit{T}}_{\mathit{p}}^{\mathit{w}}\right]}^T{}^t{\left\{\frac{\partial {\mathbf{\Phi }}_{\mathit{z}}^{\mathit{j}}\left(\xi \right)}{\partial \xi }\right\}}^T{}^t\left\{\frac{\partial {\mathbf{\Phi }}_{\mathit{z}}^{\mathit{j}}\left(\xi \right)}{\partial \xi }\right\}{\left[{\mathit{T}}_{\mathit{p}}^{\mathit{w}}\right]}^T\mathit{d\xi }\hfill \\ {}_{\left(s\times s\right)}{}^{\mathrm{ E}}\left[{\mathit{m}}_{\mathit{e}}^{\mathit{w}}\right]={\int }_0^{1}b{}^A{\rho }_e{L}_e{\left[{\mathit{T}}_{\mathit{p}}^{\mathit{w}}\right]}^T{}^t{\left\{{\overline{\mathbf{\Phi }}}_{\mathit{z}}^{\mathit{j}}\left(\xi \right)\right\}}^T{}^t\left\{{\overline{\mathbf{\Phi }}}_{\mathit{z}}^{\mathit{j}}\left(\xi \right)\right\}{\left[{\mathit{T}}_{\mathit{p}}^{\mathit{w}}\right]}^T\mathit{d\xi }\hfill \\ {}_{\left(s\times s\right)}\left[{\mathit{m}}_{\mathit{e}}^{\mathit{p}}\right]={}^A\left[{\mathit{m}}_{\mathit{e}}^{\mathit{p}}\right]-{}^B\left[{\mathit{m}}_{\mathit{e}}^{\mathit{p}}\right]-{}^C\left[{\mathit{m}}_{\mathit{e}}^{\mathit{p}}\right]+{}^D\left[{\mathit{m}}_{\mathit{e}}^{\mathit{p}}\right]+{}^E\left[{\mathit{m}}_{\mathit{e}}^{\mathit{p}}\right]\hfill \end{array}$$E80

4. Numerical examples

Example 1: A uniform axial cantilever rod (free-fixed) subjected to linear varying load $q\left(x\right)=-q_{0}x$ has a uniform cross sectional area,$A=A_0$, Young’s Modulus, $E=E_0$ and length l. The exact solution for displacement at a particular point x can be obtained by solving $u\left(x\right)=\frac{1}{\mathit{EA}}\underset{o}{\overset{x}{\int }}P\left(x\right)\mathit{dx}=\frac{1}{E_0{A}_0}\int q_0\frac{x^2}{2}\mathit{dx}$ [40]. One WFE is used to represent the rod using Daubechies and BSWI WFEM approaches and the results are compared with the exact, h-FEM and p-FEM formulations. The governing equation of the system for FEM and WFEM is

$$\left[{\mathit{K}}_{\mathit{r}}\right]\left\{{\mathit{U}}_{\mathit{r}}\right\}=\left\{{\mathit{F}}_{\mathit{r}}\right\}$$E81

where $\left[{\mathit{K}}_{\mathit{r}}\right]$ is the system stiffness matrix, $\left\{{\mathit{U}}_{\mathit{r}}\right\}$ is the system vector containing the DOFs and $\left\{{\mathit{F}}_{\mathit{r}}\right\}$ is the loading vector of the system. The axial deformation of the rod is analysed at the arbitrary point 0.1l and the rate of convergence of the different approaches is compared in Figure 5. The plot shows the absolute relative error of the axial deformation and corresponding number of DOFs. The FEM (h-FEM) solution involves increasing the number of elements, p-FEM involves increasing the order of the polynomials (one element only) and both Daubechies and BSWI WFEMs have the order and/or multiresolution scale j increased. The results show that although the rates of convergence of all the methods are similar, the WFEM approaches have a slightly improved rate with only one element employed.

Example 2: A simply supported two-stepped beam of length 2l has non-uniform flexural stiffness represented by the unequal cross sections; the bending stiffness of the right and left half is given as $E_1{I}_1=E_0{I}_0$ and $E_2{I}_2=4E_0{I}_0$ respectively. The entire beam is subjected to a uniformly distributed load q(x) = 1. The flexural stiffness function is expressed as [44]:

$$E\left(x\right)I\left(x\right)=E_0{I}_0\left[1-\gamma \widehat{H}\left(x-x_0\right)\right]$$E82

where $\gamma =0.75$ is defined as the decrement of discontinuity intensity and satisfies the condition $0\le \gamma \le 1$to ensure positivity of the flexural stiffness. $\widehat{H}\left(x-x_0\right)$ is the Heaviside function for $0\le x_0\le 2l$. The general analytical governing equation is

$${\left\{E_0{I}_0\left[1-\gamma \widehat{H}\left(x-x_0\right)\right]v^{\prime \prime }\left(x\right)\right\}}^{\prime \prime }=q\left(x\right)$$E83

The FEM and WFEM governing equation is summarised as:

$$\left[{\mathit{K}}_{\mathit{b}}\right]\left\{{\mathit{V}}_{\mathit{b}}\right\}=\left\{{\mathit{F}}_{\mathit{b}}\right\}$$E84

The vector $\left\{{\mathit{V}}_{\mathit{b}}\right\}$ contains the system DOFs within the entire beam, $\left[{\mathit{K}}_{\mathit{b}}\right]$ is the beam stiffness matrix and $\left\{{\mathit{F}}_{\mathit{b}}\right\}$ is the equivalent system load vector. The h-FEM (FEM-8; 8 elements), p-FEM of order 9 (p-FEM-9; 2 elements), Daubechies WFEM of order $L=10$ and scale $j=1$ (D10₁; 2 elements) and the BSWI WFEM of order $m=3$ and scale $j=3$ (BSWI3₃; 2 elements) are selected for comparison with the exact solution governed by Eq. (83). Each approach has a total of 18 DOFs within the beam. The deflection and rotation across the beam is presented in Figure 6(a) and (b) respectively. The percentage errors of the deflections are compared for the different approaches and presented in Figure 6(c). All numerical approaches describe the deflection and rotation across the beams very accurately. However, given that both the Daubechies and BSWI WFEM deflection solutions have a maximum error of 1.28% in comparison to 3.82% from the h-FEM and p-FEM approaches, the WFEMs exhibit better convergence to the exact solution. Furthermore, improved accuracy is attained with fewer elements implemented than the h-FEM and p-FEM and this results in reduced computational time.

Example 3: A steel-alumina functionally graded beam of length l and uniform cross-sectional area $A=0.36$ m² (height $h=0.9$ m and width $b=0.4$ m) is fully alumina at the upper surface and fully steel at the lower surface with material properties; $E_u=3.9\times {10}^{11}$ Pa, ${\rho }_u=3.96\times {10}^3$ kg·m⁻³ and $E_l=2.1\times {10}^{11}$ Pa, ${\rho }_l=7.8\times {10}^3$ kg·m⁻³ respectively ($E_{\text{ratio}}=\frac{E_u}{E_l}$; ${\rho }_{\text{ratio}}=\frac{{\rho }_u}{{\rho }_l}$). E and ρ denote the Young’s modulus and density respectively. The slenderness ratio for the beam is $\raisebox{1ex}{$l$}\!\left/ \!\raisebox{-1ex}{$h$}\right.=100$. The free vibration of the steel-alumina beam is analysed for the boundary conditions pinned-pinned (PP), pinned-clamped (PC), clamped-clamped (CC) and clamped-free (CF), for different values of n in Eq. (70). The free vibration of the functionally graded beam is governed by [45]

$$\left[\left[\mathit{K}\right]-{\omega }^2\left[\mathit{M}\right]\right]\left\{\acute{\mathit{U}}\right\}=0$$E85

The matrices $\left[\mathit{K}\right]$ and $\left[\mathit{M}\right]$ are the stiffness and mass matrices for the functionally graded beam, $\omega$ is the natural frequency and $\left\{\acute{\mathit{U}}\right\}$ is the vector containing the DOFs within the entire beam. The i^th non-dimensional frequency ${\lambda }_i$ of the FGM beam is evaluated from the relation ${{\lambda }_i}^2={\omega }_i{l}^2{\left(\frac{12{\rho }_l}{E_l{h}^2}\right)}^{\frac{1}{2}}$. The functionally graded beam is modelled for the different approaches using 2 Daubechies WFEs ($L=12$; $j=0$; 37 DOFs); one BSWI ($m=5$; $j=4$; 38 DOFs) WFE and 12 h-FEM elements (39 DOFs). The results of the first 3 non-dimensional natural frequencies of the beam are presented in Table 1 for different boundary conditions and material distributions. It is observed that all approaches give highly accurate results with respect to the reference (BSWI5₅), particularly for the fundamental frequencies. Furthermore, the BSWI WFEM solution exhibits better levels of accuracy than the Daubechies WFEM and h-FEM solutions for the higher frequencies. Both WFEM solutions achieve high levels of accuracy with the described layout of having the rotation DOFs present at elemental and nodes and using fewer elements that the h-FEM approach.

Figure 5.

Comparison of the convergence of the axial deformation at point x = 0.1l.

Figure 6.

(a) Deflection and (b) rotation (c) comparison of the deflection percentage error across a simply supported stepped beam subjected to a uniformly distributed load q(x) = 1.

			n = 0	n = 0.1	n = 0.5	n = 1	n = 5	n = 10	n = 10^4
${\mathit{\lambda }}_1$	PP	BSWI55	4.34462	4.1943	3.84903	3.65795	3.37139	3.29504	3.33251
		FEM	4.34463	4.19431	3.84912	3.65811	3.3715	3.2951	3.33258
		D120	4.34462	4.1943	3.84903	3.65795	3.37139	3.29504	3.33251
		BSWI54	4.34462	4.1943	3.84903	3.65795	3.37139	3.29504	3.33251
	PC	BSWI55	5.43022	5.24233	4.81079	4.57197	4.21381	4.11838	3.92681
		FEM	5.43024	5.24238	4.81112	4.57253	4.21419	4.11856	3.92682
		D120	5.43023	5.24234	4.8108	4.57197	4.21382	4.11839	3.92681
		BSWI54	5.43022	5.24233	4.81079	4.57197	4.21381	4.11839	3.92681
	CC	BSWI55	6.54131	6.31498	5.79514	5.50745	5.07601	4.96105	4.73028
		FEM	6.54137	6.31509	5.79585	5.50867	5.07685	4.96145	4.73028
		D120	6.54132	6.31498	5.79514	5.50745	5.07601	4.96106	4.73028
		BSWI54	6.54131	6.31498	5.79514	5.50745	5.07601	4.96105	4.73028
	CF	BSWI55	2.59318	2.50345	2.29737	2.18333	2.01229	1.96671	1.87523
		FEM	2.59318	2.50346	2.2974	2.18337	2.01232	1.96673	1.87523
		D120	2.59318	2.50345	2.29737	2.18333	2.01229	1.96671	1.87523
		BSWI54	2.59318	2.50345	2.29737	2.18333	2.01229	1.96671	1.87523
${\mathit{\lambda }}_2$	PP	BSWI55	8.68871	8.38806	7.69754	7.31541	6.74237	6.58969	6.66461
		FEM	8.68894	8.38834	7.69844	7.31684	6.74338	6.59024	6.66537
		D120	8.68968	8.389	7.6984	7.31623	6.74313	6.59043	6.66535
		BSWI54	8.68871	8.38806	7.69754	7.31541	6.74238	6.58969	6.66461
	PC	BSWI55	9.77473	9.4365	8.65966	8.22977	7.58512	7.41335	7.06849
		FEM	9.77513	9.43702	8.66145	8.23263	7.58713	7.41442	7.06879
		D120	9.7765	9.43821	8.66124	8.23127	7.5865	7.4147	7.06977
		BSWI54	9.77473	9.4365	8.65967	8.22977	7.58512	7.41335	7.06849
	CC	BSWI55	10.8597	10.4839	9.62083	9.14322	8.42702	8.23619	7.85305
		FEM	10.8604	10.4848	9.62387	9.14808	8.43044	8.238	7.85355
		D120	10.8636	10.4877	9.62433	9.14655	8.43009	8.23918	7.8559
		BSWI54	10.8597	10.4839	9.62083	9.14322	8.42702	8.23619	7.85305
	CF	BSWI55	6.49133	6.26671	5.75083	5.46534	5.03722	4.92315	4.69413
		FEM	6.49138	6.2668	5.75131	5.46615	5.03778	4.92342	4.69417
		D120	6.49134	6.26673	5.75084	5.46535	5.03723	4.92316	4.69415
		BSWI54	6.49133	6.26671	5.75083	5.46534	5.03722	4.92315	4.69413
${\mathit{\lambda }}_3$	PP	BSWI55	13.0317	12.5808	11.545	10.9719	10.1125	3.29504	3.33251
		FEM	13.0334	12.5826	11.5489	10.9774	10.1166	3.2951	3.33258
		D120	13.0461	13.5947	11.5578	10.984	10.1237	3.29504	3.33251
		BSWI54	13.0317	12.5808	11.545	10.9719	10.1125	3.29504	3.33251
	PC	BSWI55	14.1176	13.629	12.507	11.8861	10.9551	10.7071	10.209
		FEM	14.1201	13.6318	12.5131	11.895	10.9617	10.711	10.2108
		D120	14.1444	13.655	12.5308	11.9088	10.9761	10.7275	10.2284
		BSWI54	14.1176	13.629	12.507	11.8861	10.9551	10.7071	10.209
	CC	BSWI55	15.2034	14.6773	13.4689	12.8003	11.7977	11.5306	10.9942
		FEM	15.2071	14.6813	13.4779	12.8135	11.8074	11.5364	10.9968
		D120	15.2662	14.7379	13.5247	12.8533	11.8466	11.5783	11.0396
		BSWI54	15.2034	14.6773	13.469	12.8003	11.7978	11.5306	10.9942
	PC	BSWI55	10.8611	10.4853	9.62205	9.14437	8.42814	8.2373	7.85409
		FEM	10.8618	10.4861	9.62437	9.14797	8.4307	8.23872	7.85458
		D120	10.8667	10.4907	9.62702	9.1491	8.43249	8.24155	7.85814
		BSWI54	10.8611	10.4853	9.62206	9.14438	8.42815	8.2373	7.8541

Table 1.

The non-dimensional frequencies of a steel-alumina FG beam for different transverse varying distributions and boundary conditions.

Assume the same beam, with simply supported boundary conditions and length $l=20$ m, is subjected to a moving load of magnitude $P=1\times {10}^5$ N travelling across at c m·s⁻¹. The behaviour of the beam is described using Euler Bernoulli beam theory and is assumed to be undamped. The governing equation describing the dynamic behaviour of the system is given by [45]:

$$\left[\mathbf{M}\right]\left\{\ddot{\mathbf{U}}\left(t\right)\right\}+\left[\mathbf{K}\right]\left\{\mathbf{U}\left(t\right)\right\}=\left\{\mathbf{F}\left(t\right)\right\}$$E86

where $\left\{\ddot{\mathbf{U}}\left(t\right)\right\}$ and $\left\{\mathbf{U}\left(t\right)\right\}$ represent the system acceleration and displacement vectors at time t. $\left\{\mathbf{F}\left(t\right)\right\}$ is the moving load vector. The deflection of the beam $v\left(x,t\right)$, as the moving load travels across, is normalised as a non-dimensional parameter $v\left(x,t\right)/v_0$ where $v_0=\frac{P{l}^3}{48E_{l}I}$ is the deflection at the centre of the simply supported functionally graded beam when subjected to a static load of magnitude P at the centre. The maximum normalised deflection mid-span of the beam is analysed over a moving load velocity range $0<c\le 300$ m·s⁻¹ at increments of 1 m·s⁻¹ to identify the critical velocity for the different variations of the constituent materials as illustrated in Figure 7. The results present are obtained from the BSWI (2 element; $m=4$; $j=3$; 37 DOFs) WFEM solution. The h-FEM (12 elements; 39 DOFs) and Daubechies (2 elements; $L=12$; $j=0$; 37 DOFs) WFEM solution gives similar results. The values of the critical moving load velocity and corresponding maximum non-dimensional displacement are presented in Table 2 for the different values of n for all approaches. The results are compared with those presented in [46]. Both the Daubechies and BSWI WFE M solutions very accurately yield the correct values.

Figure 7.

Variation of the maximum non-dimensional vertical displacement at the centre of a simply supported steel-alumina beam with respect to moving load velocities, for different n.

n	Critical velocity c m∙s^−1				Max$\left[\frac{v\left(\frac{l}{2},t\right)}{v_0}\right]$
	Ref. [46]	FEM	D120	BSWI43	Ref. [46]	FEM	D120	BSWI43
0	252	252	252	252	0.9328	0.9322	0.9323	0.9322
0.1	–	235	235	235	–	0.9863	0.9864	0.9863
0.2	222	222	222	222	1.0344	1.0340	1.0340	1.0340
0.5	198	198	198	198	1.1444	1.1435	1.1437	1.1436
1	179	178	178	178	1.2503	1.2491	1.2495	1.2493
2	164	164	164	164	1.3376	1.3363	1.3368	1.3365
3	–	157	158	158	–	1.3747	1.3751	1.3748
5	–	151	151	152	–	1.4217	1.422	1.4218
7	–	148	148	148	–	1.4567	1.4570	1.4568
10	–	145	145	145	–	1.4974	1.4976	1.4974
10^4	132	132	132	132	–	1.7308	1.7309	1.7308

Table 2.

The non-dimensional frequencies of a steel-alumina FG beam for different transverse varying distributions and boundary conditions.

5. Conclusions

A generalised formulation framework for the construction of an axial rod, Euler Bernoulli beam and functionally graded two-dimensional wavelet based finite elements is presented. The Daubechies and BSWI families are selected due to their desirable properties, particularly compact support, ‘two-scale’ relation and multiresolution. It is illustrated via a set of numerical examples that the WFEMs perform exceptionally well when compared to conventional h-FEM and p-FEM where high levels of accuracy are achieved with fewer elements required and the approaches converge more rapidly to the exact solution. Furthermore, the methods are able to accurately describe the behaviour of static and dynamic systems with singularities, variation in material properties and loading conditions present. This exhibits the vast potential of the method in the analysis of more complicated systems and the ability to alter the multiresolution scales without affecting the original mesh allows effective and efficient avenues solution accuracy improvement.