





































Jtam-A4.dvi


JOURNAL OF THEORETICAL

AND APPLIED MECHANICS

53, 1, pp. 179-194, Warsaw 2015
DOI: 10.15632/jtam-pl.53.1.179

THE EFFECTIVE STRATEGY FOR CALCULATING STRAINS AND

STRESSES IN ELASTO-PLASTIC TORSION OF A BAR

Eugeniusz Zieniuk, Agnieszka Bołtuć

University of Bialystok, Faculty of Mathematics and Computer Science, Białystok, Poland

e-mail: ezieniuk@ii.uwb.edu.pl; aboltuc@ii.uwb.edu.pl

The paper presents the general approximation strategy for derivatives of solutions obtained by
the parametric integral equation system (PIES). The proposed approach is applied to calculating
derivatives in linear and nonlinear torsion of a bar. The nonlinear problem is solved iteratively
and requires repeated calculating of integrals in order to determine stresses or strains. Due to
elimination of the repeated integration, the proposed strategy increases the efficiency of calcula-
tions. Furthermore, it ensures high accuracy of solutions to the considered problem, which has
been tested in comparison to the analytical and/or numerical solutions reported in the literature
obtained by other methods.

Keywords: PIES, torsion of a bar, approximation, stresses, strains

1. Introduction

The parametric integral equation system (PIES) is a method developed as an alternative to
other well-known numerical methods as the finite element method (FEM) (Ameen, 2005; Zien-
kiewicz, 1977), boundary element method (BEM) (Aliabadi, 2002; Ameen, 2005) or so-called
meshless methods (Liu, 2002). This method has confirmed its advantages in solving bounda-
ry value problems modelled by various linear partial differential equations: Laplace’s (Zieniuk,
2001), Helmholtz’s (Zieniuk and Bołtuć, 2006a), Navier-Lamé’s (Zieniuk and Bołtuć, 2006b;
Zieniuk, 2013; Bołtuć and Zieniuk, 2011, 2013). Successfully verified benefits are as follows: no
discretization of the boundary and the domain which is replaced by effective modelling using
curves and surfaces, high accuracy of obtained solutions, lower computational complexity. These
advantages are very encouraging to generalize the method to other problems, including more
complex (e.g. nonlinear).

The initial attempts to generalize andapply thePIESmethod for solvingnonlinearboundary
valueproblems inmechanics showed theauthors the typeand scale of problems that occurduring
their application. The first is the necessity of repeated calculating of integrals, often singular, in
order to determine the results (often derivatives of the solutions) atmany points of the domain.
Accurate computation of singular integrals (even with a strong singularity, e.g. in the elasto-
plastic mechanical problems) is a serious problem, which determines the effectiveness of the
whole PIES method (Zieniuk, 2013). It should bementioned that these problems are not more
significant than in BEM. Moreover, the numerically computed solutions and especially their
derivatives obtained at points close to the boundary are less accurate because the singularities
occur in thevicinity of theboundary.Errors generated in these cases, accumulated in an iterative
process used for solving the nonlinear problems, affect the accuracy of the final results.
Therefore, in the research carried out by the authors, the new strategy has been searched,

which can eliminate the above mentioned drawbacks. The authors have developed the strategy
based on the interpolation of solutions obtained directly from PIES and the integral identity,
and differentiation of the so-obtained interpolation polynomial. On the basis of the obtained



180 E. Zieniuk, A. Bołtuć

interpolation polynomial, we can determine derivatives of any order with respect to any variable
in a simple manner using a single set of values (solutions on the boundary and in the domain).
The thus obtained derivatives of the polynomial allow calculating derivatives in a continuous
way at any point of the domain, and also the boundary. Verification of the mentioned strategy
on the linear elasticity problems is presented by Bołtuć and Zieniuk (2013). It should be em-
phasized that the proposed strategy can be applied also in BEM, but it does not guarantee the
highest efficiency. About the effectiveness of the strategy firstly decides the way of modelling
the considered domain and the resulting number of algebraic equations to be solved in order to
achieve the interpolation. On the other hand, on the accuracy of the approximation the greatest
impact has the optimal arrangement of interpolation nodes at which the solutions are obtained.
The effectiveness of performing the interpolation determines freedom in automatic generation
of an optimal nodes arrangement. Such an arrangement in combination with PIES is shown in
the following paper.
Due to the complexity of nonlinear problems, we decided to verify PIES with the appro-

ximation strategy on the most elementary problem in the field of nonlinear solid mechanics.
Therefore, the verification of the presented strategy has been done on nonlinear torsion of the
bar. The elementariness of this issue as compared to other nonlinear problems ensures proper
verification of themathematical apparatus applied to the proposed strategy. Numerical solution
of the elasto-plastic bar torsion has been frequently studied (including inverse problems) using
methods such as the finite differencemethod (FDM) (Mendelson, 1968), FEM (Mamedov, 1995;
Yamada et al., 1972), BEM (Mendelson, 1975; Sapountzakis and Tsipiras, 2009) and recently
using the method of fundamental solutions (MFS) (Kołodziej and Gorzelańczyk, 2012), which
significantly facilitates the verification of the results obtained by the applied strategy.
Thepurposeof thispaper is to apply the strategy for approximationofderivatives of solutions

