B spline-based method for 2-D large deformation analysis

Engineering Analysis with Boundary Elements 35 (2011) 761–767

Contents lists available at ScienceDirect

Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound

B spline-based method for 2-D large deformation analysis Yanan Liu a,n, Liang Sun a, Feng Xu a, Yinghua Liu b, Zhangzhi Cen b a b

China Special Equipment Inspection and Research Institute, Beijing 100013, China School of Aerospace, Tsinghua University, Beijing 100084, China

a r t i c l e i n f o

abstract

Article history: Received 7 July 2010 Accepted 23 December 2010 Available online 1 February 2011

In this paper, quadratic cardinal B spline functions are used for solution of 2-D large deformation problems. Because the B spline functions are directly used in function approximation, no meshes are needed and the mesh distortion issues in nonlinear analyses are avoided in this method. Using the B spline functions, the solution formulations based on total Lagrangian (TL) approach for two dimensional large deformation problems are established. The numerical examples of 2-D large deformation problems indicate that the B spline method is effective and stable for solving complicated problems. & 2011 Elsevier Ltd. All rights reserved.

Keywords: B spline Large deformation Mesh distortion Basis function

1. Introduction The importance of the computational simulations for large deformation problems has been emphasized and its application is extended to various fields of technology. However, the numerical analysis of large deformation problems is a challenge for researchers in computational mechanics because there has been a very difficult subject due to many characteristics of non-linearities. Although the finite element method (FEM) is a well-established mesh method for modeling nonlinear problems, it often encounters difficulties for nonlinear analyses. One of the main drawbacks of the mesh-based methods like FEM is the mesh distortion issues. Meanwhile, many meshless methods (or particle methods) have been developed to solve large deformation problems [1–5]. Since these particle methods involve only a number of nodal points (particles), they are totally free from the mesh distortion issues which often caused in the FEM analysis, especially for some geometrically nonlinear problems with very large deformation, for which FEM will fail to give reasonable solutions due to the large mesh distortions. In above-mentioned meshless methods, in order to approximate unknown functions, it is key and necessary to construct the so-called shape functions, which are complicated, time consuming and even hard to realize in some special conditions. Furthermore, the complexity of shape functions results in the increase of computational cost and the decrease of computational accuracy. It is desirable to find a new method, which is simple and reasonable to construct shape functions. However, it seems to be a complicated task. So we should resort to some other mathematics tools.

n

Corresponding author. Tel.: +86 10 59068598; fax: + 86 10 59068318. E-mail address: [email protected] (Y. Liu).

0955-7997/$ - see front matter & 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.enganabound.2010.12.006

Wavelet is a powerful mathematics tool in solving many problems in science and engineering. In recent years, there has been an increasing interest in wavelet-based methods due to their successes in some applications. Wavelet-based numerical method has been developed. At present, there are mainly three kinds of wavelet-based numerical methods: the wavelet finite element method [6,7], the wavelet collocation method [8,9] and the wavelet-Galerkin method [10,11]. Most of these wavelet based methods are based on orthogonal or biorthogonal wavelets. These wavelets are not easy to work with for numerical analysis because they do not have explicit functional forms or sufficient regularity. Some special formulae about computation of these wavelets [12,13] can only be used for problems on simple domains. The DB wavelet-based method [14–16] can simulate complicated deformation. However, the computational efficiency of this method is relatively low because of the complexity of DB wavelet. As the scaling function of spline wavelet, cardinal B spline functions have been used in numerical simulation for decades years. In fact, a considerable body of literatures now exists on the application of uniform and non-uniform B spline techniques to the solution of partial differential equations (PDEs) and mechanics problems. The recent studies of B spline method can be found in some articles [17–23]. The B spline basis functions have compact support and lead to banded stiffness matrices. They can be used to construct piecewise approximations that provide higher order of continuity depending on the order of the polynomial basis. For example, quadratic B-spline basis can provide solutions with continuous gradients (C1 continuity) throughout the domain. The B-spline basis functions form a partition of unity, which is an important property for convergence of the approximate solutions. As they are polynomials, accurate integration can be performed by using the Gauss quadrature. The B-spline approximation has good reproducing properties; thus, it is able to represent constant strains

