﻿ A Study of the Elastodynamic Problem by Meshless Local Petrov-Galerkin Method Using the Laplace-Transform

World Journal of Mechanics
Vol.08 No.02(2018), Article ID:82519,13 pages
10.4236/wjm.2018.82004

A Study of the Elastodynamic Problem by Meshless Local Petrov-Galerkin Method Using the Laplace-Transform

Department of Physics, Faculty of Sciences, Moulay Ismail University, Meknes, Morocco    Received: December 19, 2017; Accepted: February 11, 2018; Published: February 14, 2018

ABSTRACT

The Meshless Local Petrov-Galerkin (MLPG) with Laplace transform is used for solving partial differential equation. Local weak form is developed using the weighted residual method locally from the dynamic partial differential equation and using the moving least square (MLS) method to construct shape function. This method is a more effective alternative than the finite element method for computer modelling and simulation of problems in engineering; however, the accuracy of the present method depends on a number of parameters deriving from local weak form and different subdomains. In this paper, the meshless local Petrov-Galerkin (MLPG) formulation is proposed for forced vibration analysis. First, the results are presented for different values of ${\alpha }_{s}$ , and ${\alpha }_{q}$ with regular distribution of nodes ${n}_{t}=55$ . After, the results are presented with fixed values of ${\alpha }_{s}$ and ${\alpha }_{q}$ for different time-step.

Keywords:

Meshless, MLPG, Weak Form, MLS, PDEs, Elastodynamic, Laplace Transform, Support Domain, Quadrature Domain, Regular Distribution 1. Introduction

Generality of physical or mechanical problems are modeled by partial differential equations (PDEs). Moreover, a high number of numerical methods and techniques have been developed for the approximation of solution of the PDEs in much research that focused mainly on the improvement of accuracy and efficiency of them. In recent years, attention has turned on the development of meshless methods, especially for the numerical solution of partial differential equations. The meshless local Petrov-Galerkin (MLPG) approach  -  has become very attractive, as a promising method for solving problems. This method based on a local weak form and Moving Least Squares (MLS) approximation     . The main advantage of this method over the widely used ﬁnite element methods (FEM) is that it does not need any mesh, either for the interpolation of the solution variables or for the integration of the weak forms   . The MLPG method and its variations have been applied for elastodynamic and elastostatic problems by several authors. For example, it is applied to free and forced vibration analysis for solids by Gu Y. T. et al.  they found that the parameter which decides the size of the sub-domain needs to be chosen carefully. Long S. Y. et al., have applied MLPG5 method used the Heaviside function as the test function for elastic dynamic problems  . They found a good agreement compared with the results obtained by (FEM). This method also applied by Ping Xia et al. in elastic dynamic analysis of moderately thick plate using meshless LRPIM  . They have used the Newmark method for solving the dynamic problem and have studied the eﬀects of the size of the quadrature subdomain and the inﬂuence domain on the dynamic properties. They found if appropriate sizes are selected, good results and stability can be obtained. A meshfree-based local Galerkin method with condensation of degree of freedom for elastic dynamic analysis by De-An Hu et al.  , they have used the standard implicit Newmark’s time integration scheme for solving the global dynamic system equations obtained by assembling all local discrete equations. Recent results founded by Moussaoui and al concerning the effects of support domain ${\alpha }_{s}$ for elastostatic problem by MLPG  .

The MLPG formulation is proposed in this paper to extend the MLPG method to dynamic analysis and for solving the problem of a thin elastodynamic homogeneous rectangular plate  . The Laplace transform  is applied to eliminate the time variable, then, the obtained equations by the local formulation becomes in function with coefficient of Laplace transform. The Stehfest inversion method is applied to obtain the time-dependent solutions  . The result presented for different values of ${\alpha }_{s}$ and ${\alpha }_{q}$ with regular distribution of nodes ${n}_{t}=55$ . After, the results are presented with fixed values of ${\alpha }_{s}$ and ${\alpha }_{q}$ for different time-step. We found large domains of ${\alpha }_{s}$ , ${\alpha }_{q}$ and time-step by using Laplace transform method. The integral equations have a very simple non-singular form. Moreover, both the contour and domain integrations can be easily carried out in rectangular sub-domains.

This paper is organized as follows: Section 2 introduces the least squares approximation (MLS) for combination of shape function. In Section 3, the Basic equations of elastodynamics and their Laplace transforms are proposed. In Section 4, the MLPG formula including the local weak form in Laplace-transformed domain are developed, using the weighted residual method locally from the dynamic partial differential equation. The numerical results and discussions for 2D problem example are given in Section 5. Finally, the paper ends with the conclusion.

2. Moving Least Square (MLS) Approximation

We Consider a sub-domain ${\Omega }_{s}$ the neighbourhood of point ${X}^{\text{T}}=\left(x,y\right)$ , which is located within the problem domain $\Omega$ . To approximate a function $u\left(X,t\right)$ in ${\Omega }_{s}$ , a finite set of monomial basis functions $P\left(X\right)$ , is considered in the space coordinates $X$ in two-dimension is given by:

${P}^{\text{T}}\left(X\right)=\left[1,x,y\right]$ (1)

The approximation function of a field variable $u\left(X,t\right)$ is defined in a subdomain ${\Omega }_{s}$ by:

${u}^{h}\left(X,t\right)={\sum }_{j=1}^{m}{p}_{j}\left(X\right){a}_{j}\left(X,t\right)={P}^{\text{T}}\left(X\right)a\left(X,t\right)$ (2)

where $P\left(X\right)$ is the monomial basis function of the spatial coordinates ${X}^{\text{T}}=\left(x,y\right)$ for two dimensional problem, and m is the number of the monomial basis functions. $a\left(X,t\right)$ is a function of point $X$ and time $t$ and is a vector of coefficients ${a}_{i}\left(X,t\right)$ given by:

$a\left(X,t\right)=\left\{{a}_{1}\left(X,t\right)\text{\hspace{0.17em}}\text{\hspace{0.17em}}{a}_{2}\left(X,t\right)\text{\hspace{0.17em}}\text{\hspace{0.17em}}{a}_{3}\left(X,t\right)\text{\hspace{0.17em}}\text{\hspace{0.17em}}\cdots \text{\hspace{0.17em}}\text{\hspace{0.17em}}{a}_{m}\left(X,t\right)\right\}$ (3)

The coefficient $a\left(X,t\right)$ is obtained at any point $X$ by minimizing a weigh- ted discrete L2 norm

$J={\sum }_{i=1}^{n}R\left(X-{X}_{i}\right){\left[{P}^{\text{T}}\left({X}_{i}\right)a\left(X,t\right)-{u}_{i}\left(t\right)\right]}^{2}$ (4)

where n is the number of nodes in the support domain of $X$ for which the weight function $R\left(X-{X}_{i}\right)\ne 0$ , and ${u}_{i}$ is the nodal parameter of $u$ at $X={X}_{i}$ .

The stationarity of $J$ with respect to $a\left(X,t\right)$ gives:

$\frac{\partial J}{\partial a}=0$ (5)

which leads to the following linear relations:

$A\left(X\right)a\left(X,t\right)=B\left(X\right)u\left(t\right)$ (6)

where

$u\left(t\right)$ is the vector that collects the nodal displacements for all nodes in the support domain:

$u\left(t\right)={\left\{{u}_{1}\left(t\right)\text{\hspace{0.17em}}\text{\hspace{0.17em}}{u}_{2}\left(t\right)\text{\hspace{0.17em}}\text{\hspace{0.17em}}\cdots \text{\hspace{0.17em}}\text{\hspace{0.17em}}{u}_{n}\left(t\right)\right\}}^{\text{T}}$ (7)

$A\left(X\right)$ is called the weighted moment matrix defined by:

$A\left(X\right)={\sum }_{i=1}^{n}{R}_{i}\left(X\right)p\left({X}_{i}\right){p}^{\text{T}}\left({X}_{i}\right)$ (8)

and the matrix B is defined by:

$B\left(X\right)=\left[{R}_{1}\left(X\right)p\left({X}_{1}\right)\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{R}_{2}\left(X\right)p\left({X}_{2}\right)\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\cdots \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{R}_{n}\left(X\right)p\left({X}_{n}\right)\right]$ (9)

Solving $a\left(X,t\right)$ from Equation (6) as:

$a\left(X,t\right)={A}^{-1}\left(X\right)B\left(X\right)u\left(t\right)$ (10)

and substituting the above Equation (10) back into Equation (2) we obtain:

${u}^{h}\left(X,t\right)={\sum }_{i=1}^{n}{\phi }_{i}\left(X\right){u}_{i}\left(t\right)={\Phi }^{\text{T}}\left(X\right)u\left(t\right)$ (11)

where ${\Phi }^{\text{T}}\left(X\right)$ is the matrix of MLS shape functions corresponding n nodes in the support domain of the point $X$ and can be written as:

${\Phi }^{\text{T}}\left(X\right)={\left\{{\phi }_{1}\left(X\right)\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\phi }_{2}\left(X\right)\text{\hspace{0.17em}}\text{\hspace{0.17em}}\cdots \text{\hspace{0.17em}}\text{\hspace{0.17em}}{\phi }_{n}\left(X\right)\right\}}_{\left(1×n\right)}={p}^{\text{T}}\left(X\right){A}^{-1}\left(X\right)B\left(X\right)$ (12)

The shape function ${\varphi }_{i}\left(X\right)$ for the ith node is defined by:

${\varphi }_{i}\left(X\right)={\sum }_{j=1}^{m}{p}_{j}\left(X\right){\left({A}^{-1}\left(X\right)B\left(X\right)\right)}_{ji}={p}^{\text{T}}\left(X\right){\left({A}^{-1}B\right)}_{i}$ (13)

The following quartic spline function is used, it has the following form of