and toverify its effectiveness on theexample of the elasto-plastic torsionproblemsolvedbyPIES.
This means that on the basis of the numerical solutions to the elasto-plastic torsion problem
obtained by the PIES method and using the proposed approximation strategy, we calculate
strains and stresses. The analytical solution in the considered example exists only for the elastic
torsion. In the case of plasticity,wearedealingonlywithnumerical solutions, but this elementary
problem has been solved by various methods, which results in better verification. First, the
problem has been solved as linear, allowing one to verify the PIES method with analytical
solutions. Then, the effectiveness of the strategy has been tested in an iterative process used to
solve the nonlinear problems, where obtained results have been compared with other numerical
solutions.

2. Formulation of problem

Consider a prismatic bar subjected to a twisting moment as shown in Fig. 1.
The elasto-plastic problem of torsion can be formulated in several ways. In this paper, we

use a formulation in terms of thewarping function, because it is physicallymostmeaningful and
most universal, and we consider bars of rectangular cross-sections. The main equations are as
follows (Mendelson, 1975):
a) The equilibrium equation
Assume in accordance with the bar torsion theory that the stress state is defined inde-

pendently of x3 by components τ13 and τ23 (Fig. 1b), and the remaining elements are zero
σ11 = σ22 = σ33 = τ12 = 0. The equilibrium differential equations in this case are reduced to
one (Sokolowski, 1957)

∂τ13
∂x1
+
∂τ23
∂x2
=0 (2.1)

where τ13 and τ23 are the shear stresses acting on the cross-section.



The effective strategy for calculating strains and stresses... 181

Fig. 1. (a) Prismatic bar under the action of a twisting moment, and (b) the cross-section with the
shear stresses

b) Stress-strain relations

τ13 =2G(ε13 −εp13) τ23 =2G(ε23−ε
p
23) (2.2)

whereG is the shear modulus, ε13 and ε23 are the total strains and ε
p
13 and ε

p
23 are the plastic

shear strains.
c) Saint-Venant’s relations

ε13 =
1

2

(

−αx2+
∂w

∂x1

)

ε23 =
1

2

(

αx1+
∂w

∂x2

)

(2.3)

where α is the angle of twist per unit length and w is the warping function. The warping
function is assumed to be the same for all cross-sections, therefore it is a function of x1 and x2
only (Ameen, 2005).
Substituting equations (2.3) into (2.2) and the result into (2.1), we obtain

∂2w

∂x21
+
∂2w

∂x22
=F(x1,x2) in Ω (2.4)

whereΩ is the region of the cross-section and

F(x1,x2)= 2
(∂ε

p
13

∂x1
+
∂ε
p
23

∂x2

)

(2.5)

Theplastic strains appearing in (2.2) and (2.5) can bepresented by the following expressions

ε
p
13 =

εp
εet
ε13 ε

p
23 =

εp
εet
ε23 (2.6)

where

εet =
2√
3

√

ε213+ε
2
23 εp =

εet− 23(1+ν)ε0
1+ 2

3
(1+ν) m

1−m

(2.7)

in the case of linear strain hardening,where ε0 is the yield strain,m is the linear strain hardening
parameter and ν Poisson’s ratio. Formula (2.7)2 holds for 0¬m< 1, wherem is equal to 0 for
a perfectly plastic case. In an elastic case, m is equal to 1 and εp =0.
Calculation of the function defined by formula (2.5) requires the derivatives of the plastic

strains. Therefore, we should differentiate expression (2.6)1 with respect to the x1-variable,
obtaining the equation

∂ε
p
13

∂x1
=
Bε13

∂εet
∂x1

Aε2et
+
(εet−B
Aεet

)∂ε13
∂x1

(2.8)



182 E. Zieniuk, A. Bołtuć

where εet and ε13 are expressed by formulas (2.7)1 and (2.3)1, and the remaining functions are
defined using the following expressions

∂εet
∂x1
=
4

3εet

(

ε13
∂ε13
∂x1
+ε23

∂ε23
∂x1

) ∂ε13
∂x1
=
1

2

∂2w

∂x21
∂ε23
∂x1
=
1

2

(

α+
∂2w

∂x1∂x2

)

(2.9)

and

A=1+
2

3
(1+ν)

m

1−m
B=
2

3
(1+ν)ε0 (2.10)

In the samemanner, we obtain the derivative of ε
p
23 (2.6)2 with respect to x2.

The boundary conditions for the described above approach and the considered in the paper
rectangular boundaryΓ parallel to one of the coordinate axes can be presented in the following
form

∂w

∂n
=α(n1x2−n2x1) on Γ (2.11)

where n1 and n2 are the direction cosines of the external normal to the boundary with respect
to the x1 and x2 axes, respectively.

3. PIES for the elasto-plastic bar torsion problem

3.1. Modeling of the boundary and the domain

