- Email: [email protected]
Automatic mesh generation for computing symmetric local stresses by 2d-viscoplastic BEM with modified transfinite method
Engineering Fracture Mechanics Vol. 52, No. 6, pp. 1139--1150, 1995
Pergamon
0013-7944(94)00287-8
Copyright © 1995 Elsevier Science Ltd Printed in Great Britain. All rights reserved 0013-7944/95 $9.50 +0.00
AUTOMATIC MESH GENERATION FOR COMPUTING SYMMETRIC LOCAL STRESSES BY 2D-VISCOPLASTIC BEM WITH MODIFIED TRANSFINITE METHOD YONG LIU Dept of Mechanical Engng, Zhejiang University of Technology, Hangzhou, ER. China and HEINZ ANTES Institute of Applied Mechanics, TU Braunschweig, Germany Abstract--This paper presents a scheme of automatic mesh generation for computing the local stresses by 2D-viscoplastic boundary element method (BEM) with modified transfinite mapping method. The advanced techniques for mesh spacing and region of transfinite mapping method in BEM are studied. The applications demonstrate that the method in BEM gives high computational efficiency and precise modelling of boundaries. This reduces the time and effort required of the analyst to set up a nonlinear model, the amount of input data required is reduced greatly. The numerical examples are given in a later section.
1. INTRODUCTION IN COMPUTINGlocal stresses by 2D-viscoplastic BEM, numerical automatic mesh generation method is an important aspect of 2D-viscoplastic problem of BEM. A well distributed grid is favourable for obtaining an accurate numerical solution with higher convergence rate and less computer storage, and a minimum of input data including adjustable parameters. Therefore, considerable effort has been spent on the development of numerical procedures for constructing suitable computational meshes [1, 2]. The transfinite mapping techniques are a class of methods for establishing curvilinear co-ordinate system in arbitrary domains. The method was primarily developed by Gordon and Hall [3,4] and Cook [5] for the approximation of complex surfaces and volumes in FEM problems. But all the above previous work of transfinite mapping method in FEM can not be used directly in the problem of nonlinear BEM because the mesh of BEM requires the line cells in the boundary for the elastic region and internal 2D-cells for the nonlinear region. In this paper, the automatic meshing schemes for computing symmetric local stresses in the 2D-viscoplastic BEM problem are considered by using the transfinite mapping method. The techniques for the mesh spacing and region controlling in the grid of BEM are studied. The programs for the automatic grid generation in BEM have been developed. The numerical examples demonstrate that the method of modified automatic grid generation in BEM is highly effective and accurate. This makes the input data procedures simple to implement and the efficiency for computing the local stress is higher. 2. THE MODIFIED TRANSFINITE METHOD FOR MESH GENERATION
The general transfinite method describes an approximate boundary or region which will match a desired, or true, boundary or region at a non-denumberable number of points. It is this property which gives rise to the term transfinite mapping [11]. For the cases of two dimensions, let us introduce a bilinear projector. The bilinear projector which operates on a region F boundary by four curves u(O,t), u(l, t), u(s, 0) and u(s,1) is shown in Fig. 1(b). Two basic projectors can be formed, one interpolating in the s direction and one in the t direction.
~l [F]-p~(s,t )=(1-t )u(s, O) +
tu(s, 1), 0~
~2[F] =-p2(s, t) = (1 - s)u(O, t) + su(1,t), 1139
0 ~< t ~< 1.
(1)
1140
YONG LIU Ibl
V.,
u[0,1]~
uP.ll
p.q
utO.q~ utO,OI
ale.hi
utl,Ol
$
)x [UI
P.01 Fig. 1. (a) The mesh after mapping, (b) the actual mesh region.
The bilinear projector is defined as shown in Fig. l(b).
{ (~,~,)[V] _={~, [Vl + ~'~[rl - ~, [~} = e.(s, t) Pn(s, t) = (1 - t)u(s, O) + tu(s, 1) + (1 - s)u(O, t) + su(1, t) - (1 - s)(1 - t) x u(O, O) - (1 - s)tu(O, 1) - s(1 - t)u(1, O) - stu(1, 1)
(2)
where u(0, t), u(1, t), u(s, 0), u(s, 1) are the boundary functions of the region F. u(0, 0), u(0, 1), u(1, 0), u(1, 1) are the corner points at the region. This projector represents a curvilinear coordinate system created by a mapping of the unit square on to F. It may be called the transfinite bilinear Lagrange interplant of F; the transfinite mappings may be used to map any mesh in any mesh in the appropriate primitive polygon to the actual region. Thus, any mesh topology may be used. However, it is generally convenient to sacrifice topological generality by choosing the intersections of the constant coordinate curves for use as node points in order to minimize user input. The bilinear projectors thus create a natural partition of the region composed of quadrilateral elements. This is very useful for the automatic meshing generation of 2D-viscoplastic BEM. But the above eq. (2) has two drawbacks. One is that it cannot give the suitable mesh space for different local stresses. The other is that it cannot contol the size of the nonlinear region of a structure. In order to alleviate the two drawbacks, we present a modified transfinite method. The method is given as follows: Let us define a new projector goM and introduce the mesh spacing controlling parameter p~, i = 1, 2, 3, 4, with the region size controlling parameter ei, i = 1, 2, 3, 4. Thus, we have ~M(s, t) = (1 -- t)u(s p~, O)el + tu(s p:, 1)e2 + (1 -- s)u(O, tP3)e3 + su(l, tP~)e, -- (1 -- s) x (1 - t)u(O, O) - (1 - s)tu(O, 1) - s(1 - t)u(1, 0) - stu(1, 1).
(3)
With Pi and ei, we can control the mesh spacing and mesh region. Equation (3) is the modified transfinite method. Three basic equations for 2D-viscoplastic BEM.
(a) The boundary integral equations The starting integral equation for the boundary element method applied to initial strain problem is as follows
Automatic mesh generation for computing symmetric local stresses
1141
where F represents the boundary of the region and f~ is interior. The tensors u*, p~ and a ~ corresponding to the fundamental solutions have been presented in previous publications [7, 8] and the expression for g~j is the following form: gi/=
G 4(1 - v ) [ 2 ~ + ( 1 -4v)iT,~6~j],
for 2D plane strain
(6)
where for plane stress v is replaced by ~, E~ is the initial strain. The application of eq. (4) in discretized form to all boundary nodes generates the following matrix relationship: Hfi = Gp + Qia
(7)
in which matrix Q stands for the initial stress integral. Similarly, the stress rates at the internal points are given by eq. (5) in discretized form # = G'0 - H'fi + Q,~a
(8)
where matrix Q' represents the initial stress integrals plus the independent terms f0. Matrices H, G, H' and G' are the same as those obtained for pure elastic analysis [7]. Taking into consideration that for a well-posed problem a sufficient number of tractions and boundary displacement is prescribed, eqs (7) and (8) can be further written as A~ = f + Qi ~
(9)
6 = - A ' ~ + [" + Q'/~
(10)
and
where vector ~ is formed by the unknown tractions and boundary displacements and the contribution of the prescribed values is included in vectors [r and i". From the multiplication of eq. (9) by A-' comes
{
x
= R ~ a --I- lil
R =
A-IQ,
ill =
A-if.
(1 1)
Substituting eq. (11) into eq. (10) yields d = Vi ~ + fi where
V = Q' -
A'R,
(12)
fi = i" - A'fit.
Note that the elastic solution to the rate problem is given by the vectors ~h and ~. From the above it is seen that eq. (12) represents a system of ordinary differential equations for stresses at selected boundary nodes and internal points which can be solved by standard methods, producing a unique solution for rate dependent problem. (b ) Viscoplasticity relations
It is known that one explicit model which has wide applicability for isotropic material is offered by the following viscoplastic flow rule [6, 10]. i s = i p = ~, < 4 , ( f )
>
~-~. dQ
(13)
We restrict ourselves to associated plasticity situations, in which case F = Q and a sufficiently general expression for viscoplastic strain rate is given by a power law for ~b(F) such as ~b(F)
= ~ F _ FoT¢ ,
L F0 ]
r-to or ~b(F)=eS2(-~-*) - 1
(14)
where N,, N2 are prescribed numbers, F is the yield function, the onset of viscoplastic behaviour is governed by a scalar yield condition of the form F ( a ) - Fo = 0
(15)
1142
YONG LIU
in which F0 can be the uniaxial yield stress, F0 may itself be a function of an effective strain gp and for a linear strain hardening response, we have F0 = try + n ' g " = try + f(ffP)
(16)
where try is the initial uniaxial yield stress, H ' is the viscoplastic strain hardening parameter, ¢P is the current effective viscoplastic strain and f(¢P) can be determined by experiment. (c ) The computational algorithms The initial distribution of internal stress and boundary traction and displacements is obtained from the solution of the corresponding elastic problem using eqs (11) and (12) where i s = i p = O. The load increments are applied as discrete steps and thus Af. = 0, Af~ = 0 for all time steps other than the first within an increment. (i) Starting from known values of X,, a,, d at a time and compute the rate of viscoplastic strain i p by the constitutive relations (13-16). (ii) Calculate the ~P as a~{ = i ~ a t .
(17)
¢,P+t = ¢~ + A ~ .
(18)
(iii) Compute the value of the boundary unknowns AX. = R A ~ + Am
(19)
X,, + 1 = X. + AX,.
(20)
(iv) Determine the internal stresses Aa. = VA~ + An
(21)
O'n + l = O'n "q- Ao'n.
(22)
We arrive at all the starting values of the next step for which the same process can be started in the next interval. In order to improve the accuracy of this explicit time step algorithm which is equivalent to the Euler extrapolation procedure, an inner iterative loop can be imposed on the steps (i-iv) where step (ii) is replaced by calculating A~P.= (iP. + iP.+~)at./2.
(23)
i p. + I = iP
(24)
Taking initially and repeating the process to (iv) until some convergence criterion is reached using the second viscoplastic strain invariant for instance. The iP,+t is defined by using a limited Taylor series expansion such as i t + , = i p. + H.Aa. where
(25)
a[i~] T H.
= - -
a[~.]
=
[H,,(G,)].
The matrices [l-I. ] depend on the stress level o. and can be easily evaluated with the constitutive laws (13). The viscoplastic strain increment can be written as A ~ = At,[(1 - 0)i~ + 0i~+ ~]
(26)
where for 0 = 0 represents the Euler time integration scheme, the case 0 = 1/2 results in the so-called implicit trapezoidal scheme i.e. eq. (23) and 0 = 1 gives the fully implicit scheme, we get A ~ = ~ A t , + 0[I-I,]At, Aa,.
(27)
Automatic mesh generation for computing symmetric local stresses
1143
Using the incremental from (21) and substituting for A~ from (27), then (27) becomes Aa. = Vi.PAt. + 0VH.At.A~r. + An.
(28)
Equation (28) can be solved iteratively for the increment Aa. and the iterations go on until some convergence is achieved using the ratio ] ,]Aa. ,1~ - ,]Act. 1,~- ' ]
ep
(29)
where Ep is the convergent tolerance. Following the method used by Zienkiewicz and Cormeau the magnitude of this time step is controlled by a factor z which limits the maximum effective viscoplastic strain increment Ag~ as a fraction of the total effective strain ~,, so that Ag.p = ~ 2 .p -p At. ~<"re,,
(30)
where AgP. must be computed to satisfy (30) at each internal point and the least value taken for computation. The above time step limiting method is basically empirical and z in the range of 0.01-0.15 is found to be effective for explicit algorithms and for implicit scheme values of T up to 1 have been found to be stable with BEM. A further limit is generally imposed where the change in the time step length between any two intervals is limited according to Atn+ l ~< 1.5At..
(31)
3. NUMERICAL EXAMPLES In order to show the advantages of automatic meshing generation for 2D-viscoplastic BEM by using the transfinite mapping method and to demonstrate its effectiveness, the following examples are used (a) a plate with hole under uniform tension (b) a plate with V-shaped notch under uniform tension (c) a plate with U-shaped notch under uniform tension.
Example (a). A plate with hole under uniform tension The problem considered is a square 10 x 10 mm 2 plate with a circular hole under uniform tension, as shown in Fig. 2 [9]. The radius of the hole was taken as 5.0 mm, owing to symmetry only a quarter of the plate need to be taken for analysis, as shown in Fig. 3. Considering the plate with the hole under uniform tension, the viscoplastic deformation takes place firstly near the hole of the plate, so it needs taking the AFGE as for generating the 2D grids of the viscoplastic region near the hole, Four functions in eq. (3) are taken by the expressions
J
R(1 + e)cos(2 s )
COS
U(Sp2, 1)=
R(1 + e)sin(2 s )
(32)
u(1, tp4)= IR oeRtP 1 where Pt -'-P2 = 1, P3 = P4 = P. e is the parameter of controlling the viscoplastic region, p is grid space regulating parameter. When p > 1, the generated meshing graph is denser near the AE curve. When 0 < p < 1, the generating meshing graph is sparse near the AE curve, as shown in Fig. 4. In order to calculate the actual stress field, the parameter e is taken as e = 0.9 and the meshing space regulating term p = 1.2. The final generating meshing graph is shown as Fig. 5. The problem to consider is a 20 x 36 mm 2 plate with a V-shaped notch under uniform tension, as shown in Fig. 6(a), owing to symmetry only a quarter of the plate needs to be taken for analysis as shown in Fig. 6(b).
1144
Y O N G LIU
I
0 lOmm 4 O"a
!
Q
i - - 4
36 m m Fig. 2. The plate with the hole under uniform tension.
D
C
E
F
I Fig. 3. Quarter of the plate.
Let us consider the plate with a V-shaped notch under uniform tension. The initial viscoplastic region is taken as ABFG for computation. Four functions in eq. (3) are taken by the expression
u(s p2,1)
=~ es'(L--t,)+tle ] LsP(a+tle)+tl--tle (33)
Let"1
I[l[(O,tP3)~-Itlet__.:ltPtll,
u(l, tp,) = F k a + tl_l
(a)
(b) I
i
I
I
I
I
I
I
I
f
I
I
I
I
Fig. 4. The effect of grid space regulating term p. (a) p = 1.3. (b) p = 0.4.
q
I
I
f
I
I
I
I
I
Automatic mesh generation for computing symmetric local stresses
1145
(a)
(b)
Fig. 5. The meshes generated by bilinear projector. (a) Mesh 1. (b) Mesh 2.
where Pm= P2 = P3 = P 4 = P, e is used to define the viscoplastic region, p is the spacing regulating parameter. Its effect is the same as in example (a). The parameter e = 2/3 and p = 1.2, and the generating meshing graph is shown in Fig. 7. The material and load parameters are taken as [10] E = 70,000 MPa, ay = 243 MPa, v = 0.2, h' = 0.0, the tension load tra is taken in 2a,/ay = 1.3 for plane stress. The viscoplastic constants [6]: NI = 1, y = 0.001 day, 0 = 0.0, and 0 = 1.0, • = 0.01. By using the above 2D-viscoplastic BEM program together with the mesh in Fig. 7(b), the calculation results are showed in Figs 8 and 9. Figure 8 shows that as the days increase, the viscoplastic displacement with two computing schemes (0 = 0 and 0 = 1.0) at point C approaches the steady state solution, i.e. the elastoplastic solution. Figure 9 shows the comparison between the results of viscoplastic solution and elastoplastic solution by BEM. As the steady state solution o f viscoplastic BEM is reached, good agreement has been obtained between the two methods. This shows the effectiveness of the method in this paper.
Example (c). A plate with U-shaped notch under tension The problem considered is a plate with U-shaped notch under uniform tension, as shown in Fig. 10(a). Owing to symmetry only a quarter of the domain needs to be taken for analysis as shown in Fig. 10(b). The specified boundary conditions are traction normal to the boundary equal to along CD. Zero displacement normal to the boundary along AB and BC and zero traction otherwise.
1146
Y O N G LIU
(a)
90"
"llt-'--
o'a
(lra
"lit"---
I_
36ram
_
vI
(b)
B G
D
L ).
Fig. 6. (a) A plate with V-shaped notch. (b) Quarter of the plate.
Considering the plate with U-shaped notch under tension, the nonlinear region may be firstly taken place near the U-shaped notch. Hence the region A F G E should be taken as the computational nonlinear region. The four functions of eq. (3) are taken as in eq. (32) according to the actual specimen of Fig. 10(a). The material constants are taken as follows: R = 0.0425 m, plate thickness b = 0.0041 m A F = E G = e . R = 0.06375 m, E = 70 GPa, v = 0.27, ay = 320 MPa, N = 1, v = 0.001/day, 0 = 0.0, = 0.01. By using meshes of Fig. ll(a) and (b) and viscoplastic BEM program[15], the stress concentration factor K at root A in Fig. 10 can be calculated in Table 1 together with results in refs [12-14] Table 1. Comparison of stress concentration factor K Result of
Result of
Result o f
Result of
Result of
Fig. 10(a)
Fig. 10(b)
rcf. [12]
ref. [131
rcf. [141
2.14
2.07
1.94
2.11
2.0
. ~
I
I
I
I
I
(
1 I
I
!
I
I'
I
I
I
\ \ \ \\,
I
I
I
I
Fig. 7. Two meshes generated by bilinear projector of modified transfinite mapping method. (a) Mesh 1. (b) Mesh 2. 38 37.9 37,8 E E
37.7 37.6 37,5
1~ 37,4 ! 37,3 37.2 37,1 37
0
I
I
I
I
2
4
6
8
I
I
I
I
I
I
I
I
I
I
10 12 14 16 18 20 22 24 26 28 30 Days
-e-VPBEM--1 (O = 0.0)
~
VPBEM--2 ~ (O = 1.0)
Elastoplastic
Fig. 8. The displacement at point C vs times. Loading-disp. curves 1.2
-
f./.o
I.C
~" O
/f
0.8j
"~
0.~
~
0.4
"~
0.2
~ 1 ' 5
I 10
I 15
I 20
I 25
30
uc/lxE-4 ~'" Elutoplastic B E M - - - o - - - V i ~ o p l u t i ¢ BEbl Fig. 9. Comparison of the viscoplastic steady state solution with elastoplastic solution.
1148
YONG
LIU
(a) q •
q
D
0
f
q
L I-
im
.I
0.890 m
(b)
~XX~XX ~ X ~ XX
E Fig. 10. (a) and (b).
(a) I
I
I
I
I
I
I
I
I
I
I
I
I
I
(b) l
I
l
l
I
l
I
I
I
Fig. 1 I. (a) and (b).
l
l
I
I
l
Ilr I
~
Automatic mesh generation for computing symmetric local stresses
1149
5. CONCLUSIONS The modified transfinite mapping method is proposed in 2D-nonlinear BEM for the automatic mesh generation of specimens with symmetric notches. Although the computational efficiency of general BEM programs has steadily increased over the past years, the time consuming process of accurate data preparation is still a burden for most BEM users in engineering. The modified transfinite mappings based on boundary curve descriptions provided an effective basis for automatic generation of visco-plastic BEM. They are generally, simple to implement and computationally efficient. Grids produced by this method have predictable forms wherein the interior of the mesh accurately reflects boundary curvature information. The meshing spacing regulating parameter proposed in this paper can control the grid space and mesh size effectively. The numerical examples show that the method in this paper has the advantages of high computational efficiency, precise modelling of boundary minimum input data and is good for computing the local stress of 2D specimens. Therefore, the modified transfinite mapping method can have applications in the problem of automatic mesh generation for nonlinear BEM more and more. It is useful for the 2D-viscoplastic BEM analysis for preliminary research and Alexander yon Humboldt foundation. Acknowledgements--The project was supported by Zhejiang Science Foundation.
REFERENCES [I] S. W. Chae and K. J. Bathe, On automatic mesh construction and mesh refinement in finite element analysis. Comput. Structures 32, (3/4) 911-936 (1989). [2] J. Zhu, A hybrid differential-algebraic method for three-dimensional grid generation. Int. J. numer. Meth. Engng 29, 1271-1279 (1990). [3] W. J. Gordon and C. A. Hall, Construction of curvilinear coordinate system and applications to mesh generation. Int. J. numer. Meth. Engng 7, 461-477 (1973). [4] C. A. Hall, Transfinite interpolation and applications to engineering problems, in Theory of Approximation (Edited by Law and Sahney) pp. 308-331. Academic Press (1976). [5] W. A. Cook, Body oriented (natural) co-ordinates for generating three dimensional meshes. Int. J. numer. Meth. Engng 8, 27-43 (1974). [6] O. C. Zienkiewicz and I. C. Cormeau, Visco-plasticity, plasticity and creep in solids, A unified numerical solution approach, Int. J. Numer. Meth. Engng 8, 821-845 (1974). [7] C. A. Brebbia, J. C. F. Telles and L. C. Wrobel, Boundary Element Techniques: Theory and Applications in Engineering, Springer, Berlin (1983). [8] Liu Yong and Shen Najie, Residual stress analysis for actual material model of autofrettage tube by nonlinear BEM. Int. J. Pres. Ves. Piping 48, 21-35 (1991). [9] P. S. Theocaris and E. Marketos, Elastic-plastic analysis of perforated thin strips of a strain hardening material. J. Mech. Phys. Solids 12, 377-396 (1964). [10] G. C. Nayak and O. C. Zienkiewicz, Elastoplastic stress a generalization for various constitutive relations including generalization for various constitutive relations including strain softening. Int. J. numer. Meth. Engng 5, 1113-1135 (1972). [I 1] R. Haber, M. S. Shephard, J. F. Abol, R. H. GaUagher and D. P. Groenberg, A general two-dimensional, graphical finite element preprocessor utilizing discrete transfinite mappings. Int. J. numer. Meth. Engng 17, 1015-1044 (1981). [12] Yong Liu, Naijie Shen, and Wai Jiang, An effective method for calculating the local elastoplastic stress of a tension specimen with symmetric U shaped grooves. Engng Fracture Mech. 43, 1085-1091 (1992). [13] K. Axelsson and A. Samuelsson, FEM analysis of elastoplastic materials displaying mixed hardening. Int. J. numer. Meth. Engng 14, 211-225 (1979). [14] J. R. Crews, Local plastic stresses in aluminium alloy specimens with stress concentration factor and constantamplitude loading. NASA TN 3813, Washington (1965). [15] Yong Liu and N. J. Shen, An BEM implementation for elastoplastic analysis by using the steady state solutions of viscoplastic problem. The computational mechanics and its applications in engineering. Proe. 3rd-Computational Mechanics Conf. East China. China University of Science & Technology Ltd. pp. 41-46 (1992). (Received 17 February 1993)
EFM 52/6-~L
Pergamon
0013-7944(94)00287-8
Copyright © 1995 Elsevier Science Ltd Printed in Great Britain. All rights reserved 0013-7944/95 $9.50 +0.00
AUTOMATIC MESH GENERATION FOR COMPUTING SYMMETRIC LOCAL STRESSES BY 2D-VISCOPLASTIC BEM WITH MODIFIED TRANSFINITE METHOD YONG LIU Dept of Mechanical Engng, Zhejiang University of Technology, Hangzhou, ER. China and HEINZ ANTES Institute of Applied Mechanics, TU Braunschweig, Germany Abstract--This paper presents a scheme of automatic mesh generation for computing the local stresses by 2D-viscoplastic boundary element method (BEM) with modified transfinite mapping method. The advanced techniques for mesh spacing and region of transfinite mapping method in BEM are studied. The applications demonstrate that the method in BEM gives high computational efficiency and precise modelling of boundaries. This reduces the time and effort required of the analyst to set up a nonlinear model, the amount of input data required is reduced greatly. The numerical examples are given in a later section.
1. INTRODUCTION IN COMPUTINGlocal stresses by 2D-viscoplastic BEM, numerical automatic mesh generation method is an important aspect of 2D-viscoplastic problem of BEM. A well distributed grid is favourable for obtaining an accurate numerical solution with higher convergence rate and less computer storage, and a minimum of input data including adjustable parameters. Therefore, considerable effort has been spent on the development of numerical procedures for constructing suitable computational meshes [1, 2]. The transfinite mapping techniques are a class of methods for establishing curvilinear co-ordinate system in arbitrary domains. The method was primarily developed by Gordon and Hall [3,4] and Cook [5] for the approximation of complex surfaces and volumes in FEM problems. But all the above previous work of transfinite mapping method in FEM can not be used directly in the problem of nonlinear BEM because the mesh of BEM requires the line cells in the boundary for the elastic region and internal 2D-cells for the nonlinear region. In this paper, the automatic meshing schemes for computing symmetric local stresses in the 2D-viscoplastic BEM problem are considered by using the transfinite mapping method. The techniques for the mesh spacing and region controlling in the grid of BEM are studied. The programs for the automatic grid generation in BEM have been developed. The numerical examples demonstrate that the method of modified automatic grid generation in BEM is highly effective and accurate. This makes the input data procedures simple to implement and the efficiency for computing the local stress is higher. 2. THE MODIFIED TRANSFINITE METHOD FOR MESH GENERATION
The general transfinite method describes an approximate boundary or region which will match a desired, or true, boundary or region at a non-denumberable number of points. It is this property which gives rise to the term transfinite mapping [11]. For the cases of two dimensions, let us introduce a bilinear projector. The bilinear projector which operates on a region F boundary by four curves u(O,t), u(l, t), u(s, 0) and u(s,1) is shown in Fig. 1(b). Two basic projectors can be formed, one interpolating in the s direction and one in the t direction.
~l [F]-p~(s,t )=(1-t )u(s, O) +
tu(s, 1), 0~
~2[F] =-p2(s, t) = (1 - s)u(O, t) + su(1,t), 1139
0 ~< t ~< 1.
(1)
1140
YONG LIU Ibl
V.,
u[0,1]~
uP.ll
p.q
utO.q~ utO,OI
ale.hi
utl,Ol
$
)x [UI
P.01 Fig. 1. (a) The mesh after mapping, (b) the actual mesh region.
The bilinear projector is defined as shown in Fig. l(b).
{ (~,~,)[V] _={~, [Vl + ~'~[rl - ~, [~} = e.(s, t) Pn(s, t) = (1 - t)u(s, O) + tu(s, 1) + (1 - s)u(O, t) + su(1, t) - (1 - s)(1 - t) x u(O, O) - (1 - s)tu(O, 1) - s(1 - t)u(1, O) - stu(1, 1)
(2)
where u(0, t), u(1, t), u(s, 0), u(s, 1) are the boundary functions of the region F. u(0, 0), u(0, 1), u(1, 0), u(1, 1) are the corner points at the region. This projector represents a curvilinear coordinate system created by a mapping of the unit square on to F. It may be called the transfinite bilinear Lagrange interplant of F; the transfinite mappings may be used to map any mesh in any mesh in the appropriate primitive polygon to the actual region. Thus, any mesh topology may be used. However, it is generally convenient to sacrifice topological generality by choosing the intersections of the constant coordinate curves for use as node points in order to minimize user input. The bilinear projectors thus create a natural partition of the region composed of quadrilateral elements. This is very useful for the automatic meshing generation of 2D-viscoplastic BEM. But the above eq. (2) has two drawbacks. One is that it cannot give the suitable mesh space for different local stresses. The other is that it cannot contol the size of the nonlinear region of a structure. In order to alleviate the two drawbacks, we present a modified transfinite method. The method is given as follows: Let us define a new projector goM and introduce the mesh spacing controlling parameter p~, i = 1, 2, 3, 4, with the region size controlling parameter ei, i = 1, 2, 3, 4. Thus, we have ~M(s, t) = (1 -- t)u(s p~, O)el + tu(s p:, 1)e2 + (1 -- s)u(O, tP3)e3 + su(l, tP~)e, -- (1 -- s) x (1 - t)u(O, O) - (1 - s)tu(O, 1) - s(1 - t)u(1, 0) - stu(1, 1).
(3)
With Pi and ei, we can control the mesh spacing and mesh region. Equation (3) is the modified transfinite method. Three basic equations for 2D-viscoplastic BEM.
(a) The boundary integral equations The starting integral equation for the boundary element method applied to initial strain problem is as follows
Automatic mesh generation for computing symmetric local stresses
1141
where F represents the boundary of the region and f~ is interior. The tensors u*, p~ and a ~ corresponding to the fundamental solutions have been presented in previous publications [7, 8] and the expression for g~j is the following form: gi/=
G 4(1 - v ) [ 2 ~ + ( 1 -4v)iT,~6~j],
for 2D plane strain
(6)
where for plane stress v is replaced by ~, E~ is the initial strain. The application of eq. (4) in discretized form to all boundary nodes generates the following matrix relationship: Hfi = Gp + Qia
(7)
in which matrix Q stands for the initial stress integral. Similarly, the stress rates at the internal points are given by eq. (5) in discretized form # = G'0 - H'fi + Q,~a
(8)
where matrix Q' represents the initial stress integrals plus the independent terms f0. Matrices H, G, H' and G' are the same as those obtained for pure elastic analysis [7]. Taking into consideration that for a well-posed problem a sufficient number of tractions and boundary displacement is prescribed, eqs (7) and (8) can be further written as A~ = f + Qi ~
(9)
6 = - A ' ~ + [" + Q'/~
(10)
and
where vector ~ is formed by the unknown tractions and boundary displacements and the contribution of the prescribed values is included in vectors [r and i". From the multiplication of eq. (9) by A-' comes
{
x
= R ~ a --I- lil
R =
A-IQ,
ill =
A-if.
(1 1)
Substituting eq. (11) into eq. (10) yields d = Vi ~ + fi where
V = Q' -
A'R,
(12)
fi = i" - A'fit.
Note that the elastic solution to the rate problem is given by the vectors ~h and ~. From the above it is seen that eq. (12) represents a system of ordinary differential equations for stresses at selected boundary nodes and internal points which can be solved by standard methods, producing a unique solution for rate dependent problem. (b ) Viscoplasticity relations
It is known that one explicit model which has wide applicability for isotropic material is offered by the following viscoplastic flow rule [6, 10]. i s = i p = ~, < 4 , ( f )
>
~-~. dQ
(13)
We restrict ourselves to associated plasticity situations, in which case F = Q and a sufficiently general expression for viscoplastic strain rate is given by a power law for ~b(F) such as ~b(F)
= ~ F _ FoT¢ ,
L F0 ]
r-to or ~b(F)=eS2(-~-*) - 1
(14)
where N,, N2 are prescribed numbers, F is the yield function, the onset of viscoplastic behaviour is governed by a scalar yield condition of the form F ( a ) - Fo = 0
(15)
1142
YONG LIU
in which F0 can be the uniaxial yield stress, F0 may itself be a function of an effective strain gp and for a linear strain hardening response, we have F0 = try + n ' g " = try + f(ffP)
(16)
where try is the initial uniaxial yield stress, H ' is the viscoplastic strain hardening parameter, ¢P is the current effective viscoplastic strain and f(¢P) can be determined by experiment. (c ) The computational algorithms The initial distribution of internal stress and boundary traction and displacements is obtained from the solution of the corresponding elastic problem using eqs (11) and (12) where i s = i p = O. The load increments are applied as discrete steps and thus Af. = 0, Af~ = 0 for all time steps other than the first within an increment. (i) Starting from known values of X,, a,, d at a time and compute the rate of viscoplastic strain i p by the constitutive relations (13-16). (ii) Calculate the ~P as a~{ = i ~ a t .
(17)
¢,P+t = ¢~ + A ~ .
(18)
(iii) Compute the value of the boundary unknowns AX. = R A ~ + Am
(19)
X,, + 1 = X. + AX,.
(20)
(iv) Determine the internal stresses Aa. = VA~ + An
(21)
O'n + l = O'n "q- Ao'n.
(22)
We arrive at all the starting values of the next step for which the same process can be started in the next interval. In order to improve the accuracy of this explicit time step algorithm which is equivalent to the Euler extrapolation procedure, an inner iterative loop can be imposed on the steps (i-iv) where step (ii) is replaced by calculating A~P.= (iP. + iP.+~)at./2.
(23)
i p. + I = iP
(24)
Taking initially and repeating the process to (iv) until some convergence criterion is reached using the second viscoplastic strain invariant for instance. The iP,+t is defined by using a limited Taylor series expansion such as i t + , = i p. + H.Aa. where
(25)
a[i~] T H.
= - -
a[~.]
=
[H,,(G,)].
The matrices [l-I. ] depend on the stress level o. and can be easily evaluated with the constitutive laws (13). The viscoplastic strain increment can be written as A ~ = At,[(1 - 0)i~ + 0i~+ ~]
(26)
where for 0 = 0 represents the Euler time integration scheme, the case 0 = 1/2 results in the so-called implicit trapezoidal scheme i.e. eq. (23) and 0 = 1 gives the fully implicit scheme, we get A ~ = ~ A t , + 0[I-I,]At, Aa,.
(27)
Automatic mesh generation for computing symmetric local stresses
1143
Using the incremental from (21) and substituting for A~ from (27), then (27) becomes Aa. = Vi.PAt. + 0VH.At.A~r. + An.
(28)
Equation (28) can be solved iteratively for the increment Aa. and the iterations go on until some convergence is achieved using the ratio ] ,]Aa. ,1~ - ,]Act. 1,~- ' ]
ep
(29)
where Ep is the convergent tolerance. Following the method used by Zienkiewicz and Cormeau the magnitude of this time step is controlled by a factor z which limits the maximum effective viscoplastic strain increment Ag~ as a fraction of the total effective strain ~,, so that Ag.p = ~ 2 .p -p At. ~<"re,,
(30)
where AgP. must be computed to satisfy (30) at each internal point and the least value taken for computation. The above time step limiting method is basically empirical and z in the range of 0.01-0.15 is found to be effective for explicit algorithms and for implicit scheme values of T up to 1 have been found to be stable with BEM. A further limit is generally imposed where the change in the time step length between any two intervals is limited according to Atn+ l ~< 1.5At..
(31)
3. NUMERICAL EXAMPLES In order to show the advantages of automatic meshing generation for 2D-viscoplastic BEM by using the transfinite mapping method and to demonstrate its effectiveness, the following examples are used (a) a plate with hole under uniform tension (b) a plate with V-shaped notch under uniform tension (c) a plate with U-shaped notch under uniform tension.
Example (a). A plate with hole under uniform tension The problem considered is a square 10 x 10 mm 2 plate with a circular hole under uniform tension, as shown in Fig. 2 [9]. The radius of the hole was taken as 5.0 mm, owing to symmetry only a quarter of the plate need to be taken for analysis, as shown in Fig. 3. Considering the plate with the hole under uniform tension, the viscoplastic deformation takes place firstly near the hole of the plate, so it needs taking the AFGE as for generating the 2D grids of the viscoplastic region near the hole, Four functions in eq. (3) are taken by the expressions
J
R(1 + e)cos(2 s )
COS
U(Sp2, 1)=
R(1 + e)sin(2 s )
(32)
u(1, tp4)= IR oeRtP 1 where Pt -'-P2 = 1, P3 = P4 = P. e is the parameter of controlling the viscoplastic region, p is grid space regulating parameter. When p > 1, the generated meshing graph is denser near the AE curve. When 0 < p < 1, the generating meshing graph is sparse near the AE curve, as shown in Fig. 4. In order to calculate the actual stress field, the parameter e is taken as e = 0.9 and the meshing space regulating term p = 1.2. The final generating meshing graph is shown as Fig. 5. The problem to consider is a 20 x 36 mm 2 plate with a V-shaped notch under uniform tension, as shown in Fig. 6(a), owing to symmetry only a quarter of the plate needs to be taken for analysis as shown in Fig. 6(b).
1144
Y O N G LIU
I
0 lOmm 4 O"a
!
Q
i - - 4
36 m m Fig. 2. The plate with the hole under uniform tension.
D
C
E
F
I Fig. 3. Quarter of the plate.
Let us consider the plate with a V-shaped notch under uniform tension. The initial viscoplastic region is taken as ABFG for computation. Four functions in eq. (3) are taken by the expression
u(s p2,1)
=~ es'(L--t,)+tle ] LsP(a+tle)+tl--tle (33)
Let"1
I[l[(O,tP3)~-Itlet__.:ltPtll,
u(l, tp,) = F k a + tl_l
(a)
(b) I
i
I
I
I
I
I
I
I
f
I
I
I
I
Fig. 4. The effect of grid space regulating term p. (a) p = 1.3. (b) p = 0.4.
q
I
I
f
I
I
I
I
I
Automatic mesh generation for computing symmetric local stresses
1145
(a)
(b)
Fig. 5. The meshes generated by bilinear projector. (a) Mesh 1. (b) Mesh 2.
where Pm= P2 = P3 = P 4 = P, e is used to define the viscoplastic region, p is the spacing regulating parameter. Its effect is the same as in example (a). The parameter e = 2/3 and p = 1.2, and the generating meshing graph is shown in Fig. 7. The material and load parameters are taken as [10] E = 70,000 MPa, ay = 243 MPa, v = 0.2, h' = 0.0, the tension load tra is taken in 2a,/ay = 1.3 for plane stress. The viscoplastic constants [6]: NI = 1, y = 0.001 day, 0 = 0.0, and 0 = 1.0, • = 0.01. By using the above 2D-viscoplastic BEM program together with the mesh in Fig. 7(b), the calculation results are showed in Figs 8 and 9. Figure 8 shows that as the days increase, the viscoplastic displacement with two computing schemes (0 = 0 and 0 = 1.0) at point C approaches the steady state solution, i.e. the elastoplastic solution. Figure 9 shows the comparison between the results of viscoplastic solution and elastoplastic solution by BEM. As the steady state solution o f viscoplastic BEM is reached, good agreement has been obtained between the two methods. This shows the effectiveness of the method in this paper.
Example (c). A plate with U-shaped notch under tension The problem considered is a plate with U-shaped notch under uniform tension, as shown in Fig. 10(a). Owing to symmetry only a quarter of the domain needs to be taken for analysis as shown in Fig. 10(b). The specified boundary conditions are traction normal to the boundary equal to along CD. Zero displacement normal to the boundary along AB and BC and zero traction otherwise.
1146
Y O N G LIU
(a)
90"
"llt-'--
o'a
(lra
"lit"---
I_
36ram
_
vI
(b)
B G
D
L ).
Fig. 6. (a) A plate with V-shaped notch. (b) Quarter of the plate.
Considering the plate with U-shaped notch under tension, the nonlinear region may be firstly taken place near the U-shaped notch. Hence the region A F G E should be taken as the computational nonlinear region. The four functions of eq. (3) are taken as in eq. (32) according to the actual specimen of Fig. 10(a). The material constants are taken as follows: R = 0.0425 m, plate thickness b = 0.0041 m A F = E G = e . R = 0.06375 m, E = 70 GPa, v = 0.27, ay = 320 MPa, N = 1, v = 0.001/day, 0 = 0.0, = 0.01. By using meshes of Fig. ll(a) and (b) and viscoplastic BEM program[15], the stress concentration factor K at root A in Fig. 10 can be calculated in Table 1 together with results in refs [12-14] Table 1. Comparison of stress concentration factor K Result of
Result of
Result o f
Result of
Result of
Fig. 10(a)
Fig. 10(b)
rcf. [12]
ref. [131
rcf. [141
2.14
2.07
1.94
2.11
2.0
. ~
I
I
I
I
I
(
1 I
I
!
I
I'
I
I
I
\ \ \ \\,
I
I
I
I
Fig. 7. Two meshes generated by bilinear projector of modified transfinite mapping method. (a) Mesh 1. (b) Mesh 2. 38 37.9 37,8 E E
37.7 37.6 37,5
1~ 37,4 ! 37,3 37.2 37,1 37
0
I
I
I
I
2
4
6
8
I
I
I
I
I
I
I
I
I
I
10 12 14 16 18 20 22 24 26 28 30 Days
-e-VPBEM--1 (O = 0.0)
~
VPBEM--2 ~ (O = 1.0)
Elastoplastic
Fig. 8. The displacement at point C vs times. Loading-disp. curves 1.2
-
f./.o
I.C
~" O
/f
0.8j
"~
0.~
~
0.4
"~
0.2
~ 1 ' 5
I 10
I 15
I 20
I 25
30
uc/lxE-4 ~'" Elutoplastic B E M - - - o - - - V i ~ o p l u t i ¢ BEbl Fig. 9. Comparison of the viscoplastic steady state solution with elastoplastic solution.
1148
YONG
LIU
(a) q •
q
D
0
f
q
L I-
im
.I
0.890 m
(b)
~XX~XX ~ X ~ XX
E Fig. 10. (a) and (b).
(a) I
I
I
I
I
I
I
I
I
I
I
I
I
I
(b) l
I
l
l
I
l
I
I
I
Fig. 1 I. (a) and (b).
l
l
I
I
l
Ilr I
~
Automatic mesh generation for computing symmetric local stresses
1149
5. CONCLUSIONS The modified transfinite mapping method is proposed in 2D-nonlinear BEM for the automatic mesh generation of specimens with symmetric notches. Although the computational efficiency of general BEM programs has steadily increased over the past years, the time consuming process of accurate data preparation is still a burden for most BEM users in engineering. The modified transfinite mappings based on boundary curve descriptions provided an effective basis for automatic generation of visco-plastic BEM. They are generally, simple to implement and computationally efficient. Grids produced by this method have predictable forms wherein the interior of the mesh accurately reflects boundary curvature information. The meshing spacing regulating parameter proposed in this paper can control the grid space and mesh size effectively. The numerical examples show that the method in this paper has the advantages of high computational efficiency, precise modelling of boundary minimum input data and is good for computing the local stress of 2D specimens. Therefore, the modified transfinite mapping method can have applications in the problem of automatic mesh generation for nonlinear BEM more and more. It is useful for the 2D-viscoplastic BEM analysis for preliminary research and Alexander yon Humboldt foundation. Acknowledgements--The project was supported by Zhejiang Science Foundation.
REFERENCES [I] S. W. Chae and K. J. Bathe, On automatic mesh construction and mesh refinement in finite element analysis. Comput. Structures 32, (3/4) 911-936 (1989). [2] J. Zhu, A hybrid differential-algebraic method for three-dimensional grid generation. Int. J. numer. Meth. Engng 29, 1271-1279 (1990). [3] W. J. Gordon and C. A. Hall, Construction of curvilinear coordinate system and applications to mesh generation. Int. J. numer. Meth. Engng 7, 461-477 (1973). [4] C. A. Hall, Transfinite interpolation and applications to engineering problems, in Theory of Approximation (Edited by Law and Sahney) pp. 308-331. Academic Press (1976). [5] W. A. Cook, Body oriented (natural) co-ordinates for generating three dimensional meshes. Int. J. numer. Meth. Engng 8, 27-43 (1974). [6] O. C. Zienkiewicz and I. C. Cormeau, Visco-plasticity, plasticity and creep in solids, A unified numerical solution approach, Int. J. Numer. Meth. Engng 8, 821-845 (1974). [7] C. A. Brebbia, J. C. F. Telles and L. C. Wrobel, Boundary Element Techniques: Theory and Applications in Engineering, Springer, Berlin (1983). [8] Liu Yong and Shen Najie, Residual stress analysis for actual material model of autofrettage tube by nonlinear BEM. Int. J. Pres. Ves. Piping 48, 21-35 (1991). [9] P. S. Theocaris and E. Marketos, Elastic-plastic analysis of perforated thin strips of a strain hardening material. J. Mech. Phys. Solids 12, 377-396 (1964). [10] G. C. Nayak and O. C. Zienkiewicz, Elastoplastic stress a generalization for various constitutive relations including generalization for various constitutive relations including strain softening. Int. J. numer. Meth. Engng 5, 1113-1135 (1972). [I 1] R. Haber, M. S. Shephard, J. F. Abol, R. H. GaUagher and D. P. Groenberg, A general two-dimensional, graphical finite element preprocessor utilizing discrete transfinite mappings. Int. J. numer. Meth. Engng 17, 1015-1044 (1981). [12] Yong Liu, Naijie Shen, and Wai Jiang, An effective method for calculating the local elastoplastic stress of a tension specimen with symmetric U shaped grooves. Engng Fracture Mech. 43, 1085-1091 (1992). [13] K. Axelsson and A. Samuelsson, FEM analysis of elastoplastic materials displaying mixed hardening. Int. J. numer. Meth. Engng 14, 211-225 (1979). [14] J. R. Crews, Local plastic stresses in aluminium alloy specimens with stress concentration factor and constantamplitude loading. NASA TN 3813, Washington (1965). [15] Yong Liu and N. J. Shen, An BEM implementation for elastoplastic analysis by using the steady state solutions of viscoplastic problem. The computational mechanics and its applications in engineering. Proe. 3rd-Computational Mechanics Conf. East China. China University of Science & Technology Ltd. pp. 41-46 (1992). (Received 17 February 1993)
EFM 52/6-~L