${R}_{i}\left(r\right)=\left\{\begin{array}{l}1-6{r}_{i}^{2}+8{r}_{i}^{3}-3{r}_{i}^{4}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{r}_{i}\le 1\\ 0\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{ }\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{r}_{i}\ge 1\end{array}$ (14)

${r}_{i}=\frac{|X-{X}_{i}|}{{d}_{i}}$

In which $|X-{X}_{i}|$ is the distance from node ${X}_{i}$ to point $X$ , and ${d}_{i}$ is the size of the influence domain for the weight function.

3. Basic Equations of Elastodynamics

The governing equations for a linear two-dimensional elastodynamic problem on a domain Ω bounded by a boundary Γ are:

${\sigma }_{ij,j}+{b}_{i}=\rho {\stackrel{¨}{u}}_{i}+c{\stackrel{˙}{u}}_{i}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{in}\text{\hspace{0.17em}}\Omega$ (15)

${\epsilon }_{ij}=\frac{\left({u}_{i,j}+{u}_{j,i}\right)}{2}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{in}\text{\hspace{0.17em}}\Omega$ (16)

${\sigma }_{ij}={D}_{ijkl}{\epsilon }^{ijkl}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{in}\text{\hspace{0.17em}}\Omega$ (17)

where $\rho$ the mass density, c is the damping coefficient, ${\stackrel{¨}{u}}_{i}=\frac{{\partial }^{2}{u}_{i}}{\partial {t}^{2}}$ is the acceleration, ${\stackrel{˙}{u}}_{i}=\frac{\partial {u}_{i}}{\partial t}$ the velocity, ${\sigma }_{ij}$ the stress tensor, ${b}_{i}$ the body force tensor, and ${\left(\text{ }\right)}_{,j}$ denotes $\frac{\partial }{\partial {X}_{j}}$

The initial and boundary conditions are given as follows:

${u}_{i}={\stackrel{^}{u}}_{i}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{on}\text{\hspace{0.17em}}{\Gamma }_{u}$ (18)

${t}_{i}={\sigma }_{ij}{n}_{j}={\stackrel{^}{t}}_{i}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{on}\text{\hspace{0.17em}}{\Gamma }_{t}$ (19)

${u}_{1}\left(X,{t}_{0}\right)={u}_{i0}\left(X\right)\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}X\in \Omega$ (20)

${\stackrel{˙}{u}}_{2}\left(X,{t}_{0}\right)={v}_{i0}\left(X\right)\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}X\in \Omega$ (21)

where ${\stackrel{^}{u}}_{i}$ and ${\stackrel{^}{t}}_{i}$ are the prescribed displacement and traction on the boundary ${\Gamma }_{u}$ and ${\Gamma }_{t}$ respectively, ${u}_{i0}$ and ${v}_{i0}$ denote the initial displacement and initial velocity, ${n}_{j}$ is the unit outward normal to the boundary $\Gamma$ . ${\Gamma }_{u}$ and ${\Gamma }_{t}$ are complementary subsets of $\Gamma$ .

The Laplace-transform  of a function $f\left(x,t\right)$ is deﬁned by:

$L\left(f\left(x,t\right)\right)={\int }_{0}^{+\infty }f\left(x,t\right){\text{e}}^{-st}\text{d}\tau =\stackrel{¯}{f}\left(x,s\right)$ (22)

$L\left({f}^{\prime }\left(x,t\right)\right)=sL\left(f\left(x,s\right)\right)-f\left(x,0\right)$ (23)

$L\left({f}^{″}\left(x,t\right)\right)={s}^{2}L\left(f\left(x,s\right)\right)-sf\left(x,0\right)-{f}^{\prime }\left(x,0\right)$ (24)

Then the Laplace-transform of the basic Equation (15) gives the general expression of the partial differential equation:

${\stackrel{¯}{\sigma }}_{ij,j}\left(X,s\right)-sc{\stackrel{¯}{u}}_{i}\left(X,s\right)-\rho {s}^{2}{\stackrel{¯}{u}}_{i}\left(X,s\right)=-{\stackrel{¯}{F}}_{i}\left(X,s\right)$ (25)

where

${\stackrel{¯}{F}}_{i}\left(X,s\right)=c{\stackrel{¯}{u}}_{i}\left(X,0\right)+\rho s{\stackrel{¯}{u}}_{i}\left(X,0\right)+\rho {\stackrel{¯}{\stackrel{˙}{u}}}_{i}\left(X,0\right)+{\stackrel{¯}{b}}_{i}\left(X,s\right)$ (26)

4. The MLPG Weak Formulation in Laplace-Transformed Domain

The MLPG (meshless local Petrov-Galerkin) method constructs the weak form over local subdomain such as ${\Omega }_{s}$ , which is a small region taken for each node inside the global domain. The local subdomains overlap, and cover the whole global domain $\Omega$ . In the present paper, the local subdomains are taken to be of a quadrature shape. The local weak form   of the governing Equation (25) can be written as:

${\int }_{{\Omega }_{s}}\left[{\stackrel{¯}{\sigma }}_{ij,j}\left(X,s\right)-\left(sc+\rho {s}^{2}\right){\stackrel{¯}{u}}_{i}\left(X,s\right)+{\stackrel{¯}{F}}_{i}\left(X,s\right)\right]{\Theta }_{i}\left(X\right)\text{d}\Omega =0,\text{\hspace{0.17em}}\text{\hspace{0.17em}}i\text{\hspace{0.17em}}\text{and}\text{\hspace{0.17em}}j=1,2$ (27)