Developed by the authors PIES is an analytical modification of the boundary integral equ-
ation (BIE). This modification involves direct inclusion of the boundary shape defined in a
general way in the formalism of the equation. Such an approach allows separation of the appro-
ximation of the boundary shape from the approximation of boundary functions, and thusmakes
it possible to independently increase the accuracy and efficiency of these two elements. From
the very beginning, the main advantage of the PIES method is the elimination of the classical
discretization by elements understood in the manner known from FEM or BEM methods. To
create the shape of the boundary we use any parametric curves, while in the case of modeling
also the domain (as in the present problem) surface patches are used. Figure 2b presents an

Fig. 2. Modelling of the boundary: (a) in BEMby boundary elements using nodes •, (b) in PIES by
four Bézier curves using control points •, and the domain: (c) in BEMby cells using nodes • and ◦,

(d) in PIES by Bézier bilinear surface using control points •

exemplary way of boundary modeling by curves of the first degree, while the modeling of the
domain using the bilinear surface is presented in Fig. 2d. Considering curvilinear shapes, we
have to use curves and patches of higher degrees. Both curves and patches are well-known and



The effective strategy for calculating strains and stresses... 183

efficient computer graphics tools (Foley et al., 1994), the advantages of which are used in PIES.
It shouldbe noted, however, that these curves or patches are not used inPIESas a different kind
of elements known from the element methods. The modeling of the shape in PIES should be
treated as global, because the number of used curves or patches depends only on the complexity
of the shape. Considering, for example, any quadrangle, the definition of the domain requires
only one surface. The type of the surface depends on the type of the boundary that surrounds
the domain. In the case of the quadrangle, we use the bilinear surface, while the area with the
curvilinear boundary is defined efficiently using the bicubic patch. Dealing with more complex
shapes (e.g. polygons) we should join the surfaces. An example of themodeling of the boundary
and the domain in PIES is presented in Fig. 2b,d. For comparison Fig. 2a,c presents the way of
modeling in the BEMmethod.

3.2. Fundamentals of the PIES method

As mentioned earlier, PIES for various differential equations have been obtained by the
analytical modification of BIE, like it is shown in (Zieniuk, 2001, 2013; Zieniuk and Bołtuć,
2006a,b). PIES applied to solve the elasto-plastic torsion of a bar (as multiple solution of the
Poisson equation) is obtained in the same way, and for a smooth boundary can be presented as

1

2
ul(s1)=

n
∑

j=1

sj
∫

sj−1

[

U
∗
lj(s1,s)pj(s)−P

∗
lj(s1,s)j(s)

]

Jj(s) ds+

∫

Ω

U
∗

l (s1,y)F(y) dΩ(y) (3.1)

where y= [y1,y2]∈Ω, Jj(s) is the Jacobian of the transformation

Jj(s)=

√

√

√

√

(∂Γ
(1)
j (s)

∂s

)2
+
(∂Γ

(2)
j (s)

∂s

)2

l=1,2, . . . ,nandn is thenumberof segments thatmake theboundary.Theparametric variables
s,s1 satisfy the following relations sl−1 ¬ s1 ¬ sl, sj−1 ¬ s ¬ sj, where sl−1 and sj−1 are the
beginnings of the segments l,j, respectively, while sl,sj are their ends.
Using equation (3.1) we can solve equation (2.4) and obtain solutions on the boundary.

Therefore, ul(s1) (l = 1,2, . . . ,n) in (3.1) corresponds to the function w(x1,x2) in (2.4) for
x1,x2 ∈Γ .
Thefirst and second integrands in (3.1) are theboundary fundamental solution andboundary

singular solution to the Laplace equation. They can be presented as (Zieniuk, 2001)

U
∗
lj(s1,s)=

1

2π
ln

1
√

η21 +η
2
2

P
∗
lj(s1,s)=

1

2π

η1n1(s)+η2n2(s)

η21 +η
2
2

(3.2)

where η1 = Γ
(1)
l
(s1)− Γ

(1)
j (s) and η2 = Γ

(2)
l
(s1)− Γ

(2)
j (s), while n1(s) and n2(s) are the

direction cosines of the external normal to the boundary. The shape of the boundary defined
by any parametric curves Γ(s) = [Γ(1)(s),Γ(2)(s)]T is included into mathematical formalism of
formulas (3.2). The types of curves and formulas describing them can be found by Foley et al.
(1994).
The integrand in the integral over the domain is presented by the following expression, also

known as the Poisson equation

U
∗

l (s1,y) =
1

2π
ln

1
√

η21+η
2
2

(3.3)

where η1 =Γ
(1)
l
(s1)−y1, η2 =Γ

(2)
l
(s1)−y2.



184 E. Zieniuk, A. Bołtuć

In order to obtain the solutions u in the domainΩ (at pointx= [x1,x2]) we use the integral
identity

u(x) =
n
∑

j=1

sj
∫

sj−1

[

Û
∗

j (x,s)pj(s)− P̂
∗

j (x,s)uj(s)
]

Jj(s) ds+

∫

Ω

ˆ
U
∗

(x,y)F(y) dΩ(y) (3.4)