762

Y. Liu et al. / Engineering Analysis with Boundary Elements 35 (2011) 761–767

exactly. Compared to other orthogonal or biorthogonal wavelets scaling functions and the shape functions constructed by meshless methods, B spline functions are more simple and easy to work with for numerical analysis. In addition, because B spline functions can be directly used as basis functions, no meshes are needed in function approximation. Therefore, B spline-based method is also free from the mesh distortion issues and is potential in simulation of large deformation problems. However, the studies for the nonlinear analyses by the B spline method are very few. The main aim of this paper is to introduce the cardinal B spline functions to solve 2-D large deformation problems. The cardinal B spline functions which are uniform B spline basis can be directly used to approximate the unknown field functions instead of constructing shape functions as done in FEM and conventional meshless methods. Using the total Lagrangian (TL) approach which refers all variables to the initial configuration, we build formulations for solution of 2-D large deformation problems with B spline functions. In addition, Lagrange multipliers can be introduced to impose essential boundary condition in problems with complicated domain. The numerical examples of 2-D geometrically nonlinear analysis are given to illustrate the stability and the effectiveness of the present method in large deformation problems.

2. Approximation of 2-D functions by B spline The m degree cardinal B spline is defined as Z 1 Nm ðxÞ ¼ Nm1 N1 ¼ Nm ðxtÞdt, m Z 2

ð1Þ

0