where ${\Theta }_{i}\left(X\right)$ is a test function:

by using

${\stackrel{¯}{\sigma }}_{ij,j}{\Theta }_{i}={\left({\stackrel{¯}{\sigma }}_{ij}{\Theta }_{i}\right)}_{,j}-{\stackrel{¯}{\sigma }}_{ij}{\Theta }_{i,j}$ (28)

and applying the Gauss divergence theorem we can write:

$\begin{array}{l}\underset{\partial {\Omega }_{s}}{\int }{\stackrel{¯}{\sigma }}_{ij}\left(X,s\right){n}_{j}\left(X\right){\Theta }_{i}\left(X\right)\text{d}\Gamma -\underset{{\Omega }_{s}}{\int }{\stackrel{¯}{\sigma }}_{ij}\left(X,s\right){\Theta }_{i,j}\text{d}\Omega \\ -q\underset{{\Omega }_{s}}{\int }{\Theta }_{i}\left(X\right)\text{d}{\Omega }_{s}\text{ }{\stackrel{¯}{u}}_{i}^{h}\left(X,s\right)=-\underset{{\Omega }_{s}}{\int }{\stackrel{¯}{F}}_{i}\left(X,s\right){\Theta }_{i}\left(X\right)\text{d}\Omega \end{array}$ (29)

where

$q=sc+\rho {s}^{2}$ (30)

$\partial {\Omega }_{s}$ is the boundary of the local subdomain, which is consisted of three parts, i.e.

$\partial {\Omega }_{s}={\Gamma }_{si}\cup {\Gamma }_{st}\cup {\Gamma }_{su}$ (See Figure 1)

${\Gamma }_{si}$ is the local boundary that is totally inside the global domain,

${\Gamma }_{st}$ is the part of the local boundary, which lies on the global boundary with prescribed tractions, i.e. ${\Gamma }_{st}=\partial {\Omega }_{s}\cap {\Gamma }_{t}$

${\Gamma }_{su}$ is the part of the local boundary that lies on the global boundary with prescribed displacements, i.e. ${\Gamma }_{su}=\partial {\Omega }_{s}\cap {\Gamma }_{u}$

and considering:

${\stackrel{¯}{t}}_{i}\left(X,s\right)={\stackrel{¯}{\sigma }}_{ij}\left(X,s\right){n}_{j}\left(X\right)$ (31)

Figure 1. The support domain ΩS and integration domain Ωq for node I.

The local weak form in Equation (29) is leading to the following local integral equation:

$\begin{array}{l}-\underset{{\Gamma }_{si}}{\int }{\stackrel{¯}{t}}_{i}\left(X,s\right){\theta }_{i}\left(X\right)\text{d}\Gamma -\underset{{L}_{su}}{\int }{\stackrel{¯}{t}}_{i}\left(X,s\right){\theta }_{i}\left(X\right)\text{d}\Gamma +\underset{{\Omega }_{s}}{\int }{\stackrel{¯}{\sigma }}_{ij}\left(X,s\right){\Theta }_{i,j}\text{d}\Omega \\ +q\underset{{\Omega }_{s}}{\int }{\Theta }_{i}\left(X\right)\text{d}{\Omega }_{s}\text{ }{\stackrel{¯}{u}}^{h}{}_{i}\left(X,s\right)=\underset{{\Gamma }_{st}}{\int }{\stackrel{¯}{\stackrel{^}{t}}}_{i}\left(X,s\right){\theta }_{i}\left(X\right)\text{d}\Gamma +\underset{{\Omega }_{s}}{\int }{\stackrel{¯}{F}}_{i}\left(X,s\right){\Theta }_{i}\left(X\right)\text{d}\Omega \end{array}$ (32)

The strains can be obtained using the approximated displacements:

${\stackrel{¯}{\epsilon }}_{\left(3×1\right)}\left(X,s\right)={L}_{\left(3×2\right)}{\stackrel{¯}{u}}_{\left(2×1\right)}^{h}\left(X,s\right)$ (33)

Considering Equation (11) with Laplace transform:

${\stackrel{¯}{u}}^{h}\left(X,s\right)={\sum }_{i=1}^{n}{\phi }_{i}\left(X\right){\stackrel{¯}{u}}_{i}\left(s\right)={\Phi }^{\text{T}}\left(X\right)\stackrel{¯}{u}\left(s\right)$ (34)

${\stackrel{¯}{\epsilon }}_{\left(3×1\right)}\left(X,s\right)={L}_{\left(3×2\right)}{\sum }_{i=1}^{n}{\phi }_{i}\left(X\right){\stackrel{¯}{u}}_{i}\left(s\right)={L}_{\left(3×2\right)}{\Phi }_{\left(2×2n\right)}^{\text{T}}\left(X\right){\stackrel{¯}{u}}_{\left(2n×1\right)}\left(s\right)$ (35)