where the integrands are presented by

Û
∗

j (x,s)=
1

2π
ln

1
√

r̈21 + r̈
2
2

P̂
∗

j (x,s)=
1

2π

r̈1n1(s)+ r̈2n2(s)

r̈21 + r̈
2
2

ˆ
U
∗

(x,y) =
1

2π
ln

1
√

ˆ̈r21 +
ˆ̈r22

(3.5)

and r̈1 = x1−Γ
(1)
j (s), r̈2 = x2−Γ

(2)
j (s),

ˆ̈r1 = x1−y1, ˆ̈r2 = x2−y2, x,y∈Ω and x 6=y. The
form of the Jacobian Jj(s) is the same as in formula (3.1).
Using equation (3.4) we can obtain the solution to equation (2.4) in the domain. Therefore,

u(x) in (3.4) corresponds to the functionw(x1,x2) for x1,x2 ∈Ω.

3.3. Numerical solution of PIES

The solution of PIES represented by formula (3.1) consists of finding unknown boundary
functions uj(s) or pj(s) on particular segments of the boundary. These functions are approxi-
mated by the approximating series

uj(s)=
N−1
∑

k=0

u
(k)
j T

(k)
j (s) pj(s)=

N−1
∑

k=0

p
(k)
j T

(k)
j (s) j=1, . . . ,n (3.6)

where u
(k)
j , p

(k)
j are unknown coefficients, T

(k)
j (s) are the Chebyshev polynomials of the first

kind,N is the number of terms in the series on the segment j and n is the number of segments.
After imposing the boundary conditions, we can obtain the unknown boundary functions for
the individual segments using the approximating form of PIES. This form is obtained after
substituting series (3.6) to analytical form of PIES (3.1)

1

2

N−1
∑

k=0

u
(k)
l
T
(k)
l
(s1)=

n
∑

j=1

N−1
∑

k=0

[

p
(k)
j

sj
∫

sj−1

U
∗
lj(s1,s)−u

(k)
j

sj
∫

sj−1

P
∗
lj(s1,s)

]

Jj(s)T
(k)
j (s) ds

+

∫

Ω

U
∗

l (s1,y)F(y) dΩ(y)

(3.7)

Writing down equation (3.7) at the collocation points si1 (i = 1,2, . . . ,I), I = n×N, we
obtain the system of algebraic equations

Hu=Gp+D (3.8)

whereH= [H]I×I,G= [G]I×I are the square matrices of elements expressed by integrals over
boundary from (3.7), u andp are vectors which contain the coefficients of approximation series
(3.6), whileD is the vector of elements expressed by the integral over the domain from equation
(3.7).
After imposing the boundary conditions and some transformations, equation (3.8) takes the

following form

AX=B+D (3.9)



The effective strategy for calculating strains and stresses... 185

where the vector X contains unknown coefficients from (3.7), while the vector B is known and
depends on the given boundary conditions.
The only problem is the vectorDwhose components require calculation of the functionF(y)

in order to make the integration over the domainΩ. The value of this function is unknown and
depends on the plastic strains which are a part of the solution. For that reason it is necessary
to apply an iterative process and at its first step to solve equation (3.9) with an arbitrary
distribution of the function F(y). We have assumed that the function at the first step of the
iterative process (according to Fig. 3) takes a constant value (zero) in the entire domain.

Fig. 3. The diagram of solving the elasto-plastic torsion of the bar

Finally, the solution to equation (3.9) also requires calculation of the integral over thedomain
which containsF(y). Thementioned at the beginning of this Sectionway of themodelling gives
the possibility for integrating globally over the entire domain without dividing it into cells as it
is done inBEM(Aliabadi, 2002;Ameen, 2005). InBEM, the domain is divided into sub-domains
of elementary shapes over which integrals are calculated with a small number of weights in the
quadrature of numerical integration. Then the values of the so determined integrals are summed.
In PIES, the integrals are calculated globally over the entire domain using the Gauss-Legendre



186 E. Zieniuk, A. Bołtuć

quadrature with a large number of weights. The results of the global integration in PIES on
the example of elasticity problems with body forces can be found among others by Bołtuć and
Zieniuk (2011).
The conceptual diagram of the process of solving the elasto-plastic torsion problem is pre-

sented in Fig. 3.
After solving (3.9) with the adopted initial value of F(y) we have to compute the warping

function at the points of the domain. This is done using integral identity (3.4) after substituting
series (3.6) with known or calculated boundary functions. Then the total strains defined by
formulas (2.3) are calculated. However, it is clear from formulas (2.3) that the calculations do
not require the value of the function u at selected points of the area, but rather the value of
its first and second derivative with respect to x1 and x2. These derivatives can be determined
by analytical differentiation of integral identity (3.4). However, such an approach will result in
multiple calculations of integrals from identity (3.4) and its derivatives, which is particularly
troublesome in an iterative process. Therefore, it has been decided to obtain these derivatives
in a more efficient way, using the proposed and described in the next Section approximation
strategy. This strategy is also included in the diagram presented in Fig. 3.
Having calculated total strains (2.3) we should determine the plastic strains using formulas