where ( N1 ðxÞ ¼

1,

x A ½0,1Þ

0,

else

ð10Þ

k,l j,l i,k where, Nm ðxÞNm ðyÞ are 2-D B spline basis functions, i and j are, respectively, the scales of x direction and y direction in approximation. For 2-D problems in general domains O, the all 2-D B spline basis functions which meet the following condition should be used in function approximation: i,k j,l SuppNm ðxÞNm ðyÞ \ O a0

ð11Þ

Fig. 1 shows the domain overlaid by 2-D quadratic B spline functions used for approximation of function defined in O. The domain consists of small rectangles and with the length and height of every small rectangle are, respectively, i and j. In practical computations, we only need to consider the part of 2-D B spline functions lie in small rectangles inside or partially inside the domain of analysis. As shown in Fig. 2(a), the region for practical computations consists of inner small rectangles and boundary small rectangles. In order to obtain high accuracy, the integration in total computation process should be evaluated in every small rectangle. Especially, the regions of analysis that are inside boundary small rectangles should be triangulated or quadrangulated and the integration is computed by performing quadrature on the triangles or quadrangles. This approach leads to accurate integration computation in boundary small rectangles. The subdividing regions for integration are shown in Fig. 2(b). Similar to finite element method and meshless methods, the approximation equation can be written as uðx,yÞ ¼

n X

h Nm ðx,yÞch

ð12Þ

h¼1

ð2Þ

The major properties of cardinal B spline are

the domain 2-D B spline functions overlay

SuppNm ¼ ½0,m

ð3Þ

u ðxÞ ¼ Nm1 ðxÞNm1 ðx1Þ Nm

ð4Þ

Nm ðxÞ ¼

by B spline function can be expressed as X i,k j,l uðx,yÞ ¼ ck,l Nm ðxÞNm ðyÞ

x mx Nm1 ðxÞ þ Nm1 ðx1Þ m1 m1

Ω

ð5Þ

From the above definition, the quadratic cardinal B spline is expressed as 8 2 x =2 x A ½0,1Þ > > > < 3=2 þ 3xx2 x A ½1,2Þ ð6Þ N3 ðxÞ ¼ 2 > > > 9=23x þ x =2 x A ½2,3 : 0 else B spline functions can be used as basis functions to approximate the function u defined on interval [a,b]. X i,k uðxÞ ¼ ck Nm ðxÞ ð7Þ

Fig. 1. The domain overlaid by 2-D quadratic B spline functions used for approximation in O.

boundary small rectangle

k i,k ðxÞ ¼ Nm ð1=ixkÞ and i denotes the scale in approximawhere, Nm i,k tion. According to properties of B spline, the support of Nm ðxÞ is i,k ¼ ½ik,iðm þ kÞ SuppNm

In approximation Eq. (7), the B spline functions satisfy the following condition: i,k \ ½a,b a 0 SuppNm

Ω

ð8Þ i,k ðxÞ Nm

should

ð9Þ

The basis functions for the higher-dimensional problems are constructed by taking the product of the basis functions for 1-D B-spline. In this case, the approximation of 2-D function u(x,y)

inner small rectangle

boundary regions should be triangulated or quadrangulated

Fig. 2. The domain for practical computation: (a) inner small rectangles and boundary small rectangles and (b) domain subdivided for integration.

Y. Liu et al. / Engineering Analysis with Boundary Elements 35 (2011) 761–767 j,l h i,k where, Nm ðx,yÞ ¼ Nm ðxÞNm ðyÞ is similar to shape functions in finite element method and meshless methods, ch are the generalized h displacement related to Nm ðx,yÞ and n is the number of 2-D B spline functions used in approximation.

3.1. Basic formulations Assuming that the analysis starts at time 0 and all state variables that satisfy equilibrium are already known up to time t, further loading and deformation will satisfy the equilibrium at time t + Dt. Using virtual displacement principle which is equal to equilibrium equation, we have Z t þ Dt sij dt þ D teij t þ Dt dV ¼ t þ Dt W ð13Þ V

where, t þ Dt sij is the Cauchy stress tensor, t þ Dt W are the external virtual work from body forces and tractions, dt þ D teij is the variation of strain tensor consistent with the virtual displacement dui, t þ Dt V is the volume of the body, and the quantities are all measured at time t+ Dt. In this paper, the total Lagrange format is used for solution of large deformation. The reference configuration in the variational equation is selected as the initial configuration. In this case, Eq. (13) can be transformed to the initial configuration Z t þ Dt t þ Dt Eij 0 dV ¼ t þ Dt W ð14Þ 0 Sij d0 0V

t þ Dt 0 Sij

t þ Dt 0 Eij

and are, respectively, the second Piola– where, Kirchhoff stress and the Green strain tensor at time t + Dt, the left subscript 0 denotes that the initial configuration is used as reference configuration. According to the assumption, we have the following incremental relationships: t þ Dt 0 Sij

¼ t0 Sij þ 0 Sij

ð15Þ

t þ Dt 0 Eij

¼ t0 Eij þ 0 Eij

ð16Þ

where, t0 Sij is the second Piola–Kirchhoff stress tensor at time t, 0 Sij and 0 Eij are, respectively, the increments of second Piola– Kirchhoff stress tensor and Green strain tensor from time t to t + Dt. Because t0 Sij and t0 Eij are known, we have

dt þ D0t Eij Z 0V

¼ d0 Eij

0V

t 0 BL1

0V

t 0 Sij d0 eij 0 dV

ð20Þ

t 0 ui

n X

¼

0V

t 0 cij Nj

0 ui

¼

j¼1

n X

0 cij Nj

ði ¼ 1,2,3Þ

ð21Þ

j¼1

where, t0 ui is the displacement at time t, 0 ui is the displacement increment from time t to t + Dt, Nj consists of B spline functions, t 0 cij is the coefficient at time t, 0 cij is the increment coefficient and n is the number of B spline functions used in approximation. Introducing the approximation (21) into Eq. (20), we can obtain a group of algebra equations ¼ t þ Dt Q t0 F

t 0K 0c

ð22Þ

The computation of above formulae is as follows: t 0K

¼ t0 K L þ t0 K NL Z

t 0K L

¼ 0V

t 0F

t T t 0 BL 0 D 0 BL 0 dV

Z

t 0 K NL

¼ 0V

Z ¼ 0V

ð23Þ ð24Þ

t T t t 0 BNL 0 S 0 BNL 0 dV

ð25Þ

t Tt ^ 0 BL 0 S 0 dV

ð26Þ

where, t0 BL can be expressed as t 0 BL

¼ t0 BL0 þ t0 BL1

ð27Þ

In 2-D plain problems 9 8 n > >X t > > (t ) > N ðx,yÞ c i > > 0 xi > = < u ðx,yÞ x 0 i¼1 t ¼ t 0u ¼ n X > > 0 uy ðx,yÞ > > > > > Ni ðx,yÞt0 cyi > ; :

ð28Þ

i¼1

( 0u

¼

0 ux ðx,yÞ

) ¼

0 uy ðx,yÞ

8 n X > > > Ni ðx,yÞ0 cxi > < i¼1

9 > > > > =

ð29Þ

n X > > > > > > > Ni ðx,yÞ0 cyi > ; :

t 0 Sij d0 eij 0 dV

2 @N

ð18Þ t 0 BL0

6 6 ¼6 4

@x

0



@Nn ðx,yÞ @x

0

@N1 ðx,yÞ @y



0

@f1 ðx,yÞ @y

@N1 ðx,yÞ @x



@Nn ðx,yÞ @y

1 ðx,yÞ

L21 @N1@xðx,yÞ



L11 @Nn@xðx,yÞ

L21 @Nn@xðx,yÞ

L12 @N1@yðx,yÞ

L22 @N1@yðx,yÞ



L12 @Nn@yðx,yÞ

L22 @Nn@yðx,yÞ

L11 @N1@yðx,yÞ þL12 @N1@xðx,yÞ

L21 @N1@yðx,yÞ þL22 @N1@xðx,yÞ



L11 @Nn@yðx,yÞ þ L12 @Nn@xðx,yÞ

L21 @Nn@yðx,yÞ þ L22 @Nn@xðx,yÞ

stress–strain relation in terms of the second Piola–Kirchhoff stress increment and the Green strain increment may be defined as

7

@Nn ðx,yÞ 7 7 @y @Nn ðx,yÞ @x

5

ð30Þ

7 7 7 5

ð31Þ

where

L11 ¼

N X @Ni ðx,yÞ i¼1

ð19Þ

where, 0 Dijkl is the tangent constitutive tensor. With further linearization of Eq. (18), the approximate expression of virtual

3

0

3

L11 @N1@xðx,yÞ

0 Sij ¼ 0 Dijkl 0 Elk

0V

i¼1

where, 0 eij and 0 Zij are, respectively, the linear part and the nonlinear part of Green strain increment from time t to t + Dt. The

6 6 ¼6 4

0V

ð17Þ

In this case, Eq. (14) can be rewritten as Z Z t t þ Dt W 0 Sij d0 Zij 0 dV ¼ 0 Sij d0 Eij 0 dV þ

2

displacement principle is Z Z Z t t þ Dt W 0 Dijkl 0 eij d0 eij 0 dV þ 0 Sij d0 Zij 0 dV ¼

Using B spline method to solve large deformation problems, the approximation of displacement functions and increment of displacement functions can be expressed as

3. Numerical implementation

t þ Dt

763

L21 ¼

@x

N X @Ni ðx,yÞ i¼1

@x

t b 0 cxi

L12 ¼

N X @Ni ðx,yÞ i¼1

t b 0 cyi

L22 ¼

@y

N X @Ni ðx,yÞ i¼1

@y

t b 0 cxi

t b 0 cyi

ð32Þ

764

Y. Liu et al. / Engineering Analysis with Boundary Elements 35 (2011) 761–767

2 t 0 BNL

@N1 ðx,yÞ 6 @N @xðx,yÞ 6 1 6 @y

¼6 6 6 4

0 0

2t

0 S11 6t 6 0S 6 21 t 0S ¼ 6 6 0 4 0



0



@Nn ðx,yÞ @x @Nn ðx,yÞ @y

@N1 ðx,yÞ @x @N1 ðx,yÞ @y



0



0

t 0 S12 t 0 S22

0

0

t 0 S11 t 0 S21

0

0

2t

t ^ 0S

0

0

0

3

7 0 7 7 7 @Nn ðx,yÞ 7 7 @x 5 @N ðx,yÞ

ð33Þ

n

and I1, I2 and I3 are the first, second and third invariants of Green 1=2 deformation tensor G, respectively, J ¼ I3 and K is the bulk modulus to be measured experimentally. The second Piola– Kirchhoff stress calculated from the strain energy density function is     @W 1 2 1=3 2=3 þ 2q2 I3 þ PJG1 ¼ 2q1 I3 dij  I1 G1 I1 dij Gij  I2 G1 ij ij ij @Eij 3 3

@y

Sij ¼

3

7 0 7 7 7 t 0 S12 7 5 t 0 S22

ð46Þ ð34Þ

where P¼

3

0 S11 6 7 ¼ 4 t0 S22 5

ð35Þ

~ @W ¼ KðJ1Þ @J

qn ¼

t 0 S12

@W , @In

ð47Þ

n ¼ 1,2

ð48Þ

3.4. Followed loads 3.2. Imposition of essential boundary condition For some problems with complicated boundary, Lagrange multiplier method is efficient for the realization of essential boundary condition. Introducing the Lagrange multiplier into Eq. (20), we have Z Z Z 0 0 t Dijkl 0 eij d0 eij dV þ Sij d0 Zij dV þ dðli ð0 ui 0 ui ÞdS 0 0 0 0 0 V V Su Z 0 t ¼ t þ Dt W  ð36Þ 0 Sij d0 eij dV 0

V

where, li is Lagrange multiplier and expressed as X j L li ¼ Nm cij

ð37Þ

j¼1

The stiffness matrix and load vector in Eq. (22) should be changed "t # KL 0K t ð38Þ 0K B ¼ K TL 0 " 0 cB

¼

0c

# ð39Þ

0 cL

2

t þ Dt 0Q t þ Dt 4 0QB ¼ t þ Dt 0QL

" t þ Dt 0 FB

¼

t þ Dt 0F

5

ð40Þ

ð41Þ

3.3. The constitutive model In this paper, the so-called hyperelastic material is considered in some examples. A strain energy density function that decouples the distortional and dilatational deformation [24] is ~ ðJÞ WðI1 ,I2 ,JÞ ¼ WðI1 ,I2 Þ þ W

ð42Þ

~ ðJÞ are the distortional and dilatational where, WðI1 ,I2 Þ and W strain energy density functions given by Apq ðI1 3Þp ðI2 3Þ

q

ð43Þ

pþq ¼ 1

~ ðJÞ ¼ 1 KðJ1Þ2 W 2

ð44Þ

1=3

I1 ¼ I1 I3

,

2=3

I2 ¼ I2 I3

ð45Þ

dS

ð49Þ

t þ Dt 0 t þ Dt 1 0 0 pFik vi dS 0 t k dS ¼ J

ð50Þ

where, 0 vk is the component of unit outward normal on initial configuration, Fik is the deformation gradient and J ¼det(Fik). 3.5. Method of solution In this work, N-R iteration [25] is used to solve Eq. (22). In this approach, Eq. (22) can be rewritten as t þ Dt ðlÞ ¼ t þ Dt Q B 0t þ Dt F ðlÞ FL B 0

ðl ¼ 0,1,2,. . .Þ

ð51Þ

and ðlÞ ¼ 0t þ Dt cðlÞ B þ D0 c B

ð52Þ

The superscript l denotes iteration counters and t þ D0t F ðlÞ L is the imbalance force vector related to Lagrange multiplier and can be expressed as 2 3 K L t þ D0t cðlÞ L t þ Dt ðlÞ 4 5 ð53Þ 0 FL ¼ K TL t þ D0t c ðlÞ for l¼0 t þ Dt ð0Þ 0 KB

¼ t0 K B

ð54Þ

D0 cð0Þ B ¼ 0 cB

ð55Þ

t þ Dt ð0Þ 0 cB

¼ t0 cB

ð56Þ

t þ Dt ð0Þ 0 FB

¼ t0 F B

ð57Þ

t þ Dt ð0Þ 0 FL

¼

"

where

t þ Dt

where, t þ Dt p is the pressure on surface of configuration and t þ Dt vk is the component of unit outward normal on configuration at time t+ Dt. Using the transformation relationship, the followed loads can be expressed as:

t þ Dt ðl þ 1Þ 0 cB

0

WðI1 ,I2 Þ ¼

dTk ¼ t þ Dt p t þ Dt vk

ðlÞ t þ Dt ðlÞ 0 K B D0 c B

3

#

1 X

In practical problems, the pressure loads which are dependent on deformed states should be considered. In this case, the load on tiny area of configuration at time t + Dt can be expressed as

K L t0 c L K TL t0 c

# ð58Þ

Y. Liu et al. / Engineering Analysis with Boundary Elements 35 (2011) 761–767

The load vector 0t þ Dt F ðlÞ related to follows. From below formulae Z ðlÞ t þ Dt ðlÞ 0 Sij d0 eij 0 dV 0

t þ Dt ðlÞ FB 0

is computed as

765

p/2

h

ð59Þ

V

p/2

and

d0 eij ¼

l

1 ðd u þ d0 uj,i þ 0t þ Dt u k,i d uk,j þ t þ D0t u k,j d 0 uk,i Þ 2 0 i,j

ð60Þ

Fig. 4. A cantilever beam under uniform loads on its upper and bottom boundary.

where, the superscript l shows that stress and strain are calculated by t þ D0t c ðlÞ , we have Z t þ Dt ðlÞ T t þ Dt ^ ðlÞ 0 t þ Dt ð61Þ S dV 0 BL 0 0 F ðlÞ ¼ 0

0.2 0.0

V

ðlÞ

t þ Dt ^ In 2-D problems, t þ D0t BðlÞ are, respectively, corre0S L and t ^ t sponding to 0 BL in (27) and 0 S in (35). However, they are computed from t þ D0t c ðlÞ .

initial configuration

Y

-0.2

deformation under vertical load p = 10

-0.4

4. Numerical examples

-0.6

In this part, numerical simulation of some large deformation problems are presented using B spline method. Quadratic cardinal B spline functions are used for approximation in all numerical examples. The results are compared with those calculated by finite element method or analytical results to show the validity of the proposed method. For simplification, the units are omitted in this paper.

-0.8

deformation under followed load p = 10

-1.0 0.0

0.2

0.4

0.6

0.8

1.0

X Fig. 5. Deformation of a beam under vertical load and followed load.

4.1. Tension along axial direction A plain elastic panel under tension along axial direction is analyzed by B spline method. The panel is made of linear elastic material. When v¼0, the analytical result is    u u 1 u2 p ¼ EA 1 þ þ ð62Þ L L 2 L Fig. 3 shows the comparison of analytical and numerical results computed by B spline method, it can be found that the two results agree very well.

Table 1 Comparison of deflection at middle point of free end. Load value

2

4

6

8

10

Finite element method Vertical 0.238 Followed 0.248

0.425 0.472

0.555 0.651

0.643 0.776

0.704 0.841

B spline method Vertical Followed

0.430 0.474

0.560 0.657

0.650 0.782

0.712 0.844

0.241 0.248

4.2. Large displacement of cantilever beam Fig. 4 depicts a beam under uniform loads on its upper and bottom boundary [25]. The length of beam is l ¼1, and the height is h¼0.1.

400

analytical B spline method

350

p

300

p/EA

250 200 150 100 50 0 -50 0

2

4 u/L

6

Fig. 3. The comparison of numerical and analytical results.

8

Two kinds of loads are considered in this example, one is vertical load and the other is followed load. The beam is made of linear elastic material with material constants given by E¼ 1.2  104 and v¼ 0.2. The scale used in approximation is (1/20)  (1/20) and the total number of 2-D B spline basis functions is 88. Fig. 5 depicts the deformation of beam under vertical load and followed load with the load value p¼10. Table 1 shows the comparison of the deflection at middle point of free end computed by finite element method and B spline method under followed load and vertical load. We can see the results calculated by the two methods agree well. Fig. 6 depicts a beam which is subjected to an end moment by applying linear normal traction with zero resultant normal force, on the free end of the beam. Also note that the linear normal traction moves with the beam so that it is tangent to the neutral surface. The beam is made of linearly elastic material with material constants given by E ¼100 and v ¼0. The beams of size, 1.0  0.1667, 1.0  0.1 and 1.0  0.05 are analyzed in this example. The three beams are, respectively, under bending moment M¼0.0606, 0.01309 and 0.001636. From the analytical result of 1-D large deformation beam, the three beams will turn 1/4 circle under corresponding bending moment value. Fig. 7 describes the pure bending deformation of the beam of size 1.0  0.1667. It is obvious that this figure catches the pure bending deformation state. Furthermore, the deformation of this 2-D beam slightly overrunning 1/4 circle agrees with the practical

766

Y. Liu et al. / Engineering Analysis with Boundary Elements 35 (2011) 761–767

p

Fig. 6. A beam under bending moment on right boundary.

1.0 M = 0.0606

0.8

Y

0.6 0.4 0.2

made of Neo-Hookean hyperelastic material and sandwiched between metal layers. The material constant is given by K¼1000, A10 ¼0.273. The right boundary moves leftward while the left boundary is fixed. The scale used in approximation is (1/34)  (1/ 34) and the total number of 2-D B spline basis functions is 684. The deformations of selected compression scale are given in Fig. 9. The second example is the tension of a plane strain hyperelastic bonded rubber of dimension 1.0  2.0 of which the problem statement is given in Fig. 10. The rubber is also made of Neo-Hookean hyperelastic material with material properties given by K ¼10,000, A10 ¼ 80.273. The right boundary moves rightward and the left surface is fixed. The scale used in approximation is (1/17)  (1/17) and the total number of 2-D B spline basis functions is also 684. The deformations of selected tension scale are given in Fig. 11. The two examples of plain strain hyperelastic rubber are also analyzed by the FEM. It is found that the FEM converges very slowly once the compression is more than 50% or the tension is more than 700%. It demonstrates that B spline method is more effective than the FEM for the large deformation problems. It is

0.0 0.0

0.2

0.4

0.6

0.8

1.0

X Fig. 7. Pure bending deformation of beam under bending moment.

Table 2 Comparison of displacement at middle point of free end of beams with different dimension. Size

B spline method Analytical result

1.0  0.1667

1.0  0.1

ux

uy

ux

uy

ux

uy

0.5983

0.649

0.6225 0.637

0.642

0.628

0.6416

fixed boundary

50% compression

70% compression

Fig. 9. Deformation for compression of a plain strain bonded hyperelastic rubber.

1.0  0.05

fixed boundary

free surface

free surface

2.0 0.5

1.0

free surface Fig. 8. Models of a plain stain hyperelastic rubber under compression.

free surface

1.0

Fig. 10. Models of a plain strain hyperelastic rubber under tension.

case. Table 2 demonstrates the comparison of displacement at middle point of free end of beams with different dimension. It can be found that the deformation of the beam with larger length-toheight ratio agrees better with the 1-D analytical results. These results show that the present B spline method can accurately capture the pure bending deformation of beam.

800% tension

4.3. Large deformation of bonded rubber A plain-strain bonded rubber of dimension 1.0  0.5 is analyzed at first. The problem statement is illustrated in Fig. 8. The rubber is

1059% tension Fig. 11. Deformation for tension of a plain strain bonded hyperelastic rubber.

Y. Liu et al. / Engineering Analysis with Boundary Elements 35 (2011) 761–767

obvious that the present B spline method is free from mesh distortion in large deformation analysis.

5. Conclusion In this paper, quadratic cardinal B spline functions are used for simulation of 2-D large deformation problems. Because no meshes are needed in function approximation using the B spline functions, the mesh distortion issues in nonlinear analysis do not appear in this method. Based on total Lagrangian (TL) approach, the solution formulations for two dimensional large deformation problems are established and the imposition of essential boundary condition and solution method are discussed. Numerical examples for 2-D elastic and hyperelastic large deformation problems illustrate that B spline method is effective and stable for solving complicated problems. References [1] Gu YT, Wang QX, Lam KY. A meshless local Kriging method for large deformation analyses. Computer Methods in Applied Mechanics and Engineering 2007;196:1673–84. [2] Chen J, Pan C, Wu C, Liu WK. Reproducing Kernel particle methods for large deformation analysis of non-linear structures. Computer Methods in Applied Mechanics and Engineering 1996;139:195–227. [3] Jun S, Liu WK, Belytschko T. Explicit reproducing kernel particle methods for large deformation problems. International Journal for Numerical Methods in Engineering 1998;41:137–66. [4] Liu WK, Jun S. Multi-scale reproducing kernel particle methods for large deformation problems. International Journal for Numerical Methods in Engineering 1998;4:1339–62. [5] Li S, Hao W, Liu WK. Mesh-free simulations of shear banding in large deformation. International Journal of Solids and Structures 2000;377:185–7206. [6] Ko J, Kurdila AJ, Pilant MS. A class of finite element methods based on orthonormal compactly supported wavelets. Computational Mechanics 1995;16:235–44. [7] Ma J, Xue J, Yang S, He Z. A study of the construction and application of a Daubechies wavelet-based beam element. Finite Elements in Analysis and Design 2003;39:965–75. [8] Bertoluzza S, Naldi G. A wavelet collocation method for the numerical solution of partial differential equations. Applied and Computational Harmonic Analysis 1996;3:1–9.

767

[9] Vasilyev OV, Paolucci S, Sen M. A multilevel wavelet collocation method for solving partial differential equations in a finite domain. Journal of Computational Physics 1995;120:33–47. [10] Amaratunga K, Williams JR. Wavelet-Galerkin solutions for one dimensional partial differential equation. International Journal for Numerical Methods in Engineering 1994;37(16):2703–16. [11] Diaz AR. A wavelet-Galerkin scheme for analysis of large-scale problems on simple domains. International Journal for Numerical Methods in Engineering 1999;44:1599–616. [12] Latto A, Resnikoff HK, Tenenbaum E. The evaluation of connection coefficients of compactly supported wavelets. In: Proceedings of the USA–French workshop on wavelets and turbulance, Princeton University; 1991. [13] Restrepo JM, Leaf GK. Inner product computations using periodized Daubechies wavelets. International Journal for Numerical Methods in Engineering 1997;40:3557–78. [14] Liu Y, Qin F, Liu Y, Cen ZA. Daubechies wavelet-based method for elastic problems. Engineering Analysis with Boundary Element 2010;34:114–21. [15] Liu Y, Qin F, Liu Y, Cen Z. The 2D large deformation analysis using Daubechies wavelet. Computational Mechanics 2010;45:179–87. [16] Liu Y, Liu Y, Cen Z. Multi-scale Daubechies wavelet-based method for 2-D elastic problems. Finite Element Analysis and Design 2010, doi:10.1016/ j.finel.2010.11.004. [17] Naginoa H, Mikamia T, Mizusawa T. Three-dimensional free vibration analysis of isotropic rectangular plates using the B-spline Ritz method. Journal of Sound and Vibration 2008;317:329–53. [18] Caglar N, Caglar H. B-spline method for solving linear system of second-order boundary value problems. Computers and Mathematics with Applications 2009;57:757–62. [19] Kagan P, Fischer A, Bar-Yoseph PZ. Mechanically based models: adaptive refinement for B-spline finite element. International Journal for Numerical Methods in Engineering 2003;57:1145–75. [20] Burla RK, Kumar AV. Implicit boundary method for analysis using uniform B-spline basis and structured grid. International Journal for Numerical Methods in Engineering 2008;76:1993–2028. [21] Lakestani M, Dehghan M. Numerical solution of Fokker–Planck equation using the cubic B-spline scaling functions. Numerical Methods for Partial Differential Equations 2009;25:418–29. [22] Lakestani M, Dehghan M. Numerical solution of Riccati equation using the cubic B-spline scaling functions and Chebyshev cardinal functions. Computer Physics Communications 2010;181(5):957–66. [23] Dehghan M, Lakestani M. The use of cubic B-spline scaling functions for solving the one-dimensional hyperbolic equation with a nonlocal conservation condition. Numerical Methods for Partial Differential Equations 2007;23(6):1277–89. [24] Chang TYP, Saleeb AF, Li G. Large strain analysis of rubber-like materials based on a perturbed Lagrangian variational principle. Computational Mechanics 1991;8:221–33. [25] Wang XC, Shao M. Finite element method. Beijing: Tsinghua University Press; 2003.