The constitutive equation gives the relationship between the stress and the strain:

${\stackrel{¯}{\sigma }}_{ij}\left(X,s\right)=D{\stackrel{¯}{\epsilon }}_{ij}\left(X,s\right)$

$\stackrel{¯}{\sigma }\left(X,s\right)=D{\sum }_{i=1}^{n}{B}_{i}\left(X\right){\stackrel{¯}{u}}_{i}\left(s\right)$ (36)

The traction vectors ${\stackrel{¯}{t}}_{i}\left(X,s\right)$ at a boundary point $X\in \partial {\Omega }_{s}$ are approximated by numbers of nodal values ${\stackrel{¯}{u}}_{i}\left(s\right)$ as:

${\stackrel{¯}{t}}_{i}\left(X,s\right)={\stackrel{¯}{\sigma }}_{ij}\left(X,s\right){n}_{j}\left( X \right)$

$\stackrel{¯}{t}\left(X,s\right)=N\left(X\right)D{\sum }_{i=1}^{n}{B}_{i}\left(X\right){\stackrel{¯}{u}}_{i}\left(s\right)$ (37)

where the matrix $N\left(X\right)$ related to the normal vector ${n}_{i}$ of $\partial {\Omega }_{s}$ by:

$N\left(X\right)=\left[\begin{array}{ccc}{n}_{1}& 0& {n}_{2}\\ 0& {n}_{2}& {n}_{1}\end{array}\right]$ (38)

The matrix ${B}_{i}\left(X\right)$ represented by the gradients of the shape functions as:

${B}_{i}\left(X\right)=\left[\begin{array}{cc}{\phi }_{i,1}& 0\\ 0& {\phi }_{i,2}\\ {\phi }_{i,2}& {\phi }_{i,1}\end{array}\right]$ (39)

The stress-strain matrix D for plane stress is deﬁned by:

$D=\frac{E}{1-{\nu }^{2}}\left[\begin{array}{ccc}1& \nu & 0\\ \nu & 1& 0\\ 0& 0& \frac{1-\nu }{2}\end{array}\right]$ (40)

In which E is the Young’s modulus and $\nu$ is the Poisson’s ratio.

Obeying the boundary conditions at those nodal points on the global boundary, where displacements are prescribed, and making use of the approximation formulae Equation (11), we obtains the discretized form of the displacement boundary conditions:

${\sum }_{k=1}^{n}{\phi }_{k}\left(\chi \right){\stackrel{¯}{u}}_{k}\left(s\right)=\stackrel{¯}{\stackrel{^}{u}}\left(\chi ,s\right)\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{for}\text{\hspace{0.17em}}\chi \in {\Gamma }_{u}$ (41)

Substituting equations. (36) and (37) into Equation (32) we obtain the discretized Local integral equations:

$\begin{array}{l}\left[-{\int }_{{\Gamma }_{si}}N\left(X\right)\theta \left(X\right)D{B}_{i}\left(X\right)\text{d}\Gamma -{\int }_{{\Gamma }_{su}}N\left(X\right)\theta \left(X\right)D{B}_{i}\left(X\right)\text{d}\Gamma \\ +{\int }_{{\Omega }_{s}}{W}_{i}D{B}_{i}\left(X\right)\text{d}{\Omega }_{s}+q{\int }_{{\Omega }_{s}}{\Phi }_{i}\left(X\right)\theta \left(X\right)\text{d}{\Omega }_{s}\right]\stackrel{¯}{u}\left(s\right)\\ =+{\int }_{{\Gamma }_{st}}{\stackrel{¯}{\stackrel{^}{t}}}_{i}\left(X,s\right)\theta \left(X\right)\text{d}\Gamma +{\int }_{{\Omega }_{s}}\stackrel{¯}{F}\left(X,s\right)\theta \left(X\right)\text{d}\Omega \end{array}$ (42)

where $\theta \left(X\right)$ is a matrix of weight functions given by:

$\theta \left(X\right)=\left[\begin{array}{cc}{\theta }_{i}\left(X\right)& 0\\ 0& {\theta }_{i}\left(X\right)\end{array}\right]$ (43)

${W}_{i}$ is a matrix that collects the derivatives of the weight:

${W}_{i}\left(X\right)=\left[\begin{array}{cc}{\theta }_{i,x}\left(X\right)& 0\\ 0& {\theta }_{i,y}\left(X\right)\\ {\theta }_{i,y}\left(X\right)& {\theta }_{i,x}\left(X\right)\end{array}\right]$ (44)

The cubic spline functions are used as the test functions for the local weak form:

${\theta }_{i}\left(r\right)=\left\{\begin{array}{l}\frac{2}{3}-4{r}_{j}^{2}+4{r}_{i}^{3}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{r}_{i}\le 0.5\\ \frac{4}{3}-4{r}_{i}+4{r}_{i}^{2}-\frac{4}{3}{r}_{i}^{3}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}0.5<{r}_{i}\le 1\\ 0\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{r}_{i}>1\end{array}$ (45)