(2.6). Then, in order to obtain a better approximation of F(y), we determine its value again
using expression (2.5). Again, we solve equation (3.9) and the solution becomes approximated
in successive iteration steps until fulfilling the given stop criterion. The iterative process should
be recognized as finished if the difference between two lastly obtained values at all considered
points is smaller than the convergence criterion δ.

4. The approximation strategy for derivatives

The proposed in the paper approximation algorithm initially has been based on the already
known fromPIES series used in the boundary function approximation, and in fact has been their
generalization to a series dependent on two variables. They are easily differentiable, hence there
occurs the possibility of direct calculation of the solution derivatives of any order with respect
to any variable. It should be emphasized that the derivatives are obtained in a continuous way
at any point of the area and the boundary. Tests of the approach based on the approximation
series considering the example of linear elasticity problems are presented byBołtuć and Zieniuk
(2013).
Carrying out the research, we noticed a few shortcomings. The most important is ill-

conditioning of the system of algebraic equations solved during approximation. The explicit
form of this system is presented in the paper (Bołtuć and Zieniuk, 2013). Therefore, various im-
provements have been introduced, which are described in detail in (Bołtuć and Zieniuk, 2013).
Themost important is application of the 2D Lagrange polynomial to the approximation of the
derivative of the solutions. Using it, we can eliminate the necessity of solving the ill-conditioned
system of equations. This is particularly important in problems solved iteratively, because ill-
conditioning does not guarantee the stability of solutions for even small disorders of the vector
with constant terms (the right-hand side vector).
In this paper,we focus only on general describing the proposed technique.Details and results

of the research on its efficiency and accuracy can be found in (Bołtuć and Zieniuk, 2013). The
steps of the algorithm are as follows (they are also presented in Fig. 3):

K1 – Assume the number (P ×R) of the interpolation nodes.

K2 – Assume the distribution of interpolation nodes; different node distributions in the square
(P = R = 9) are presented in Fig. 4 and their impact on the quality of the obtained
solutions is described in details in (Bołtuć and Zieniuk, 2013). Figure 4 also shows the



The effective strategy for calculating strains and stresses... 187

Fig. 4. Variants of the arrangement of interpolation nodes: (a) uniform, (b) at roots of the Chebyshev
polynomial of the first kind, (c) at roots of the Legendre polynomial, (d) at roots of the Chebyshev

polynomial of the second kind

distance of the extreme interpolation node fromtheboundaryof theunit square for various
node distributions.

K3 – Map the interpolation nodes from the area in which they are arranged (the unit square)
to the actual area usingmathematical expressions describing surfaces (Foley et al., 1994).
Such a situation is possible because of the fact that the mentioned unit area is also the
domain of polynomials describing the surfaces. An example of mapping the interpolation
nodes from the unit square to the actual shape of the considered area is shown in Fig. 5.
This is an original, efficient andautomaticway of optimal arrangement of the interpolation
nodes with the help of surface patches.

Fig. 5. Mapping the interpolation nodes from the domain of the surface to the actual area

K4 – Obtain solutions u(x) at the interpolation nodes using the integral identity (3.4).

K5 –Create the 2D Lagrange polynomial, which is a generalization of the 1D case and has the
following form

Wu(x1,x2)=
P−1
∑

p=0

R−1
∑

r=0

u(x1p,x2r)Lpr(x1,x2) (4.1)

where

Lpr(x1,x2)=Lp(x1)Lr(x2)

Lp(x1)=
P−1
∏

g=0,g 6=p

x1−x1g
x1p−x1g

Lr(x2)=
R−1
∏

g=0,g 6=r

x2−x2g
x2r −x2g

and u(x1p,x2r) are values of the warping function at P ×R given (declared in step K2)
points of the area and the boundary; the polynomialWu(x1,x2) approximates functionw
from (2.3) and (2.4).



188 E. Zieniuk, A. Bołtuć

K6 – Direct differentiation of the resulting polynomial with respect to each variable in order
to obtain required partial derivatives. As a result, we obtained polynomials which appro-
ximate the derivatives of solutions, e.g.

∂Wu(x1,x2)

∂x1
=
P−1
∑

p=0

R−1
∑

r=0

u(x1p,x2r)
∂Lpr(x1,x2)

∂x1

∂Lpr(x1,x2)

∂x1
=
∂Lp(x1)

∂x1
Lr(x2)

∂Lp(x1)

∂x1
=

P−1
∑

h=1,h 6=p

1

x1p−x1h

P−1
∏

z=1,z 6=p,z 6=h

x1−x1z
x1p−x1z

∂2Wu(x1,x2)

∂x21
=
P−1
∑

p=0

R−1
∑

r=0

u(x1p,x2r)
∂2Lpr(x1,x2)

∂x21

∂2Lpr(x1,x2)

∂x21
=
∂2Lp(x1)

∂x21
Lr(x2)

∂2Lp(x1)

∂x21
=

P−1
∑

h=1,h 6=p

1

x1p−x1h

P−1
∑

