Two-dimensional beams in rectangular coordinates using the radial point interpolation method

Abstract

The three-dimensional Theory of Elasticity equations lead to a complex solution for most problems in engineering. Therefore, the solutions are typically developed for reduced systems, usually symmetrical or two-dimensional. In this context, computational resources allow the reduction of these simplifications. The most recognized methods of algebraic approximation of the differential equations are the Finite Differences Method and the Finite Element Method (FEM). However, they have limitations in mesh generation and/or adaptation. As follows, Meshless Methods appear as an alternative to these options. The present work uses the Radial Point Interpolation Method (RPIM) to evaluate the stress in two-dimensional beams in regions close to loading (Saint Venant's Principle). Formulations based on the Fourier Series Theory and the RPIM are presented. Multiquadrics Radial Basis Functions were used to obtain the stiffness matrix. Two numerical examples demonstrate the validity of the RPIM for the proposed theme. The results were obtained from the formulations cited and the Finite Element Method for comparison.

Keywords:
two-dimensional beams; Saint-Venant's principle; Radial Point Interpolation Method; stress analysis

1. Introduction

The analytical solution to most problems of the Theory of Elasticity is difficult due to the complexity of the equations. Therefore, the resolutions are typically designed for reduced systems, usually symmetrical or two-dimensional (Saad, 2005).

In this way, the computational analysis has considerable relevance in the solution of these problems. The most known methods of numerical analysis are the Finite Differences Method and the Finite Element Method (FEM), the latter being the most used. However, it includes limitations, mainly in mesh generation and adaptation. In this manner, Meshless Methods appear as a significant alternative to these options (Liu, 2010).

In the previous two decades, Meshless Methods have been used in several engineering areas. Silva (2012) explored the application of the Element Free Galerkin (EFG) Method in physically non-linear static structures of reinforced concrete. Asprone et al. (2014) investigate the Modified Finite Particle Method (MFPM) and propose modifications to it in the static and dynamic problems, both in the elastic range. Hu et al. (2014) developed a technique to condense the degrees of freedom to increase the computational efficiency of the meshless methods in dynamic linear elastic analysis. The equations of the Plane Theory of Elasticity can be applied to two cases of practical interest: plane stress and strain of thin plates under forces applied to their boundaries and acting in their planes. An important fact to be observed in the structure is the effect of loading in regions close to the point of application.

This effect is called the Saint-Venant's Principle. It enunciates that two statically equivalent force systems acting over a small portion Ps of the surface of a body produce (approximately) the same stress and displacement at a point sufficiently far from Ps in the body where the force systems act.

Relevant researches have been published about Saint-Venant's Theory. Genoese et al. (2014) examined a geometrically nonlinear model for homogeneous and isotropic beams including non-uniform warping due to torsion and shear derived from the Saint-Venant's rod. Genoese et al. (2013) also presented an alternative linear model for thin-walled section beams, whose formulation is based on the Hellinger-Reissner Principle. Zhao et al. (2012) proposed an approach to investigate the Saint-Venant's problem in graded beams with Young's Modulus varying exponentially in the axial direction and constant Poisson Ratio. Fatmi and Ghazouani (2011) suggested a higher-order composite beam theory, which can be viewed as an extension of the Saint-Venant's Theory. Petrolo and Casciaro (2004) investigated the use of the Saint-Venant's general rod theory for deriving the stiffness matrix in three-dimensional beam elements with a general cross-section.

The proposed research aims to demonstrate the Saint-Venant Principle for two-dimensional beams using the Radial Point Interpolation Method (RPIM). The formulations for RPIM and the analytical solution provided by the Fourier Series are presented. Two examples are demonstrated to validate the RPIM. The results are compared with the analytical solution and numerical solution of the Finite Element Method utilizing the SAP2000^(c) software.

2. Two-dimensional beams in rectangular coordinates

2.1 Solution based on fourier series theory

The biharmonic equation for the stress functions in two-dimensional problems is given by:

(1) ${\nabla }^4{\varphi }^{\prime }=0\Rightarrow \frac{{\partial }^4{\varphi }^{\prime }}{\partial {\mathit{x}}^4}+2\frac{{\partial }^4{\varphi }^{\prime }}{\partial {\mathit{x}}^2\partial {\mathit{y}}^2}+\frac{{\partial }^4{\varphi }^{\prime }}{\partial {\mathit{y}}^4}=0$

where f'=f'(x,y) is the Airy Stress Function. A general solution may be found by Separation of Variables with Fourier Series (Saad, 2005). In cartesian coordinates:

(2) ${\varphi }^{\prime }\left(\mathit{x},\mathit{y}\right)=\mathit{X}\left(\mathit{x}\right)\mathit{Y}\left(\mathit{y}\right)$

In Eq. (2), X(x)=e ^{a x} and Y(y)=e ^{b y} . Replacing in Eq. (1):

(3) $\left({\mathit{\alpha }}^4+2{\mathit{\alpha }}^2{\mathit{\beta }}^2+{\mathit{\beta }}^4\right){\mathit{e}}^{\mathit{\alpha }\mathit{x}}{\mathit{e}}^{\mathit{\beta }\mathit{y}}=0$

The term in parentheses must be zero, leading to the following characteristic equation:

(4) ${\left({\mathit{\alpha }}^2+{\mathit{\beta }}^2\right)}^2=0$

(5) $\mathit{\alpha }=\pm \mathit{i}\mathit{\beta }$

The general solution includes zero root and general roots. In the case of zero root with b=0, there are 3 additional roots (Eq. 6). For the case with a=0, the solution is given by Eq. (7):

(6) ${{\varphi }^{\prime }}_{\mathit{\beta }=0}={\mathit{C}}_0+{\mathit{C}}_1\mathit{x}+{\mathit{C}}_2{\mathit{x}}^2+{\mathit{C}}_3{\mathit{x}}^3$

(7) ${{\varphi }^{\prime }}_{\mathit{\alpha }=0}={\mathit{C}}_4\mathit{y}+{\mathit{C}}_5{\mathit{y}}^2+{\mathit{C}}_6{\mathit{y}}^3+{\mathit{C}}_7\mathit{xy}+{\mathit{C}}_8{\mathit{x}}^2\mathit{y}+{\mathit{C}}_9{\mathit{xy}}^2$

The Eqs. (6) and (7) satisfy the Eq. (1). Therefore:

(8) ${\varphi }^{\prime }\left(\mathit{x},\mathit{y}\right)={\mathit{e}}^{\mathit{i}\mathit{\beta }\mathit{x}}\left[{\mathit{A}}_1{\mathit{e}}^{\mathit{\beta }\mathit{y}}+{\mathit{A}}_2{\mathit{e}}^{-\mathit{\beta }\mathit{y}}+{\mathit{A}}_3{\mathit{ye}}^{\mathit{\beta }\mathit{y}}+{\mathit{A}}_4{\mathit{ye}}^{-\mathit{\beta }\mathit{y}}\right]+{\mathit{e}}^{-\mathit{i}\mathit{\beta }\mathit{x}}\left[{{\mathit{A}}^{\prime }}_1{\mathit{e}}^{\mathit{\beta }\mathit{y}}-{{\mathit{A}}^{\prime }}_2{\mathit{e}}^{\mathit{\beta }\mathit{y}}+{{\mathit{A}}^{\prime }}_3{\mathit{ye}}^{\mathit{\beta }\mathit{y}}+{{\mathit{A}}^{\prime }}_4{\mathit{ye}}^{-\mathit{\beta }\mathit{y}}\right]$

where Ci, Ai, and A'i are arbitrary constants determined by boundary conditions. The complete solution is given by the superposition of Eqs. (6), (7) and (8). Substituting exponentials for equivalent trigonometric and hyperbolic forms:

(9) $\begin{array}{c}{\varphi }^{\prime }\left(\mathit{x},\mathit{y}\right)=\sin \mathit{\beta }\mathit{x}\left[\left({\mathit{A}}_1+{\mathit{A}}_3\mathit{\beta }\mathit{y}\right)\sin h\mathit{\beta }\mathit{y}+\left({\mathit{A}}_2+{\mathit{A}}_4\mathit{\beta }\mathit{y}\right)\cosh \mathit{\beta }\mathit{y}\right]+\hfill \\ \cos \mathit{\beta }\mathit{x}\left[\left({{\mathit{A}}^{\prime }}_1+{{\mathit{A}}^{\prime }}_3\mathit{\beta }\mathit{y}\right)\sin h\mathit{\beta }\mathit{y}+\left({{\mathit{A}}^{\prime }}_2+{{\mathit{A}}^{\prime }}_4\mathit{\beta }\mathit{y}\right)\cosh \mathit{\beta }\mathit{y}\right]+\hfill \\ \sin \mathit{\alpha }\mathit{y}\left[\left({\mathit{A}}_5+{\mathit{A}}_7\mathit{\alpha }\mathit{x}\right)\sinh \mathit{\alpha }\mathit{x}+\left({\mathit{A}}_6+{\mathit{A}}_8\mathit{\alpha }\mathit{x}\right)\cosh \mathit{\alpha }\mathit{x}\right]+\hfill \\ \cos \mathit{\alpha }\mathit{y}\left[\left({{\mathit{A}}^{\prime }}_5+{{\mathit{A}}^{\prime }}_7\mathit{\alpha }\mathit{x}\right)\sinh \mathit{\alpha }\mathit{x}+\left({{\mathit{A}}^{\prime }}_6+{{\mathit{A}}^{\prime }}_8\mathit{\alpha }\mathit{x}\right)\cosh \mathit{\alpha }\mathit{x}\right]+{{\varphi }^{\prime }}_{\mathit{\alpha }=0}+{{\varphi }^{\prime }}_{\mathit{\beta }=0}\hfill \end{array}$

The stresses can then be obtained from differential relations:

(10) ${\mathit{\sigma }}_{\mathit{x}}=\frac{{\partial }^2{\varphi }^{\prime }\left(\mathit{x},\mathit{y}\right)}{{\partial }^2\mathit{y}}$

(11) ${\mathit{\sigma }}_{\mathit{y}}=\frac{{\partial }^2{\varphi }^{\prime }\left(\mathit{x},\mathit{y}\right)}{{\partial }^2\mathit{x}}$

(12) ${\mathit{\tau }}_{\mathit{xy}}=\frac{{\partial }^2{\varphi }^{\prime }\left(\mathit{x},\mathit{y}\right)}{\partial \mathit{x}\partial \mathit{y}}$

The applications of the Fourier solution method usually incorporate the Fourier series theory (Saad, 2005). A periodic function f(x) with period 2L can be represented on the interval (-L,L) by the Fourier trigonometric series:

(13) $\mathit{f}\left(\mathit{x}\right)=\frac{1}{2}{\mathit{a}}_0+\sum _{\mathit{n}=1}^{\infty }\left({\mathit{a}}_{\mathit{n}}\cos \frac{\mathit{n}\mathit{\pi }\mathit{x}}{\mathit{L}}+{\mathit{b}}_{\mathit{n}}\sin \frac{\mathit{n}\mathit{\pi }\mathit{x}}{\mathit{L}}\right)$

(14) ${\mathit{a}}_{\mathit{n}}=\frac{1}{\mathit{L}}\underset{-\mathit{L}}{\overset{\mathit{L}}{\int }}\mathit{f}\left(\mathit{\xi }\right)\cos \frac{\mathit{n}\mathit{\pi }\mathit{\xi }}{\mathit{L}}\mathit{d}\mathit{\xi }\mathit{n}=\mathit{0},\mathit{1},\mathit{2},\dots$

(15) ${\mathit{b}}_{\mathit{n}}=\frac{1}{\mathit{L}}\underset{-\mathit{L}}{\overset{\mathit{L}}{\int }}\mathit{f}\left(\mathit{\xi }\right)\sin \frac{\mathit{n}\mathit{\pi }\mathit{\xi }}{\mathit{L}}\mathit{d}\mathit{\xi }\mathit{n}=\mathit{1},\mathit{2},\mathit{3},\dots$

These expressions can be simplified in some cases. If f(x) is an even function, f(x)=-f(x) and Eq. (13) reduces to the Fourier cosine series (Eqs. 16 and 17). If f(x) is an odd function, f(x)=-f(-x) and Eq. (13) reduces to the Fourier sine series (Eqs. 18 and 19):

(16) $\mathit{f}\left(\mathit{x}\right)=\frac{1}{2}{\mathit{a}}_0+\sum _{\mathit{n}=1}^{\infty }{\mathit{a}}_{\mathit{n}}\cos \frac{\mathit{n}\mathit{\pi }\mathit{x}}{\mathit{L}}$

(17) ${\mathit{a}}_{\mathit{n}}=\frac{2}{\mathit{L}}\underset{0}{\overset{\mathit{L}}{\int }}\mathit{f}\left(\mathit{\xi }\right)\cos \frac{\mathit{n}\mathit{\pi }\mathit{\xi }}{\mathit{L}}\mathit{d}\mathit{\xi }\mathit{n}=\mathit{0},\mathit{1},\mathit{2},\dots$

(18) $\mathit{f}\left(\mathit{x}\right)=\sum _{\mathit{n}=1}^{\infty }{\mathit{b}}_{\mathit{n}}\sin \frac{\mathit{n}\mathit{\pi }\mathit{x}}{\mathit{L}}$

(19) ${\mathit{b}}_{\mathit{n}}=\frac{2}{\mathit{L}}\underset{0}{\overset{\mathit{L}}{\int }}\mathit{f}\left(\mathit{\xi }\right)\sin \frac{\mathit{n}\mathit{\pi }\mathit{\xi }}{\mathit{L}}\mathit{d}\mathit{\xi }\mathit{n}=\mathit{1},\mathit{2},\mathit{3},\dots$

2.2 Solution based on Radial Point Interpolation Method (RPIM)

The use of polynomials to create basis functions is advantageous for two reasons: simplicity and good numerical precision. Besides, shape functions of any order can be reproduced by increasing the number of interpolation points (field nodes). Among these advantages, the RPIM method of obtaining form functions avoids the occurrence of singularities in the moment matrix (Liu and Gu, 2005). The displacement approximation u^h at a point of interest x^T={x,y} is given by (Liu, 2010):

(20) $\mathit{u}\left(\mathbf{x}\right)={\sum _{\mathit{i}=1}^{\mathit{n}}{\mathit{R}}_{\mathit{i}}\left(\mathbf{x}\right){\mathit{a}}_{\mathit{i}}=\mathbf{R}}^{⊺}\left(\mathbf{x}\right)\mathbf{a}$

(21) ${\mathbf{a}}^{⊺}=\left\{{\mathit{a}}_1{\mathit{a}}_2{\mathit{a}}_3\dots {\mathit{a}}_{\mathit{n}}\right\}$

where R_i is the Radial Basis Function (RBF), a is a vector of unknown constants and n is the number of nodes in a support domain. The distance r between points x and x_i is obtained by:

(22) $\mathit{r}=\sqrt{{\left(\mathit{x}-{\mathit{x}}_{\mathit{i}}\right)}^2+{\left(\mathit{y}-{\mathit{y}}_{\mathit{i}}\right)}^2}$

The vector of Radial Basis Functions R has the following form:

(23) ${\mathbf{R}}^{\mathit{⊺}}\left(\mathbf{x}\right)=\left\{{\mathit{R}}_1\left(\mathbf{x}\right){\mathit{R}}_2\left(\mathbf{x}\right){\mathit{R}}_3\left(\mathbf{x}\right)\dots {\mathit{R}}_{\mathit{n}}\left(\mathbf{x}\right)\right\}$

Table 1 presents the four most often used forms of radial functions R_i (x). The parameters can be tuned for better performance.

Table 1
Radial Basis Functions and dimensionless shape parameters.

Function Type	Expression	Shape Parameter
Multiquadrics (MQ)	Ri (x,y)=(ri2+(ac dc )2 )q	ac≥ 0, q
Gaussian (EXP)	Ri (x,y) = exp (- cri2)= exp {-c[(x-xi)2+(y-yi)2]}	C
Thin plate spline (TPS)	Ri (x,y) = rih = [(x-xi)2+(y-yi )2 ]h	H
Logarithmic RBF	Ri (ri)=rih log(ri)	H

Misra and Kumar (2013) point out that Multiquadrics radial basis functions (MQ-RBF) present advantages, such as easy implementation for structural analysis and reasonable results for a small number of field nodes. Besides, its implementation is highly suitable, and no connectivity is required for arbitrarily distributed nodes. The main idea in the MQ method is to create a coefficient matrix with a significant number of zero elements for reducing the computational costs (Fallah et al., 2019). Thus, MQ-RBF was used in the present study. In Table 1, a_c is the dimensionless shape parameter, d_c is the characteristic length (usually the average nodal spacing for all the n nodes in the support domain) and q is an exponent parameter.

The interpolation at the point k has the form:

(24) ${\mathit{u}}_{\mathit{k}}=\mathit{u}\left({\mathit{x}}_{\mathit{k}},{\mathit{y}}_{\mathit{k}}\right)=\sum _{\mathit{i}=1}^{\mathit{n}}{\mathit{a}}_{\mathit{i}}{\mathit{R}}_{\mathit{i}}\left({\mathit{x}}_{\mathit{k}},{\mathit{y}}_{\mathit{k}}\right)\mathit{k}=1,2,\dots ,\mathit{n}$

In matrix form, these n equations can be written as:

(25) ${\mathbf{d}}_{\mathit{S}}={\mathbf{R}}_{\mathit{Q}}\mathbf{a}$

In the Eq. (25), d_s is the vector within the field nodal variables at the n local nodes and R_Q is the moment matrix of Radial Basis Functions:

(26) ${\mathbf{R}}_{\mathit{Q}}=\left[\begin{array}{cccc}{\mathit{R}}_1\left({\mathit{r}}_1\right)& {\mathit{R}}_2\left({\mathit{r}}_1\right)& \cdots & {\mathit{R}}_{\mathit{n}}\left({\mathit{r}}_1\right)\\ {\mathit{R}}_1\left({\mathit{r}}_2\right)& {\mathit{R}}_2\left({\mathit{r}}_2\right)& \cdots & {\mathit{R}}_{\mathit{n}}\left({\mathit{r}}_2\right)\\ ⋮& ⋮& \ddots & ⋮\\ {\mathit{R}}_1\left({\mathit{r}}_{\mathit{n}}\right)& {\mathit{R}}_2\left({\mathit{r}}_{\mathit{n}}\right)& \cdots & {\mathit{R}}_{\mathit{n}}\left({\mathit{r}}_{\mathit{n}}\right)\end{array}\right]$

(27) ${\mathit{r}}_{\mathit{k}}=\sqrt{{\left({\mathit{x}}_{\mathit{k}}-{\mathit{x}}_{\mathit{i}}\right)}^2-{\left({\mathit{y}}_{\mathit{k}}-{\mathit{y}}_{\mathit{i}}\right)}^2}$

Since the distance has no direction, then:

(28) ${\mathit{R}}_{\mathit{i}}\left({\mathit{r}}_{\mathit{j}}\right)={\mathit{R}}_{\mathit{j}}\left({\mathit{r}}_{\mathit{i}}\right)$

which indicates symmetry of the matrix R_Q .A unique solution for a is then obtained by:

(29) $\mathbf{a}={\mathbf{R}}_{\mathit{Q}!}^{-1}{\mathbf{d}}_{\mathit{s}}$

Replacing Eq. (29) into (28):

(30) $\mathit{u}\left(\mathbf{x}\right)={\mathbf{R}}^{\mathit{⊺}}\left(\mathbf{x}\right){\mathbf{R}}_{\mathit{Q}}^{-1}{\mathbf{d}}_{\mathit{s}}=\mathbf{\Phi }\left(\mathbf{x}\right){\mathbf{d}}_{\mathit{s}}$

where Φ(x) is the vector of shape functions:

(31) $\mathbf{\Phi }\left(\mathbf{x}\right)={\mathbf{R}}^{\mathit{⊺}}\left(x\right){\mathbf{R}}_{\mathit{Q}}^{-1}=\left\{{\mathit{R}}_1\left(\mathbf{x}\right){\mathit{R}}_2\left(\mathbf{x}\right){\mathit{R}}_3\left(\mathbf{x}\right)\dots {\mathit{R}}_{\mathit{n}}\left(\mathbf{x}\right)\right\}{\mathbf{R}}_{\mathit{Q}}^{-1}=\left\{{\varphi }_1\left(\mathbf{x}\right){\varphi }_2\left(\mathbf{x}\right){\varphi }_3\left(\mathbf{x}\right)\cdots {\mathit{\phi }}_{\mathit{n}}\left(\mathbf{x}\right)\right\}$

and f_k is the shape function for the nodek:

(32) ${\varphi }_{\mathit{k}}\left(\mathbf{x}\right)=\sum _{\mathit{i}=1}^{\mathit{n}}{\mathit{R}}_{\mathit{i}}\left(\mathbf{x}\right){\mathit{S}}_{\mathit{ik}}^{\mathit{a}}$

In the Eq. (33), S^a _ik is the (i,k) element of the constant matrix R^-1 _Q in the support domain. The equilibrium equation for the problem can be put in matrix form as:

(33) $\mathbf{Ku}=\mathbf{F}$

(34) $\mathbf{F}={\mathbf{F}}^{\mathit{b}}+{\mathbf{F}}^{\mathit{t}}=\underset{\Omega }{\int }{\mathbf{\Phi }}^{\mathbf{⊺}}\mathbf{b}\mathit{d}\Omega +\underset{\Gamma \mathit{t}}{\int }{\mathbf{\Phi }}^{\mathit{⊺}-}\mathbf{t}\mathit{d}\Gamma$

In the Eqs. (33) and (34), u is the displacements of field nodes, F is the global vector of forces, F^b is the global body force vector at the domain Ω, F^t is the global traction force vector at boundary domain G, b is the body force vector and t is the external traction force vector.

The global stiffness matrix K is defined as:

(35) $\mathbf{K}={\sum _{\mathit{I}}^{\mathit{N}}\sum _{\mathit{J}}^{\mathit{N}}\mathbf{K}}_{\mathit{IJ}}$

(36) ${\left({\mathbf{K}}_{\mathit{IJ}}\right)}_{2\mathit{x}2}=\underset{\Omega }{\int }{\left({\mathbf{B}}_{\mathit{I}}\right)}_{2\mathit{x}3}^{\mathit{T}}{\mathbf{D}}_{2\mathit{x}3}{\left({\mathbf{B}}_{\mathit{J}}\right)}_{3\mathit{x}2}\mathit{d}\Omega$

(37) ${\mathbf{B}}_{3\mathit{x}2\mathit{n}}=\left[\begin{array}{ccccc}\frac{\partial {\varphi }_1}{\partial \mathit{x}}& 0& \dots & \frac{\partial {\varphi }_{\mathit{n}}}{\partial \mathit{x}}& 0\\ 0& \frac{\partial {\varphi }_1}{\partial \mathit{y}}& \dots & 0& \frac{\partial {\varphi }_{\mathit{n}}}{\partial \mathit{y}}\\ \frac{\partial {\varphi }_1}{\partial \mathit{y}}& \frac{\partial {\varphi }_1}{\partial \mathit{x}}& \dots & \frac{\partial {\varphi }_{\mathit{n}}}{\partial \mathit{y}}& \frac{\partial {\varphi }_{\mathit{n}}}{\partial \mathit{x}}\end{array}\right]$

(38) ${\mathbf{D}}_{3\mathit{x}3}=\frac{\mathit{E}}{1-{\mathit{v}}^2}\left[\begin{array}{ccc}1& \mathit{v}& 0\\ \mathit{v}& 1& 0\\ 0& 0& \frac{\left(1-\mathit{v}\right)}{2}\end{array}\right]$

where K_IJ is the nodal stiffness matrix, B the strain matrix and D the matrix of elastic constants.

Liu and Gu (2005) demonstrate that the interpolation quality changes with the exponent q. However, the RPIM-MQ fails because of the singularity of the moment matrix for q=1.0, 2.0 and 3.0. According to the authors, the preferred value of parameter q is close to 1.0 or 2.0 (0.98, 1.03 or 1.99 being recommended). The same authors observed that the a_c shape parameter has less influence than q(a_c≥1.0 is recommended). Besides this, the average fitting errors of function values over the entire domain decreases when the number of interpolation points in the entire domain (N) increases.

The RPIM code was written in FORTRAN language and divided into modules to make the management of the main program easier.

3. Examples

In the present study,SAP2000^(c) was used to obtain the solution by the Finite Element Method. Shell elements were used, and the number of elementsin each example was chosen so that their nodes matched the positions of the RPIM field nodes. The other data were the same as described in the examples.

3.1 Example 1: Beam under equidistant forces P

The first example shows a beam subjected to two equal forces of P=1.2N (Fig. 1a) distant b from the middle section s-s (TIMOSHENKO and GOODIER, 1980). The beam has a height H=1.2m (c=0.6m), length L=4.8m (l=2.4m) and base B=1m. The YoungModulus is 200GPa and Poisson's ratio 0.3. The number of field nodes is 891 to represent the domain (Fig. 1b) and 800 background cells for integrations with 2 Gauss points in each one. The parameters for the radial shape functions are a_c=1.0, d_c=2.0 and q=1.03.

Figure 1
Beam subjected to two equal forces P: (a) geometry; (b) model discretized in field nodes.

The analytical solution was obtained by Fourier Series Theory. Fig. 2a presents good agreement between the present study and the analytical result in the s-s cross-section. More significant variations can be observed in y=±1.2m. The stress at centroid (y=0.6m) is the same for all curves. In Fig. 2b, it can be seen that the shear stress at the ends obtained using FEM is closer to the analytical response. However, the region between y=0.2m and y=1.0m is better described by RPIM, with asmall difference about Timoshenko and Goodier (1980).The shear stress curve for the RPIM (Fig. 2c) shows consistency compared with Timoshenko (1980). The FEM demonstrates a small variation between y=0.36m to y=0.72m.Fig. 2d shows that the numerical responses differ from the analytical response for shear stress.

Figure 2
Numerical and analytical shear stress for Example 1: (a) b/c=1/10; (b) b/c=1/5; (c) b/c=1/2; (d) b/c=∞.

3.2 Example 2: Cantilevered beam under axial and transverse load

The second example refers to a cantilevered beam under forces N=1.8N and P=1.2N (Fig. 3a). This example is proposed by Saad (2005, p.192). The Airy Stress Function presented as an analytical solution of the problem (formulated in terms of the resulting force system) is given as:

(39) ${\varphi }^{\prime }\left(\mathit{x},\mathit{y}\right)=\frac{3\mathit{P}}{4\mathit{c}}\left(\mathit{xy}-\frac{{\mathit{xy}}^3}{3{\mathit{c}}^2}\right)+\frac{\mathit{N}}{4{\mathit{c}}^2}{\mathit{y}}^2$

For the beam, 275 field nodes were used to represent the domain (Fig. 3b) and 240 background cells for integrations, with 4 Gauss points in each one. The parameters for the radial shape functions are a_c=1.0, d_c=2.0, q=1.03. The beam has length L=4.8m, height H=1.2m (c=0.6m) and base B=1m. The Young Modulus is 200GPa and Poisson's ratio 0.3.

Figure 3
Cantilevered beam under axial and transverse load: (a) geometry; (b) model discretized in field nodes.

The stress functions in the plane for the problem are obtained by the differential relationships given in Eqs. (10) to (12):

(40) ${\mathit{\sigma }}_{\mathit{x}}=\frac{{\partial }^2{\varphi }^{\prime }\left(\mathit{x},\mathit{y}\right)}{{\partial }^2\mathit{y}}=\frac{\mathit{N}}{2\mathit{c}}-\frac{3\mathit{Pxy}}{2{\mathit{c}}^3}$

(41) ${\mathit{\sigma }}_{\mathit{y}}=\frac{{\partial }^2{\varphi }^{\prime }\left(\mathit{x},\mathit{y}\right)}{{\partial }^2\mathit{x}}=0$

(42) ${\mathit{\tau }}_{\mathit{xy}}=\frac{{\partial }^2{\varphi }^{\prime }\left(\mathit{x},\mathit{y}\right)}{\partial \mathit{x}\partial \mathit{y}}=\frac{3\mathit{P}}{4\mathit{c}}\left(1-\frac{{\mathit{y}}^2}{{\mathit{c}}^2}\right)$

Fig. 4 indicates the stresses. In x=4.6m (close to loading) both methods present great normal stresses at y=0.6m (centroid), since N and P were applied to the axis of the beam in the subsequent section (x=4.8m). The Airy Stress Function presents linear distribution as it does not consider the Saint-Venant Principle. When x increases (Figs. 4b and 4c), the normal stress gradually shows proportionality with section height according to Eq. (40). In x=4.2m, RPIM and FEM are close to the analytical response, with RPIM showing small divergence at y=0.6m (~3Pa) and subsequent stress reduction (~2.4Pa). For both numerical results with the analytical response.

Figure 4
Numerical and analytical results for example 2: (a) normal stress at x=4.6m; (b) normal stress at x=4.2m; (c) normal stress at x=2.4m; (d) shear stress for all sections.

Fig. 4d shows the shear stresses. According to Equation (42), the results obtained from Airy Stress Function are independent of x in the section. The results obtained numerically consider the Saint-Venant Principle, and the curves gradually approximate the result of Saad (2005) when the section positionx decreases. It should be noted that the RPIM presents better convergence than the FEM in this case for the analytical response (see curves x=4.6m and x=4.2m).

4. Conclusions

This study presented the Radial Point Interpolation Method (RPIM) to evaluate the stress in two-dimensional beams. Formulations based on the Fourier Series Theory and the RPIM were presented. The MQ Radial Basis Functions were used. The numerical results using SAP2000 were also presented. The stress results for the RPIM end FEM were taken at the nodes and not at the Gauss points, which may have caused the difference in the analytical result. RPIM shape parameters are frequently difficult to determine, so they should be adjusted for each problem. Compared to FEM, the solution using RPIM provides satisfactory results for two-dimensional beams. However, a more precise understanding of shape parameters is required. The authors recommend performing a similar study considering the values 0.98 and 1.99 for exponent q, varying the a_c shape parameter and testing different values of N (field nodes). Besides, the authors recommend evaluating the influence of the number of elements (FEM) and the number of field nodes (RPIM) in the results.

Acknowledgments

The authors are grateful to the Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq), the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES), the Fundação de Amparo à Pesquisa do Estado de Minas Gerais (FAPEMIG) and the Liga Acadêmica de Estruturas (LAE) of the Pontifícia Universidade Católica de Minas Gerais.