where ${r}_{i}=\frac{|X-{X}_{i}|}{{d}_{i}}.$

In which $|X-{X}_{i}|$ is the distance from node ${X}_{i}$ to point $X$ , and ${r}_{i}$ is the size of the influence domain for the weight function.

Collecting the discretized local integral equations together with the discretized boundary conditions for displacements, we get the complete system of algebraic equations for computation of nodal displacements, which are the Laplace transforms of ﬁctitious parameters ${\stackrel{¯}{u}}_{k}\left(p\right)$ .

The time dependent values of the transformed variables can be obtained by an inverse Laplace transform. There are many inversion methods available for the Laplace transformation. In the present analysis, the Stehfest algorithm  is used. If $\stackrel{¯}{g}\left(p\right)$ is the Laplace-transform of $g\left(t\right)$ , an approximate value ${g}_{a}$ of ${g}_{k}\left(t\right)$ for a speciﬁc time is given by:

${g}_{a}\left(t\right)=\frac{\mathrm{ln}2}{t}{\sum }_{i=1}^{N}{v}_{i}\stackrel{¯}{g}\left(\frac{\mathrm{ln}2}{t}i\right)$ (46)

${v}_{i}={\left(-1\right)}^{N/2+i}\underset{k=\left[\left(i+1\right)/2\right]}{\overset{\mathrm{min}\left(i,N/2\right)}{\sum }}\frac{{k}^{N/2}\left(2k\right)!}{\left(N/2-k\right)!k!\left(k-1\right)!\left(i-k\right)!\left(2k-i\right)!}$

The selected number N = 10 with a single precision arithmetic is optimal to receive accurate results. It means that for each time t, it is needed to solve N boundary value problems for the corresponding Laplace parameters force $s=i\mathrm{ln}2/t$ with $i=1,2,\cdots ,N$ . If M denotes the number of the time instants in which we are interested to know g(t), the number of the Laplace transform solutions $\stackrel{¯}{g}\left({s}_{j}\right)$ is then M × N.

5. Numerical Results and Discussions

In this section, numerical results will be presented to illustrate the implementation and effectiveness of the proposed method. We present a numerical study for elastodynamic 2-D problem of a rectangular homogeneous isotropic plate  by using MLPG method, subjected to a dynamic force at the right end (Figure 2). A plane stress problem is considered, and a unit thickness is used. The dimensions of plate are length $L=48\text{\hspace{0.17em}}\text{m}$ and height $H=12\text{\hspace{0.17em}}\text{m}$ . The external excitation force $P=1.0f\left(t\right)$ , where $f\left(t\right)=\mathrm{sin}\left(wt\right)$ (simple harmonic load with $w=27\text{\hspace{0.17em}}\text{rad}/\text{s}$ ) is a function of time, the damping coefficient $c=0.4$ is fixed for all numerical computation. A total number ${n}_{t}=55$ uniformly distributed nodes is used, as shown in Figure 3, to represent the problem domain.

The dimension of the quadrature domain ${r}_{q}$ is set to ${d}_{q}={\alpha }_{q}{d}_{c}$ , where ${\alpha }_{q}$ is the size of the quadrature domain, ${d}_{c}$ is a distance to the ﬁrst nearest neighbouring point from node i. The dimension of the inﬂuence domain ${d}_{s}$ is set to ${d}_{s}={\alpha }_{s}{d}_{c}$ , where ${\alpha }_{s}$ is the size of the support domain.

The isotropic rectangular plate analyzed with the following material properties: $\left(E=200\text{\hspace{0.17em}}\text{GP};\nu =0.29;\rho =7860\text{\hspace{0.17em}}\text{kg}/{\text{m}}^{\text{3}}\right)$ in Steel. Some important parameters on the performance of the method have been investigated.

Figure 4 displays the variations of displacements ${u}_{y}$ as a function of time at point B under the harmonic load for different values of size of support domain ${\alpha }_{s}=1.5,1.8,1.9,2.0,2.5\text{\hspace{0.17em}}\text{and}\text{\hspace{0.17em}}3.0$ , where the time step is $\Delta t=0.01\text{\hspace{0.17em}}\text{s}$ . It can be seen

Figure 2. Rectangular homogeneous isotropic plate subjected to a dynamic force at the right end of the plate.

Figure 3.Conﬁguration and nodal arrangement for the plate.

Figure 4. Displacements uy at the middle point B at the free end of the plate excited by the time-step load for different values of ${\alpha }_{s}$ where $\Delta t=0.01\text{\hspace{0.17em}}\text{s}$ and ${\alpha }_{q}=2.0$ .

from this figure that the size of support domain inﬂuences on the results if ${\alpha }_{s}<1.8$ , and has a small effect on the results if ${\alpha }_{s}>1.8$ ). When ${\alpha }_{s}=3.0$ , the results obtained by the present method is very good compared with the other authors  that have used the Newmark method. In the following analysis,${\alpha }_{s}=3.0$ is employed.