o=1,o 6=h,
o 6=p

1

x1p−x1o

P−1
∏

z=1,z 6=p,
z 6=h,z 6=o

x1−x1z
x1p−x1z

(4.2)

The derivatives with respect to the remaining variables are obtained in an analogous way.
As shown in Fig. 3, the differentiated polynomials are used further to determine the warping
function derivatives, and these derivatives to calculate the derivatives of plastic strains and
finally the function F(y) as a result.

5. Analysis of results

We consider the problem of the elasto-plastic torsion of a bar of square cross-section with side
length 2a. Due to double symmetry, only a quarter of the domain has been taken into account
(Fig. 6). The shape of the considered square domain is defined using one bilinear Bézier surface.
This requires only four corner points P0,P1,P2,P3 of the surface.

Fig. 6. The quarter of cross-section of the bar

In calculations, the following values have been assumed: a = 1, ν = 0.3, whilst the dimen-
sionless angle of twist per unit lengthβ=αa/ε0 is increased in steps of one fromβ=1 toβ=6.
The strain hardening parameter m is taken as 0 (perfect plasticity), 0.1, 0.2 and 1 (elasticity).
Initially, we considered torsion of an elastic bar, and therefore β = 1 and m = 1 have

been adopted. This problem has an analytical solution for τ1 = τ13/(2Gε0) (Sokolnikoff, 1956),



The effective strategy for calculating strains and stresses... 189

which can be calculated also numerically using (2.2). Equations (2.3) contain the first derivative
of w with respect to x1, so to determine the shear stresses we could use the proposed strategy
for approximation of the derivative. In order to use it, we have chosen 64 interpolation nodes
(P =R = 8) distributed at places corresponding to the roots of 8th Chebyshev polynomial of
the second kind (Fig. 7).

Fig. 7. The arrangement of the interpolation nodes (P =R=8)

Taking into account the values of w (obtained by 3.4) at these points as well as their coor-
dinates, we can determine Lagrange polynomial (4.1) and, analogously to (4.2), differentiate it.
Then these derivatives are used to determine the stresses presented by (2.2)1. The obtained in
the describedway results at selected points of the domain and compared to analytical solutions
are presented in Table 1.

Table 1.Comparison of the exact elastic τ1 with the obtained by the proposed strategy

(x1,x2) Exact PIES
Relative
error [%]

(0.2,0.2) 0.0971 0.097081 0.01957

(0.4,0.2) 0.0839 0.083950 0.05959

(0.6,0.2) 0.0623 0.062312 0.01926

(0.8,0.2) 0.0333 0.033380 0.24024

(0.2,0.4) 0.2030 0.202945 0.02709

(0.4,0.4) 0.1760 0.176557 0.31647

(0.6,0.4) 0.1320 0.132236 0.17878

(0.8,0.4) 0.0714 0.071442 0.05882

(0.6,0.6) 0.2190 0.219405 0.18493

(0.6,0.8) 0.3380 0.338050 0.01479

(0.8,0.6) 0.1210 0.121260 0.21487

(0.8,0.8) 0.1980 0.197739 0.13181

As shown in Table 1, the solutions obtained by means of PIES and the proposed algorithm
of approximation are characterized by a very high accuracy. A very important advantage is also
the fact that if we want to obtain the value of the derivative of another order or with respect
to another variable, we can use the same Lagrange polynomial, just differentiate it according to
current needs. This avoids repeated calculation of the integrals from analytically differentiated
integral identity (3.4).

Another very interesting and useful advantage of the strategy is the ability to determine the
derivatives of solutions on the boundary, again by the previously obtainedLagrange polynomial.
In the case of the classical approach, obtaining stresses on the boundary is possible using the
derivative of the integral identity, wherein there is a singularity arising fromthe fact thatx→Γ .
Using the proposed technique, derivatives of any order with respect to any variable at any point
of the boundary and area are calculated using the same expression. The stress values at selected



190 E. Zieniuk, A. Bołtuć

points along the boundary compared with the exact solutions (Sokolnikoff, 1956) are presented
in Table 2.

Table 2.Values of τ1 obtained on the boundary

(x1,x2) Exact PIES
Relative
error [%]

(0.0,0.2) 0.101 0.101450 0.44554

(0.0,0.4) 0.212 0.211629 0.17500

(0.0,0.6) 0.339 0.339247 0.07286

(0.0,0.8) 0.492 0.492385 0.07825

(0.0,1.0) 0.675 0.674900 0.01481

(0.2,1.0) 0.658 0.658401 0.06094

(0.4,1.0) 0.605 0.605372 0.06148

(0.6,1.0) 0.507 0.506964 0.00710

(0.8,1.0) 0.342 0.342750 0.21929

As can be seen inTable 2, solutions obtained by the approximation strategy are very similar
to the analytical results, relative errors do not exceed 0.5%. In order to confirm the reliability
and accuracy of the approach, we also examined the dimensionless momentM∗ =M/(2Gε0a

3),
where M =

∫∫

Ω
(x1τ23 − x2τ13)dx1dx2 and the dimensionless maximum shear stress τ∗max =