References

• ASPRONE, D.; AURICCHIO, F.; MONTANINO, A.; REALI, A. Modified finite particle method: multi-dimensional elasto-statics and dynamics. International Journal for Numerical Methods in Engineering, v. 99, p. 1-25, 2014.
• BORESI, A.P.; CHONG, K.; LEE, J. D. Elasticity in engineering mechanics John Wiley &Sons, 2011. 656p.
• FALLAH, A.; JABBARI, E.; BABAEE, R. Development of the Kansa method for solving seepage problems using a new algorithm for the shape parameter optimization. Computers and Mathematics with Applications, v. 77, p. 815-829, 2019.
• FATMI, R. E.; GHAZOUANI, N. Higher order composite beam theory built on Saint-Venant's solution. Part-I: Theoretical developments. International Journal of Composite Structures, v. 93, p. 557-566, 2011.
• GENOESE, A.; GENOESE, A.; BILOTTA, A.; GARCEA, G. A geometrically exact beam model with non-uniform warping coherently derived from the Saint Venant rod. International Journal of Engineering Structures, v. 68, p. 33-46, 2014.
• GENOESE, A.; GENOESE, A.; BILOTTA, A.; GARCEA, G. A mixed beam model with non-uniform warpings derived from the Saint Venant rod. International Journal of Composite Structures, v. 121, p. 87-98, 2013.
• HU, D.; WANG, Y.; LI, Y.; HAN, X.; GU, Y. A meshfree-based local Galerkin method with condensation of degree of freedom for elastic dynamic analysis. Acta Mechanica Sinica, v. 30, p. 92-99, 2014.
• LIU, G.R.; GU, Y.T. An introduction to Meshfree methods and their programming Springer, 2005. 496p.
• LIU, G.R. Meshfree methods: moving beyond the finite element method. CRC Press, 2010. 773p.
• PETROLO, A. S.; CASCIARO, R. 3D beam element based on Saint Venant's rod theory. Computers and Structures, v. 82, p. 2471-2481, 2004.
• MISRA, R. K.; KUMAR, S. Multiquadric radial basis function method for boundary value and free vibration problems. Indian Journal of Industrial and Applied Mathematics, v. 4, p. 1-4, 2013.
• SAAD, M. H. Elasticity: theory, applications and numerics. Elsevier, Inc, 2005. 474p.
• SILVA, R. P. Análise não-linear de estruturas de concreto por meio do método Element Free Galerkin Orientador: Roque Luiz da Silva Pitangueira. 2012. 128f. Tese (Doutorado em Engenharia de Estruturas) - Escola de Engenharia, Universidade Federal de Minas Gerais, Belo Horizonte, 2012.
• TIMOSHENKO, S.; GOODIER, J.N. Teoria da elasticidade 3. ed. Rio de Janeiro: Guanabara Dois, 1980. 545p.
• ZHAO, L.; CHEN, W. Q.; LÜ, C. F. New assessment on the Saint-Venant solutions for functionally graded beams. International Journal of Mechanics Research Communications, v. 43, p. 1-6, 2012.

Publication Dates

• Publication in this collection
  20 Dec 2019
• Date of issue
  Jan-Mar 2020

History

• Received
  28 Dec 2018
• Accepted
  13 Aug 2019