The time variation of displacements ${u}_{y}$ are given in Figure 5 with different values of size of the quadrature domain ${\alpha }_{q}$ , where $\Delta t=0.01\text{\hspace{0.17em}}\text{s}$ . It is found that the size of the quadrature domain ${\alpha }_{q}$ inﬂuences seriously the results if ${\alpha }_{q}<0.8$ and has a small effect on the results if ${\alpha }_{q}>0.8$ .

Figure 5. The time variation of displacements uy for different values of ${\alpha }_{q}$ where $\Delta t=0.01\text{\hspace{0.17em}}\text{s}$ and ${\alpha }_{S}=3.0$ .

Figure 6. The time variation of displacements in the y direction at the point B with different time steps $\Delta t$ using the fixed values ${\alpha }_{s}=3.0,{\alpha }_{q}=2.0$ .

However, if the size of the quadrature domain is too large, the result obtained for displacements by MLPG is great. We found good results with ${\alpha }_{q}=2.0$ comparing with the results obtained by Long S. Y et al.  , they have used ${\alpha }_{q}=0.59$ .

Many time steps are used in computation to check the stability of the presented MLPG formulation. The displacement variations as a function of time and results for diﬀerent time steps are plotted in Figure 6. It can be found that when the time step $\Delta t$ is less than 0.02 s, perfect results have been obtained using the Laplace transform comparing with the results obtained by other authors   that have used the Newmark method. It also can be found that when a time step $\Delta t$ is larger than 0.02 s, the results are not convergent and not accurate.

6. Conclusion

The present method MLPG that uses the cubic spline test function is used to analyze elastodynamic problem. The equation formulation based on MLPG method in Laplace transform and time domain with MLS approximation has been successfully implemented to solve elastodynamic problems in isotropic solids, subjected to a dynamic force at the right end of the plate. We found that the amplitude of the vibration decreases with time because the effects of damping and the harmonic load, the response should converge to the static deformation. We found that when the time step $\Delta t$ is less than 0.02 s, perfect results have been obtained by using the Laplace transform and when a time step $\Delta t$ is larger than 0.02 s the results are not accurate. We found that the size of the quadrature domain ${\alpha }_{q}$ has a small effect on the results if ${\alpha }_{q}>0.8$ , and the size of the support domain ${\alpha }_{s}$ influences on the results if ${\alpha }_{s}<1.8$ . The sizing parameters ${\alpha }_{s}$ and ${\alpha }_{q}$ , which decides the size of the subdomain needs to be chosen carefully, especially, in the dynamic analysis.

Cite this paper

Eddaoudy, J. and Bouziane, T. (2018) A Study of the Elastodynamic Problem by Meshless Local Petrov-Galerkin Method Using the Laplace- Transform. World Journal of Mechanics, 8, 46-58. https://doi.org/10.4236/wjm.2018.82004

References

1. 1. Atluri, S.N. and Zhu, T.L. (1998) A New Meshless Local Petrov-Galerkin (MLPG) Approach in Computational Mechanics. Computational Mechanics, 22, 117-127. https://doi.org/10.1007/s004660050346

2. 2. Atluri, S.N. and Zhu, T.L. (2000) The Meshless Local Petrov-Galerkin (MLPG) Approach for Solving Problems in Elasto-Statics. Computational Mechanics, 25, 169-179. https://doi.org/10.1007/s004660050467

3. 3. Sladek, J., Sladek, V., Wen, P.H. and Aliabadi, M.H. (2006) Meshless Local Petrov-Galerkin (MLPG) Method for Shear Deformable Shells Analysis. CMES: Computer Modeling in Engn. & Sciences, 13, 103-118.

4. 4. Sladek, J., Sladek, V., Zhang, C. and Solek, P. (2007) Application of the MLPG to Thermopiezoelectricity. CMES: Computer Modeling in Engn. & Sciences, 22, 217-233.

5. 5. Mirzaei, D. and Dehghan, M. (2010) Meshless Local Petrov-Galerkin (MLPG) Approximation to the Two Dimensional Sine-Gordon Equation. Journal of Computational and Applied Mathematics, 233, 2737-2754. https://doi.org/10.1016/j.cam.2009.11.022

6. 6. Sladek, J., Stanak, P., Han, Z.D., Sladek, V. and Atluri, S.N. (2013) Applications of the MLPG Method in Engineering & Sciences: A Review. CMES, Computer Modeling in Engn. & Sciences, 92, 423-475.

7. 7. Zhang, T., Dong, L.T., Alotaibi, A. and Satya, A.N. (2013) Application of the MLPG Mixed Collocation Method for Solving Inverse Problems of Linear Isotropic/Anisotropic Elasticity with Simply/Multiply-Connected Domains. CMES, Computer Modeling in Engn. & Sciences, 94, 1-28.

8. 8. Stanak, P., Sladek, V., Sladek, J., Krahuleca, S. and Sator, L. (2013) Application of Patch Test in Meshless Analysis of Continuously Non-Homogeneous Piezoelectric Circular Plate. Applied and Computational Mechanics, 7, 65-76.