τmax/(2Gε0), where τmax =(
√

τ213+ τ
2
23)max. The obtained byPIES results have been compared

with the exact solutions and are presented in Table 3.

Table 3.Comparison of M∗ and τ∗max obtained by PIES with the exact results

Exact PIES
Relative
error [%]

M∗ 1.125 1.12463 0.0328

τ∗max 0.6754 0.67490 0.0740

Table 3 presents the high accuracy of the proposed approach.
The last test carried out for the elastic case concerned the comparison of the dimensionless

moment M∗ for m = 1 and different values of β. These results, together with the analytical
solutions and obtained by the methods: boundary integral method (BIM) (Mendelson, 1975)
andMFS (Kołodziej and Gorzelańczyk, 2012) are presented in Fig. 8.

Fig. 8. Comparison ofM∗ for various values of β

As can be seen in Fig. 8, the values obtained using BIM and PIES are the same as the
analytical solutions. The only slight divergence can be observed in the case of theMFSmethod
for β greater than 3.



The effective strategy for calculating strains and stresses... 191

Since the approximation strategy has been accurate in the case of the elastic problem, it has
been decided to examine its possibilities considering the plastic problem. The solution of such
problems requires calculating of the integrals present in (3.1), (3.4)with integrandswhich require
derivatives of plastic strains (2.5) at many points of the domain. Furthermore, the problem is
resolved in an iterative process, which involves repeated calculation of the integrals included in
the derivative of expression (3.4). Therefore, it has been decided to apply the approximation
strategy, which can be done in two ways:

• make the approximation of expressions ε13,ε23 (2.3) on the basis of the function u, using
them calculate plastic strains ε

p
13,ε

p
23 (2.6) and then approximate the derivatives from

(2.5) – two interpolation polynomials are generated; differentiation is made to obtain the
first-order derivatives only,

• make the approximation of expressions ε13,ε23 (2.3) andderivatives from(2.5) on the basis
of the function u, expressions (2.7)-(2.10) and similarly generated for the derivative ε

p
23

with respect to x2. One interpolation polynomial is generated; differentiation is made to
obtain the derivatives of the first and second order.

In the paper, the second approach is chosen and applied (it is also presented in Fig. 3)
because of its efficiency. It is characterized by the fact that at each step of the iterative process,
the Lagrange polynomial is obtained only once on the basis of the function u obtained at the
interpolation points. Then, only the derivatives of the polynomial with respect to the proper
variable and of the proper order are generated.

In order to perform tests, the dimensionless angle of twist per unit length β has been then
increased in unit steps from 1 to 6 for different values of the strain hardening parameter m.
Initially, the dimensionless maximum shear stresses τ∗max for m= 0.2 have been obtained. For
this purpose, we used PIES and the strategy for the approximation with 196 interpolation
nodes (P = R = 14) placed at roots of the 14th Chebyshev polynomial of the second kind.
The obtained results have been compared with solutions from other numerical methods: BIM
(Mendelson, 1975) and FDM (Mendelson, 1968) as shown in Table 4.

Table 4.Comparison of τ∗max obtained by PIES with other numerical methods form=0.2

β PIES BIM FDM

2 0.87558 0.879 0.881

3 1.01375 1.020 1.022

4 1.14824 1.150 1.156

5 1.28215 1.280 1.290

6 1.4174 1.400 1.426

As shown in Table 4, the results obtained by PIES using the approximation strategy are
very accurate comparing with other numerical methods.

Thenext examined parameter has been the dimensionlessmomentM∗ form=0.0, 0.1. This
time we compare results with two other numerical methods: BIM (Mendelson, 1975) and MFS
(Kołodziej and Gorzelańczyk, 2012). The parameters of PIES and the approximation strategy
are taken from the tests presented above. The results are shown in Fig. 9.

Again, the results obtained using PIES are very similar to those obtained by Mendelson
(1975) andKolodziej andGorzelańczyk (2012). It shouldbenoted that in theproposedapproach,
the discretization has not been applied. The whole bar cross-section has been modelled by one
surface. Therefore, it has not been necessary to calculate and sum the integrals over cells, the
global integration over the entire areawith a larger numberofweights hasbeenused instead.The
necessity of calculation of integrals from thedifferentiated integral identity for displacements has



192 E. Zieniuk, A. Bołtuć

Fig. 9. Comparison ofM∗ for various variants of β andm

been eliminated and replaced by the numerical approximation of the derivatives of the solution,
which is an elementary procedure.

The final stage of the comparative analysis is to examine the effectiveness of the proposed
method measured by the number of required iterations. Again, the data are taken from the
literature (Mendelson, 1975), therefore, we can compare only two cases:

a) m=0, β=5with δ=0.0001,

b) m=0, β=6with δ=0.0001 and δ=0.00001.

Table 5 presents the number of iterations required to obtain the solution in case a) with the
given level of the accuracy, taking into account three methods: PIES, BIM (Mendelson, 1975)
and FDM (Mendelson, 1975). Table 6 contains the number of iterations required for obtaining
the solution in case b) taking into account two different levels of the accuracy and twomethods:
PIES and BIM (Mendelson, 1975). No results are available for theMFS method.