9. 9. Atluri, S.N. and Shen, S.P. (2002) The Meshless Local Petrov-Galerkin (MLPG) Method: A Simple & Less-Costly Alternative to the Finite Element and Boundary Element Methods. CMES, 3, 11-51.

10. 10. Sladek, J., Sladek, V. and Keer, R.V. (2003) Meshless Local Boundary Integral Equation Method for 2D Elastodynamic Problems. International Journal for Numerical Methods in Engineering, 57, 235-249. https://doi.org/10.1002/nme.675

11. 11. Qian, L.-F. and Ching, H.-K. (2004) Static and Dynamic Analysis of 2-d Functionally Graded Elasticity by using Meshless Local Petrov Galerkin Method. Journal of the Chinese Institute of Engineers, 27, 491-503. https://doi.org/10.1080/02533839.2004.9670899

12. 12. Han, Z.D. and Atluri, S.N. (2004) Meshless Local Petrov-Galerkin (MLPG) Approaches for Solving 3D Problems in Elasto-Statics. Computer Modeling in Engineering & Sciences, 6, 168-188.

13. 13. Han, Z.D. and Atluri, S.N. (2004) A Meshless Local Petrov-Galerkin (MLPG) Approach for 3-Dimensional Elastodynamics. Computers Materials & Continua, 1, 129-140.

14. 14. Atluri, S.N., Han, Z.D. and Shen, S. (2003) Meshless Local Petrov-Galerkin (MLPG) Approaches for Solving the Weakly-Singular Traction & Displacement Boundary Integral Equations, CMES, 4, 507-517.

15. 15. Lancaster, P. and Salkauskas, K. (1981) Surfaces Generated by Moving Least Squares Methods. Mathematics of Computation 37, 141-158. https://doi.org/10.1090/S0025-5718-1981-0616367-1

16. 16. Li, D., Wei, A., Luo, K. and Fan, J. (2015) An Improved Moving-Least-Squares Reconstruction for Immersed Boundary Method, International Journal for Numerical Methods in Engineering, 104, 789-804. https://doi.org/10.1002/nme.4949

17. 17. Shivanian, E. (2016) Local Integration of Population Dynamics via Moving Least Squares Approximation, Engineering with Computers, 32, 331-342. https://doi.org/10.1007/s00366-015-0424-z

18. 18. Belytschko, T., Krongauz, Y., Organ, D., Fleming, M. and Krysl, P. (1996) Meshless Methods: An Overview and Recent Developments. Computer Methods in Applied Mechanics and Engineering, 139, 3-47. https://doi.org/10.1016/S0045-7825(96)01078-X

19. 19. Pandey, S.S., Kasundra, P.K. and Daxini, S.D. (2013) Introduction of Meshfree Methods and Implementation of Element Free Galerkin (EFG) Method to Beam Problem. International Journal on Theoretical and Applied Research in Mechanical Engineering, 2, 2319-3182.

20. 20. Atluri, S.N., Cho, J.Y. and Kim, H.G. (1999) Analysis of Thin Beams using the Meshless Local Petrov-Galerkin Method with Generalized Moving Least Squares Interpolations. Computational Mechanics, 24, 334-347. https://doi.org/10.1007/s004660050456

21. 21. Gu, Y.T. and Liu, G.R. (2001) A Meshless Local Petrov-Galerkin (MLPG) Method for Free and Forced Vibration Analyses for Solids. Computational Mechanics, 27, 188-198. https://doi.org/10.1007/s004660100237

22. 22. Long, S.Y., Liu, K.Y. and Hu, D.A. (2006) A New Meshless Method Based on MLPG for Elastic Dynamic Problems. Engineering Analysis with Boundary Elements, 30, 43-48. https://doi.org/10.1016/j.enganabound.2005.09.001

23. 23. Xia, P., Long, S. and Cui, H. (2009) Elastic Dynamic Analysis of Moderately Thick Plate using Meshless LRPIM. Acta Mechanica Solida Sinica, 22, 116-124. https://doi.org/10.1016/S0894-9166(09)60096-3

24. 24. Hu, D., Yi, G.W., Li, Y., Han, X. and Gu, Y.T. (2014) A Meshfree-Based Local Galerkin Method with Condensation of Degree of Freedom for Elastic Dynamic Analysis. Acta Mechanica Sinica, 30, 92-99. https://doi.org/10.1007/s10409-013-0090-6

25. 25. Moussaoui, A. and Bouziane, T. (2013) Comparative Study of the Effect of the Parameters of Sizing Data on Results by the Meshless Methods (MLPG). World Journal of Mechanics, 3, 82-87. https://doi.org/10.4236/wjm.2013.31006

26. 26. Timoshenk, S. (1959) Theory of Plates and Shells. McGraw-Hill, New York, 2.

27. 27. Davies, B. and Martin, B. (1979) Numerical Inversion of the Laplace Transform: A Survey and Comparison of Methods. Journal of Computational Physics, 33, 1-32. https://doi.org/10.1016/0021-9991(79)90025-1

28. 28. Stehfest (1970) Algorithm 368: Numerical Inversion of Laplace Transform. Communications of the ACM, 13, 47-49.