Table 5.The number of iterations in case a)

FDM BIM PIES

203 33 21

Table 6.The number of iterations in the case b)

BIM PIES

δ=0.0001 39 25

δ=0.00001 53 32

As can be noticed in Table 5 and 6, in the considered cases the proposed method requires
fewer iterations than the rest of the studiedmethods. For completeness, we should also examine
the execution time. However, it is not possible for the considered problem and included in the
manuscript results, because the cited literature does not provide such data. Such tests will be
possible to do, if for solving the considered problemweuse professional software based on known
numerical methods.We will include it in our further research plans.

6. Conclusions

The efficient approximation strategy for the derivatives of solutions has been tested on an exam-
ple of nonlinear torsion of the bar using the PIES method. The application of PIES eliminates
the necessity of dividing the considered domain into cells, as it is done in BEM in order to
calculate the integrals over the domain. The application of surfaces to define shapes gives the
possibility for implementing the global integration with a large number of weights in theGauss-
-Legendre quadrature. Moreover, the proposed in the paper approximation strategy enables us



The effective strategy for calculating strains and stresses... 193

to avoid repeated calculation of the integrals (often singular) and gives the possibility for effec-
tive determination of the derivatives of any order with respect to various variables at any point
of the domain and the boundary. It should be emphasized that the obtained results are highly
accurate, as shown in comparison to analytical solutions and these obtained by other numerical
methods.
We realize that the solved in the paper problem is quite elementary, but we had to test

the concept and the mathematical apparatus based on the PIESmethod. At the next research
stage, the gained knowledgemakes it possible to apply the strategy and its advantages to more
complex nonlinear mechanical problems. In such problems, we will eliminate the necessity of
calculating strongly singular integrals, which are present in the integral identity for stresses.

References

1. AliabadiM.H., 2002,The Boundary ElementMethod, Volume 2,Applications in Solid and Struc-
tures, JohnWiley & Sons, Ltd., Chichester

2. Ameen M., 2005,Computational Elasticity, Alpha Science International Ltd., Harrow

3. Bołtuć A., Zieniuk E., 2011, PIES in problems of 2D elasticity with body forces on polygonal
domains, Journal of Theoretical and Applied Mechanics, 49, 2, 369-384

4. Bołtuć A., Zieniuk E., 2013, The effective determination of stress by PIES method using the
generalized strategy of derivative approximation (in Polish), Modelowanie Inzynierskie, 16, 47,
41-47

5. Foley J., van Dam A., Feiner S., Hughes J., Phillips R., 1994, Introduction to Computer
Graphics, Addison-Wesley

6. Kołodziej J.A., Gorzelańczyk P., 2012, Application of method of fundamental solutions for
elasto-plastic torsion of prismatic rods,Engineering Analysis with Boundary Element, 36, 81-86

7. LiuG.R., 2002,MeshfreeMethods –Moving Beyond the Finite ElementMethod, CRCPress,Boca
Raton

8. Mamedov A., 1995, An inverse problem related to the determination of elastoplastic properties
of a cylindrical bar, International Journal of Non-Linear Mechanics, 30, 23-32

9. Mendelson A., 1968, Elastic-Plastic Torsion Problems for Strain-Hardening Materials, NASA
TND-4391

10. Mendelson A., 1975, Solution of Elastoplastic Torsion Problem by Boundary Integral Method,
NASA TND-7872

11. Sapountzakis E.J., Tsipiras V.J., 2009, Nonlinear inelastic uniform torsion of composite bars
by BEM,Computers and Structures, 87, 3/4, 151-166

12. Sokolnikoff I.S., 1956,Mathematical Theory of Elasticity, McGraw-Hill Book Co. Inc.

13. Sokolowski W. W., 1957,Theory of Plasticity (in polish), PWN,Warszawa

14. YamadaY., Nakagiri S., TakatsukaK., 1972, Elastic-plastic analysis of Saint-Venant torsion
problemby a hybrid stressmodel, International Journal for Numerical Methods in Engineering, 5,
2, 193-207

15. ZieniukE., 2001,Potential problemswith polygonal boundaries by aBEMwith parametric linear
functions,Engineering Analysis with Boundary Element, 25, 185-190

16. Zieniuk E., 2013, Computational Method PIES for Solving Boundary Value Problems, PWN,
Warsaw

17. Zieniuk E., Bołtuć A., 2006a, Bézier curves in the modeling of boundary geometry for 2D
boundary problems defined by Helmholtz equation, Journal of Computational Acoustics, 14, 3,
353-367



194 E. Zieniuk, A. Bołtuć

18. Zieniuk E., Bołtuć A., 2006b, Non-element method of solving 2D boundary problems defined
on polygonal domainsmodeled byNavier equation, International Journal of Solids and Structures,
43, 7939-7958

19. Zienkiewicz O., 1977,The Finite Element Methods, McGraw-Hill, London

Manuscript received April 16, 2014; accepted for print August 19, 2014