- Email: [email protected]
Elasto-plastic analysis using a coupled boundary element finite element technique
EngineeringAnalysiswithBoundaryElements14 (1994) 39-49 © 1994 Elsevier Science Limited Printed in Great Britain. All fights reserved 0955-7997/94/$07.00
ELSEVIER
Elasto-plastic analysis using a coupled boundary element finite element technique J. L. Wearing & M. C. Burstow Department of Mechanical and Process Engineering, University of Sheffield, Mappin Street, Sheffield, UK, S1 3JD
A combined finite element boundary element method for the analysis of elastoplastic stress analysis is discussed in this paper. The merged technique utilizes the Regular indirect discrete method which has been modified to enable coupling of the two methods to be achieved and maximizes the advantages of each method while minimizing their disadvantages. The merged technique is shown in particular to be useful for the analysis of elasto-plastic problems by modelling the plastic region, which normally only occupies a small area of the problem being analysed, by the finite element method and the large elastic region using the boundary element method. This approach avoids the problems associated with body diseretization in the boundary element method and the requirement for finite element models with large numbers of degrees of freedom when either of these methods are used on their own for elasto-plastic analysis. Results for the analysis of four two dimensional elasto-plastic problems are presented. These results compare well with the results from other sources and indicate that a successful union of the finite element method and the boundary element method can be achieved.
Key words: Elasto-plasticity, boundary elements, finite elements, regular indirect BEM, fracture mechanics.
INTRODUCTION
to faster analysis times and makes the BEM particularly attractive for the linear stress analysis of three dimensional components. Where plastic regions exist the BEM becomes less attractive as, for such an analysis, volume discretization is required, with a consequential increase in the number of degrees of freedom. The FEM and the BEM have advantages and disadvantages compared with each other but they do not have shared strengths and weaknesses and a logical step is to develop a combined B E M - F E M technique which exploits their strengths and minimizes their weaknesses, and to use the combined technique in situations where it is appropriate. The combined method is particularly advantageous for the analysis of elasto-plastic problems since the FEM only needs to be applied to the small plastic region and the BEM to the remaining large elastic region. The result is a considerable reduction in the number of degrees of freedom, which reduces the size of the resulting stiffness matrix and yields a significant improvement in the solution time compared to the FEM. Further economies are achieved as the BEM region remains elastic and it is therefore only necessary to recalculate the stiffness matrix in the finite element region at each load increment. A suitable technique for dealing with plastic regions is
The extensive industrial use of the finite element method (FEM) for engineering stress analysis and the increased size and speed of modem digital computers has been coupled with a constant demand for the stress analysis of three dimensional components of increasing complexity. These complex three dimensional components do however require a large number of elements and degrees of freedom to model them successfully, leading to extremely time consuming procedures for the preparation of the initial models and the subsequent analyses, and although the finite element method is particularly suited to the analysis of problems with plastic zones, the calculational time can be dramatically increased for such problems, despite the fact that the plastic zone only occupies a small part of the component. The requirement for increased efficiency for the analysis of complex three dimensional components has led to the emergence of the boundary element method (BEM) as an alternative method for engineering stress analysis. As the boundary element method only requires the surface of the component to be analysed, fewer degrees of freedom are required for the initial model and this leads 39
40
J. L. Wearing, M.C. Burstow
particularly important in fracture analysis as the presence of cracks in components significantly reduced their strength and is therefore an important consideration in assessing their structural integrity. Fracture analysis is governed by linear elastic fracture mechanics (LEFM) and elasto-plastic fracture mechanics (EPFM). The former is ideally suited to analysis by the BEM and there are numerous examples of the application of the BEM for LEFM. 1'2 Although the BEM is superior to the FEM for LEFM, the BEM is less attractive for EPFM, which accounts for the plastic zone, and therefore requires body discretization in the plastic region. Consequently the combined FEM-BEM method is ideally suited for EPFM. Various techniques have been used to achieve coupling of the FEM with the BEM. Kelly et al.3 have described various methods of linking the FEM with the direct and indirect versions of the BEM. The coupling technique used in the work discussed in this paper is based on an indirect approach, as discussed by Scholfield, 4 known as the regular indirect discrete method (RIDM). In the normal version of the indirect BEM, as discussed by Bannerjee and Butterfield,5 fictitious sources are assumed to be distributed continuously over the boundary of the problem being analysed. However this approach leads to singular integrals when the source and field points coincide. Additionally, poor results can be obtained at structural discontinuities. The problems with singular integrals are avoided by using the regular indirect boundary element method 6 in which the fictitious sources are assumed to be distributed on a boundary, which is placed outside the actual boundary of the problem being analysed. Singular integrals are eliminated when the regular method is used and improved results are obtained at discontinuities. The regular indirect discrete method brings additional benefits by replacing the continuously distributed sources with sources distributed at discrete points on the source boundary and consequently avoids the complex integrals normally associated with BEMs. The work discussed by the authors uses the RIDM to create a boundary element region which can be coupled to a finite element region. This combined technique has been used successfully by Wearing eta/. 7'8 for the analysis of two and three dimensional linear elastic problems. One of the principle advantages of the combined method compared to other numerical techniques, is in the analysis of problems with plastic regions and in this paper the combined method has been applied to the solution of two dimensional elasto-plastic problems and for the analysis of two dimensional EPFM problems. Case studies, illustrating the use of the combined approach, are presented and the results, which are obtained from the combined method, compare favourably with those obtained using alternative analytical and numerical techniques.
T H E COMBINED FINITE ELEMENT BOUNDARY
ELEMENT METHOD The two main numerical techniques which are currently used for stress analysis are the boundary element method and the finite element method. Each method, however, has advantages and disadvantages compared to the other. The BEM on the one hand requires simpler initial models than the FEM, it is ideal for linear elastic stress analysis, it is more appropriate for problems extending to infinity but is less suitable for problems with complex geometries and for the analysis of components with plastic zones. On the other hand, the FEM requires more complex initial models than the BEM but it is more convenient for problems with complex geometries and for analyses which require data at a large number of internal points. Obviously the two methods do not have shared strengths and weaknesses and an ideal approach, where appropriate, is a combined finite element - boundary element method to enable the strengths of each method to be maximized and their weaknesses minimized. Due to their distinct mathematical formulations, the FEM and the BEM cannot be directly linked. The basic variables in the FEM are nodal displacements and forces and in the BEM the variables are surface tractions and displacements. It is obvious therefore that the basic formulation of one of the methods must be altered to make it compatible with the other to allow a satisfactory union of the two methods to be achieved. In the combined method, adopted for the work discussed in this paper, the boundary element equations have been modified to make them compatible with the finite element equations. This approach has the potential for wide industrial appeal as the FEM is more established than the BEM in the industrial environment and enables the combined approach to be used with minimal modifications to existing finite element software. The finite element method
The FEM for stress analysis is based on the minimization of the potential energy function which is given by the expression (e.g. Hinton and Owen9):
pE= Iv W(U)dV- IsUbdV- IsupdS
(1)
in which W(u) is the strain energy density, u is the deflections, b is the distributed body forces, p is the distributed surface forces and S and V denote surface and volume integrals respectively. On minimizing eqn (1), the following expression relating nodal forces {F} and nodal displacement {u} by the structural stiffness matrix [K] is obtained. {F} = [K]{u}
(2)
where the stiffness matrix is given in terms of the
Elasto-plastic analysis using a coupled boundary element finite element technique
the actual boundary of the problem are obtained from the RIDM at collocation points relative to the source points. Hence
/
j=p ui(xf) = Z Gu(xfxs)cpj(xs) j=i
(a)
j=p ti(xf) = Z Fij(xrxs)qSj(Xs) j=l
(b)
....
Fig. 1. Source distribution in indirect boundary element methods.
strain-displacement matrix [B] and the stress-strain matrix [D] as (3)
The regular indirect discrete method
The regular indirect discrete method (RIDM) is a development of the indirect boundary element method (IBEM). The normal formulation of the IBEM for linear elastic stress analysis is based on the use of the fundamental (Kelvin) solution 5 which satisfies the governing differential equation. For a homogeneous, isotropic domain, the relevant fundamental solutions are related to the displacements ui(xf) and the surface tractions ti(xf) at a point xf on the boundary of the components being analysed by a fictitious source q~distributed on the boundary as shown in Fig. l(a). Hence
(8)
{t} = [F]{~b}.
(9)
From eqns (8) and (9) the boundary tractions and displacements can be related by the following expression: {t} = [F][O] -l{u}.
(10)
The combined method
It is obvious, on examination of eqns (2) and (10), that the FEM cannot be linked directly with the RIBM. The latter method can however be modified to make it compatible with the FEM by using Betti's second law which is:
JvuA(u) = I v W ( u ) d V - JsUtdS
(11)
in which A is a differential operator. If eqn (11) is substituted into eqn (1) the following expression is obtained
aE=:Jvua(u)dV+:IsutS-IvubdV-JsupdS1
1
(12)
(4) In the absence of body forces eqn (12) becomes:
and
ti (Xf) : IS Fij (XfXs)~ dS.
(7)
and
(c)
Ui(Xf) = IS Gij(XfXs)O dS
(6)
in which p is the number of source points. Matrix equations for the deflections and tractions at all points on the boundary can be obtained from eqns (6) and (7). Hence
Node
{u} = [a]{~}
[rl = aIv [Blw[DI[BIdV
41
PE=~
(5)
When the field point, xf and the source point, Xs, in eqns (4) and (5) coincide the integrals become singular. However Wearing et aL6 showed that the problems associated with the singular integrals can be overcome by distributing the source on a boundary which is placed outside the actual boundary of the problem as shown in Fig. l(b). A further development 4 has resulted in the RIDM in which the sources q~j(xs) are distributed at discrete points around the source boundary as shown in Fig. l(c). By removing the integrals in eqns (4) and (5), the displacements and tractions at each point on
'Is u t d S - IsupdS.
(13)
Substituting eqn (10) into eqn (13) and minimizing the resulting expression gives the following equation:
Is[F][GI-I {u} dS = Is{q} dS
(14)
which can be written in matrix form as:
K{u}
= {e}
(15)
Equation (15) now has the same form as eqn (3) which means that the RIDM can be linked directly with the FEM. When using the combined F E M - R I D M it is convenient
J. L. Wearing,M.C. Burstow
42 to express the finite element equations as:
{F}f¢----[K]fe{u}f e
(16)
and the modified boundary element equations as: {F}be = [r]~{U}be
(17)
in which the subscripts fe and be refer to the finite element and boundary element regions respectively. When eqns (16) and (17) are combined a system of equations is obtained for the complete component. Hence {F} = [r]{u}
(18)
in which [K] is the merged stiffness matrix obtained from eqns (16) and (17) and may be written as: [K] ---- [Klfe + [K]be
(19)
and {F} and {u} are the merged force and displacement vectors.
ELASTO-PLASTIC ANALYSIS In an elasto-plastic analysis eqn (18) is solved initially to give a global elastic solution to the problem. An incremental process is then employed whereby the applied loads are incremented according to specified load factors. As the plastic behaviour is confined to the finite element region, only [K]fe requires to be recalculated in eqn (19) after each load increment and, as the boundary element region remains elastic, [K]be remains unaltered and does not therefore require to be calculated after each load increment. [K]fe can be rewritten as: [K]f e =
Iv[B]TtD]ep[B]d g
(20)
where [D]~p is the elasto-plastic stress-strain matrix and is given by the expression (e.g. Owen and Hintonm). [D]ep = [ D ] -
{do}{dD}V H ' + {dD}r{a}
(21)
in which {dn} = [Dl{a}
(22)
and {a} is the plastic flow vector given by: {a}T
0(_~] FDF OF OF OF]
(23)
In eqn (23) F(tr, K) governs the yield criterion, with K being a parameter which depends on the specific yield criterion which is used in the elasto-plastic analysis.
CASE STUDIES The following four case studies are presented to test the
(a)
(b)
Fig. 2. Finite element and boundary element for two dimensional analyses. effectiveness of the F E M - R I D M technique for the analysis of two dimensional problems with plastic zones. (a) A rectangular plate with a central circular hole. (b) A rectangular plate with symmetrical vee notches on opposite edges. (c) A centre cracked rectangular plate. (d) A four point bending specimen. For each case study, with the initiation of plastic deformation, the load was applied in increments and, since plastic deformation is confined to the finite element region, only [K]fe in eqn (19) needs to be recalculated after each load increment. As the boundary element region remains elastic throughout the analysis [K]be retains its initial values. For all case studies the finite element regions were modelled using eight noded isoparametric elements as shown in Fig. 2(a) and the boundary element regions were modelled using three noded isoparametric line elements as shown in Fig. 2(b). As the elements in each region are compatible, merging of the two regions can be achieved. The source boundary is placed outside the actual boundary of the component being analysed when the RIDM is used and previous studies, 6 on the application of the regular boundary element method for stress analysis, have indicated that there is a range of possible positions of the source boundary which yield good, stable results. These studies have revealed that when the source and component boundaries are too close or too far apart unsatisfactory results are obtained. The position of the source boundary is defined by a nondimensional factor A, which is the ratio of the distance of the source boundary from the boundary of the component to the length of the boundary element. This means that for boundary element models with elements of different lengths, the source boundary does not require to be at a fixed distance from the boundary of the component being analysed. For the case studies considered in this paper good stable results were obtained for values of ), between 1.0 and 10"0. In all analyses the value chosen for ), was 4.0. Due to the symmetrical nature of the problems discussed in the paper it was only necessary to analyse a quarter of the plate in each of the first three case studies and in the four point bending specimen the analysis was confined to half of the beam. All analyses were plane strain with a Von Mises yielding criterion and the first
Elasto-plastic analysis using a coupled boundary element finite element technique
43
--,--D.-
(b)
(a)
t
0.4~ 00
o.'s:4~
O.561
0/56 (d)
(c)
Fig. 3. Plate with central circular hole. Poisson's ratio = Uniaxial yield stress = Work hardening parameter =
three case studies had the following material properties: Young's modulus Poisson's ratio Uniaxial yield stress Work hardening parameter
= = = ---
7000 N/ram 2 0.33 24-3 N / m m 2 0-0 (perfectly plastic)
The material properties for the four point bending specimen were: Young's modulus = 200 000 N / m m 2
With respect to the fracture analyses an L E F M and E P F M analysis was undertaken for each case study. Quarter point elements were used in the finite element mesh to model the singularity at the crack tip. The results for the fracture analysis, using the combined method, were used in conjunction with a number
v,,k
.--.-I)-
,
(b)
(a)
O. 7 2 5 \
o.6
!: ....
0"28 400 N/ram 2 0"0 (perfectly plastic)
-"~-~
!. (d)
(c)
Fig. 4. Notched plate.
44
J. L. Wearing, M.C. Burstow 10mm
--
r
:
"
'
-/
/:
))
', Z . . . .
\-- ...........
-/ . . . .
I I I I t
I I 1 I I
15ram
/
/
Crack
J - Path
J - Path 8
h u'-
) . I Ill
Fig. 7. J integral paths for centre crack plate. 10ram v I
II I l l
Fig. 5. Centre cracked plate. of approaches to determine the stress intensity factors, the results of which are given in Tables 1 and 2 (later). Rectangular plate with central hole The plate has a length of 36 mm, a width of 20 mm and
_
)!
Y
l l l l l l l
the diameter of the central hole is 10ram. Details of the loading and the boundary conditions for the quarter plate model, used in the analysis, are given in Fig. 3(a). An initial tensile load which ensured an elastic stress distribution throughout the plate was applied initially as shown in Fig. 3(a). This load was scaled to a value which caused the initiation of yield at the edge of the hole and was followed by a number of load increments up to a value of 13.78 N/mm 2, which is 0"56 times the value of the yield stress of the material. The results showing the growth of the plastic zone, adjacent to the hole, are shown in Fig. 3(d) and indicate the extent of the plastic zone for each load increment, which is indicated as a percentage of the yield stress used as the applied load. The region adjacent to the hole containing the plastic region was modelled using 20 finite elements, having 130 nodes and the remainder of the plate was modelled using 16 boundary elements with a total of 32 nodes. Results were also obtained from a finite element model with 42 elements and 160 nodes. Details of the meshes are shown in Figs 3(b) and
3(c). Notched plate Details of the loading and boundary conditions for the quarter plate model are shown in Fig. 4(a). The length of the plate is 36 mm, its width is 18 mm and the depth of the 90 ° notch is 5 mm. The problem was analysed using two boundary element regions with a total of 34 Table 1. Stress intensity factors for centre cracked plate Method of calculating Ki (1)
.......
i i i i i
Fig. 6. Combined RIDM-FEM model for centre cracked plate.
(2)
LEFM analysis (a) Combined method Displacement extrapolation Strain energy release rate J-integrals (b) Rooke and Cartwrightll (c) J integrals using FEM EPFM analysis (a) J integrals using combined method (b) Theoretical plastic z o n e c o r r e c t i o n 12
Ki
2.31 2.33 2.32 2.31 2.28 34.85 34"26
Elasto-plastic analysis using a coupled boundary element finite element technique 140
45
loads being indicated as a percentage of the yield stress used as the applied load.
-
130 •
Centre cracked plate II0
9O 8O
~
50
Sl~-ess(N/ram2) ~0 lied
30
4.471
6.901
20
io 0
3.136 "~T 0.00
0.05
I
I
I
I
I
I
I
0.10
O. I S
0.20
0.~
0.30
0.35
0.40
0.486
Distance Across Ligament From Crack Tip (nma)
Fig. 8. Stress distribution across ligament of centre cracked plate. elements and a finite element region with 28 elements giving a total of 139 nodes as shown in Fig. 4(b). It was also analysed using 52 finite elements with a total of 185 nodes as shown in Fig. 4(c). The plate was initially loaded along its top edge with an elastic tensile load as shown in Fig. l(a). The load was then scaled to induce yielding at the notch and increased in six increments until the applied load was 0"725 times the yield stress of the material. Figure 4(d) shows the increase in the yield surface for each load increment with the
The dimensions and loading of the plate are shown in Fig. 5 with the combined R I D M - F E M mesh, having a total of 299 nodes, which was used to analyse the quarter plate model, being indicated in Fig. 6. Details of the finite element mesh and the J integral paths used in the LEFM and EPFM analysis are shown in Fig. 7. The problem was also analysed using a finite element mesh having 1000 nodes. For the LEFM analysis a tensile load was applied to the plate and this load was sufficiently small to ensure that the plastic zone, which must exist in the vicinity of the crack tip, was not large enough to encroach on any of the Gauss points close to the crack tip, ensuring a purely elastic analysis. The KI stress intensity factors obtained by a number of methods for the LEFM analysis are given in Table 1. For the EPFM analysis the KI values were obtained for a number of tensile loads ranging from 0-486 to 7.205N/mm 2 with the results for a tensile load of 7.205N/mm 2 being given in Table 1. The effect of increasing the load on the stress distribution along the plate ligament is shown in Fig. 8 and the growth of yield surface with increasing load is shown in Fig. 9.
Four point bending specimen As in the previous example, this case study presents the problem of a crack opening under mode I deformation.
Applied Stress = 4,471 N / r a m 2
Applied S t r e ~ = 5.686 N / r a m 2
Applied Stress = 6.598 N / r a m 2
Applied S t r ¢ ~ = 7.205 N/ram2
Fig. 9. Plastic zone growth for centre cracked plate. • FEM; - -
RIDM-FEM.
46
J. L. Wearing, M.C. Burstow 450
mm Crack
R = 0.254 mm ~ l / - -
19.05 mm 4
--I12.7
mm
d 31.75 mm -I
Fig. 10. Four point bending specimen. = :
= :
IIIFig. 11. Combined RIDM-FEM model for four point bending specimen. In this case however a bending load is applied and the crack is emanating from a stress concentration. The geometry of the specimen is shown in Fig. 10. The position of the supports ensures that the bending moment between the supports is constant. The crack extends 0"5 mm along the vertical centre line from the root of the notch. Figure 11 gives details of the combined R I D M - F E M mesh, which comprised four boundary element regions, with a total of 48 elements, and a finite element region for the plastic zone, with a total of 90 elements, giving a total of 361 nodes for the analysis. A finite element model comprising 180 elements with 850 nodes was also used to analyse the problem. The finite element mesh for the combined method and the J integral paths used in the LEFM and EPFM analyses are shown in Fig. 12. A deflection of 0.1 × 10-3mm was applied at the load points of the specimen for the LEFM analysis. This was small enough to ensure that the yield surface
adjacent to the crack tip was not large enough to be detected by the Gauss points. The KI stress intensity factors for the LEFM analysis, calculated from various methods are shown in Table 2. The elasto-plastic solution to the problem was obtained by scaling the elastic solution until the yield stress of the most highly stressed Gauss point was reached and the load was then applied incrementally until the deflection at the load points was 0.07ram. Figure 13 shows the growth of the plastic zone with increasing deflection and the distribution of the effect Von Mises' stress through the plate below the crack tip is shown in Fig. 14. DISCUSSION OF RESULTS The effectiveness of combining the finite element method and the regular indirect discrete method for the analysis of two dimensional elasto-plastic problems has been investigated using four case studies. Figures 3(d) and 4(d) show the results for the growth of the
/ J - Path 7
~__~J-Path 1
Table 2. Stress intensity factors for four point bending problem Method of calculating KI (1)
(2)
LEFM analysis (a) Combinedmethod Displacement extrapolation J integrals (b) Rooke and Cartwrightll (c) J integrals using FEM EPFM analysis (a) J integrals using combined method (b) Theoretical plastic zone correction12
KI
1' 15 1.17 1.16 1.17 821.5 816.7
/ ......5 ........L .......i Fig. 12. J integral paths for four point bending spcomen.
Elasto-plastic analysis using a coupled boundary element finite element technique
J/ /
/
/
AppliedDisplacement= 0.035 mm
/
47
/
AppliedDisplacement= 0.045 mm
AppliedDisplacement= 0.07 mm
Fig. 13. Plastic zone growth rate for four point bending specimen. • FEM; - -
plastic zone with increasing tensile load for the plane strain analysis of a rectangular plate with a central circular hole and a rectangular plate with symmetrical vee notches on opposite edges. The plates have been analysed using both the combined method and the FEM and the results from both analyses are in close correspondence. Due to the reduction in the number of degrees of freedom in the combined method compared to the FEM the analysis time using the R I D M - F E M approach was 30% less than the time taken for the FEM analysis. The results for the stress intensity factors from the LEFM and EPFM analyses are given in Tables 1 and 2. The results from the R I D M - F E M analyses have been obtained using a number of approaches as indicated in Tables 1 and 2. These results are compared with
RIDM-FEM.
the results from finite element analyses and from other sources and close agreement has been achieved. Figures 9 and 13 show the growth of the plastic zones, with increasing load, for the centre cracked plate and the four point bend specimen respectively and in each case good agreement with the results obtained from finite element analyses has been achieved. Considerable savings in model size were achieved when using the combined method compared to the FEM resulting in improved efficiency at the calculational phase. The results shown in Fig. 8 for the variation in stress across the ligament of the centre cracked plate show a smooth transition in the results between the boundary element region and the finite element region showing that the FEM and the RIDM can be successfully coupled.
48
J. L. Wearing, M.C. Burstow
400350
Applied displacements (ram)
300~
/
0.07
oIlll -~ 200 1so
0.035
'~ 100
o
0"3 0
I
2
3
4
5
6
7
8
Depth Through Beam Below Crack Tip (ram)
Fig. 14. Von Mises' stress distribution below crack tip for four point bending specimen.
CONCLUSIONS The results which have been obtained for the analysis of the case studies discussed in the paper clearly indicate that the modification to the R I D M to allow it to be combined with the FEM has proved successful for the analysis of two dimensional elasto-plastic problems and overcome the major disadvantages of the BEM and the FEM when the methods are used on their own for elasto-plastic analysis. The results from the combined analyses are in excellent agreement with finite element results and with results from other sources. H'12 The models for the combined approach require considerably fewer nodes than for the finite element models, which leads to a significant reduction in the size of the stiffness matrix for the R I D M - F E M analysis, compared to the stiffness matrix for the FEM. Further savings in analysis times are achieved from the combined method for elasto-plastic analysis since the boundary element region remains elastic and the stiffness matrix for this region retains its initial values. Only the stiffness matrix for the fmite element region needs to be recalculated since the plastic zone is confined to this region. Although the case studies discussed in this paper are all two-dimensional, there is clear evidence that the combined R I D M - F E M would be particularly advantageous for three-dimensional elasto-plastic analysis, since considerable savings in model size will be achieved for the combined method compared to the FEM and removes the problems associated with internal discretization when the BEM is used for such analyses.
REFERENCES 1. Aliabadi, M. H., Cartwright, D. J., Naehring, D. W. & Daney, W. A. Stress intensity factor weight functions for cracks at holes and half plane. Proceedings of the llth International Conference on Boundary Element Methods, Vol. 3, eds C. A. Brebbia & J. J. Connor. Springer-Verlag, Berlin, 1989, pp. 83-98. 2. Stok, B. & Bukovec, B. Stress and stress intensity factor analysis in cracked elastic bars under torsion by the boundary element method. Proceedings of the llth International Conference on Boundary Element Methods, Vol. 3, ed. C. A. Brebbia & J. J. Connor. Springer-Verlag, Berlin, 1989, pp. 111-22. 3. Kelly, D. W., Mustoe, G. G. W. & Zienkiewicz, O. C. Coupling boundary elements with other numerical methods. In Developments in Boundary Element Methods, Vol. 1, eds P. K. Baunerjee & R. Butterfield. Applied Science Publishers, London, 1979. 4. Scboliield, R. P. Development of the indirect discrete boundary method and its application to three dimensional design analysis. PhD thesis, The University of Sheffield, 1986. 5. Bannerjee, P. K. & Butterfield, R. Boundary Element Methods in Engineering Science. McGraw-Hill, Chichester, 1981. 6. Wearing, J. L., Sheikh, M. A., Patterson, C. & Abdul Rahman, A. G. A regular indirect boundary element method for stress analysis. Proceedings of the 9th International Conference on Boundary Element Methods, ¢ds C. A. Brebbia, W. L. Wendland & G. Khun. Computational Mechanics Publications, Southampton and Boston, 1987, pp. 183-98. 7. Wearing, J. L., Sheikh, M. A. & Hickson, A. J. On the application of a combined finite element boundary element method for design analysis. Proceedings of the International Conference on Boundary Element Technology - BETECH 89, eds C. A. Brebbia & N. G. Zamani.
Elasto-plastic analysis using a coupled boundary element finite element technique Computational Mechanics Publications, Southampton, 1989, pp. 99-110. 8. Wearing, J. L., Sheikh, M. A. & Rahmani, O. A combined finite element boundary element approach for three dimensional stress analysis., Proceedings of the International Symposium on Boundary Elements in Mechanical and Electrical Engineering, ed. C. A. Brebbia et al. SpringerVerlag, Ikdin, 1990, pp. 31-41. 9. Hinton, E. & Owen, D. R. J. Finite Element Programming. Academic Press, London, 1977.
49
10. Owen, D. R. J. & Hinton, E. Finite Elements in Plasticity: Theory and Practice. Pineridge Press, Swansea, 1980. 1I. Rooke, D. P. & Cartwright, D. J. A Compendium of Stress Intensity Factors. Her Majesty's Stationery Ofl~c¢, London, 1976. 12. Brock, D. Elementary Engineering Fracture Mechanics. Sijhoff and Noordhoff, Holland, 1978.
ELSEVIER
Elasto-plastic analysis using a coupled boundary element finite element technique J. L. Wearing & M. C. Burstow Department of Mechanical and Process Engineering, University of Sheffield, Mappin Street, Sheffield, UK, S1 3JD
A combined finite element boundary element method for the analysis of elastoplastic stress analysis is discussed in this paper. The merged technique utilizes the Regular indirect discrete method which has been modified to enable coupling of the two methods to be achieved and maximizes the advantages of each method while minimizing their disadvantages. The merged technique is shown in particular to be useful for the analysis of elasto-plastic problems by modelling the plastic region, which normally only occupies a small area of the problem being analysed, by the finite element method and the large elastic region using the boundary element method. This approach avoids the problems associated with body diseretization in the boundary element method and the requirement for finite element models with large numbers of degrees of freedom when either of these methods are used on their own for elasto-plastic analysis. Results for the analysis of four two dimensional elasto-plastic problems are presented. These results compare well with the results from other sources and indicate that a successful union of the finite element method and the boundary element method can be achieved.
Key words: Elasto-plasticity, boundary elements, finite elements, regular indirect BEM, fracture mechanics.
INTRODUCTION
to faster analysis times and makes the BEM particularly attractive for the linear stress analysis of three dimensional components. Where plastic regions exist the BEM becomes less attractive as, for such an analysis, volume discretization is required, with a consequential increase in the number of degrees of freedom. The FEM and the BEM have advantages and disadvantages compared with each other but they do not have shared strengths and weaknesses and a logical step is to develop a combined B E M - F E M technique which exploits their strengths and minimizes their weaknesses, and to use the combined technique in situations where it is appropriate. The combined method is particularly advantageous for the analysis of elasto-plastic problems since the FEM only needs to be applied to the small plastic region and the BEM to the remaining large elastic region. The result is a considerable reduction in the number of degrees of freedom, which reduces the size of the resulting stiffness matrix and yields a significant improvement in the solution time compared to the FEM. Further economies are achieved as the BEM region remains elastic and it is therefore only necessary to recalculate the stiffness matrix in the finite element region at each load increment. A suitable technique for dealing with plastic regions is
The extensive industrial use of the finite element method (FEM) for engineering stress analysis and the increased size and speed of modem digital computers has been coupled with a constant demand for the stress analysis of three dimensional components of increasing complexity. These complex three dimensional components do however require a large number of elements and degrees of freedom to model them successfully, leading to extremely time consuming procedures for the preparation of the initial models and the subsequent analyses, and although the finite element method is particularly suited to the analysis of problems with plastic zones, the calculational time can be dramatically increased for such problems, despite the fact that the plastic zone only occupies a small part of the component. The requirement for increased efficiency for the analysis of complex three dimensional components has led to the emergence of the boundary element method (BEM) as an alternative method for engineering stress analysis. As the boundary element method only requires the surface of the component to be analysed, fewer degrees of freedom are required for the initial model and this leads 39
40
J. L. Wearing, M.C. Burstow
particularly important in fracture analysis as the presence of cracks in components significantly reduced their strength and is therefore an important consideration in assessing their structural integrity. Fracture analysis is governed by linear elastic fracture mechanics (LEFM) and elasto-plastic fracture mechanics (EPFM). The former is ideally suited to analysis by the BEM and there are numerous examples of the application of the BEM for LEFM. 1'2 Although the BEM is superior to the FEM for LEFM, the BEM is less attractive for EPFM, which accounts for the plastic zone, and therefore requires body discretization in the plastic region. Consequently the combined FEM-BEM method is ideally suited for EPFM. Various techniques have been used to achieve coupling of the FEM with the BEM. Kelly et al.3 have described various methods of linking the FEM with the direct and indirect versions of the BEM. The coupling technique used in the work discussed in this paper is based on an indirect approach, as discussed by Scholfield, 4 known as the regular indirect discrete method (RIDM). In the normal version of the indirect BEM, as discussed by Bannerjee and Butterfield,5 fictitious sources are assumed to be distributed continuously over the boundary of the problem being analysed. However this approach leads to singular integrals when the source and field points coincide. Additionally, poor results can be obtained at structural discontinuities. The problems with singular integrals are avoided by using the regular indirect boundary element method 6 in which the fictitious sources are assumed to be distributed on a boundary, which is placed outside the actual boundary of the problem being analysed. Singular integrals are eliminated when the regular method is used and improved results are obtained at discontinuities. The regular indirect discrete method brings additional benefits by replacing the continuously distributed sources with sources distributed at discrete points on the source boundary and consequently avoids the complex integrals normally associated with BEMs. The work discussed by the authors uses the RIDM to create a boundary element region which can be coupled to a finite element region. This combined technique has been used successfully by Wearing eta/. 7'8 for the analysis of two and three dimensional linear elastic problems. One of the principle advantages of the combined method compared to other numerical techniques, is in the analysis of problems with plastic regions and in this paper the combined method has been applied to the solution of two dimensional elasto-plastic problems and for the analysis of two dimensional EPFM problems. Case studies, illustrating the use of the combined approach, are presented and the results, which are obtained from the combined method, compare favourably with those obtained using alternative analytical and numerical techniques.
T H E COMBINED FINITE ELEMENT BOUNDARY
ELEMENT METHOD The two main numerical techniques which are currently used for stress analysis are the boundary element method and the finite element method. Each method, however, has advantages and disadvantages compared to the other. The BEM on the one hand requires simpler initial models than the FEM, it is ideal for linear elastic stress analysis, it is more appropriate for problems extending to infinity but is less suitable for problems with complex geometries and for the analysis of components with plastic zones. On the other hand, the FEM requires more complex initial models than the BEM but it is more convenient for problems with complex geometries and for analyses which require data at a large number of internal points. Obviously the two methods do not have shared strengths and weaknesses and an ideal approach, where appropriate, is a combined finite element - boundary element method to enable the strengths of each method to be maximized and their weaknesses minimized. Due to their distinct mathematical formulations, the FEM and the BEM cannot be directly linked. The basic variables in the FEM are nodal displacements and forces and in the BEM the variables are surface tractions and displacements. It is obvious therefore that the basic formulation of one of the methods must be altered to make it compatible with the other to allow a satisfactory union of the two methods to be achieved. In the combined method, adopted for the work discussed in this paper, the boundary element equations have been modified to make them compatible with the finite element equations. This approach has the potential for wide industrial appeal as the FEM is more established than the BEM in the industrial environment and enables the combined approach to be used with minimal modifications to existing finite element software. The finite element method
The FEM for stress analysis is based on the minimization of the potential energy function which is given by the expression (e.g. Hinton and Owen9):
pE= Iv W(U)dV- IsUbdV- IsupdS
(1)
in which W(u) is the strain energy density, u is the deflections, b is the distributed body forces, p is the distributed surface forces and S and V denote surface and volume integrals respectively. On minimizing eqn (1), the following expression relating nodal forces {F} and nodal displacement {u} by the structural stiffness matrix [K] is obtained. {F} = [K]{u}
(2)
where the stiffness matrix is given in terms of the
Elasto-plastic analysis using a coupled boundary element finite element technique
the actual boundary of the problem are obtained from the RIDM at collocation points relative to the source points. Hence
/
j=p ui(xf) = Z Gu(xfxs)cpj(xs) j=i
(a)
j=p ti(xf) = Z Fij(xrxs)qSj(Xs) j=l
(b)
....
Fig. 1. Source distribution in indirect boundary element methods.
strain-displacement matrix [B] and the stress-strain matrix [D] as (3)
The regular indirect discrete method
The regular indirect discrete method (RIDM) is a development of the indirect boundary element method (IBEM). The normal formulation of the IBEM for linear elastic stress analysis is based on the use of the fundamental (Kelvin) solution 5 which satisfies the governing differential equation. For a homogeneous, isotropic domain, the relevant fundamental solutions are related to the displacements ui(xf) and the surface tractions ti(xf) at a point xf on the boundary of the components being analysed by a fictitious source q~distributed on the boundary as shown in Fig. l(a). Hence
(8)
{t} = [F]{~b}.
(9)
From eqns (8) and (9) the boundary tractions and displacements can be related by the following expression: {t} = [F][O] -l{u}.
(10)
The combined method
It is obvious, on examination of eqns (2) and (10), that the FEM cannot be linked directly with the RIBM. The latter method can however be modified to make it compatible with the FEM by using Betti's second law which is:
JvuA(u) = I v W ( u ) d V - JsUtdS
(11)
in which A is a differential operator. If eqn (11) is substituted into eqn (1) the following expression is obtained
aE=:Jvua(u)dV+:IsutS-IvubdV-JsupdS1
1
(12)
(4) In the absence of body forces eqn (12) becomes:
and
ti (Xf) : IS Fij (XfXs)~ dS.
(7)
and
(c)
Ui(Xf) = IS Gij(XfXs)O dS
(6)
in which p is the number of source points. Matrix equations for the deflections and tractions at all points on the boundary can be obtained from eqns (6) and (7). Hence
Node
{u} = [a]{~}
[rl = aIv [Blw[DI[BIdV
41
PE=~
(5)
When the field point, xf and the source point, Xs, in eqns (4) and (5) coincide the integrals become singular. However Wearing et aL6 showed that the problems associated with the singular integrals can be overcome by distributing the source on a boundary which is placed outside the actual boundary of the problem as shown in Fig. l(b). A further development 4 has resulted in the RIDM in which the sources q~j(xs) are distributed at discrete points around the source boundary as shown in Fig. l(c). By removing the integrals in eqns (4) and (5), the displacements and tractions at each point on
'Is u t d S - IsupdS.
(13)
Substituting eqn (10) into eqn (13) and minimizing the resulting expression gives the following equation:
Is[F][GI-I {u} dS = Is{q} dS
(14)
which can be written in matrix form as:
K{u}
= {e}
(15)
Equation (15) now has the same form as eqn (3) which means that the RIDM can be linked directly with the FEM. When using the combined F E M - R I D M it is convenient
J. L. Wearing,M.C. Burstow
42 to express the finite element equations as:
{F}f¢----[K]fe{u}f e
(16)
and the modified boundary element equations as: {F}be = [r]~{U}be
(17)
in which the subscripts fe and be refer to the finite element and boundary element regions respectively. When eqns (16) and (17) are combined a system of equations is obtained for the complete component. Hence {F} = [r]{u}
(18)
in which [K] is the merged stiffness matrix obtained from eqns (16) and (17) and may be written as: [K] ---- [Klfe + [K]be
(19)
and {F} and {u} are the merged force and displacement vectors.
ELASTO-PLASTIC ANALYSIS In an elasto-plastic analysis eqn (18) is solved initially to give a global elastic solution to the problem. An incremental process is then employed whereby the applied loads are incremented according to specified load factors. As the plastic behaviour is confined to the finite element region, only [K]fe requires to be recalculated in eqn (19) after each load increment and, as the boundary element region remains elastic, [K]be remains unaltered and does not therefore require to be calculated after each load increment. [K]fe can be rewritten as: [K]f e =
Iv[B]TtD]ep[B]d g
(20)
where [D]~p is the elasto-plastic stress-strain matrix and is given by the expression (e.g. Owen and Hintonm). [D]ep = [ D ] -
{do}{dD}V H ' + {dD}r{a}
(21)
in which {dn} = [Dl{a}
(22)
and {a} is the plastic flow vector given by: {a}T
0(_~] FDF OF OF OF]
(23)
In eqn (23) F(tr, K) governs the yield criterion, with K being a parameter which depends on the specific yield criterion which is used in the elasto-plastic analysis.
CASE STUDIES The following four case studies are presented to test the
(a)
(b)
Fig. 2. Finite element and boundary element for two dimensional analyses. effectiveness of the F E M - R I D M technique for the analysis of two dimensional problems with plastic zones. (a) A rectangular plate with a central circular hole. (b) A rectangular plate with symmetrical vee notches on opposite edges. (c) A centre cracked rectangular plate. (d) A four point bending specimen. For each case study, with the initiation of plastic deformation, the load was applied in increments and, since plastic deformation is confined to the finite element region, only [K]fe in eqn (19) needs to be recalculated after each load increment. As the boundary element region remains elastic throughout the analysis [K]be retains its initial values. For all case studies the finite element regions were modelled using eight noded isoparametric elements as shown in Fig. 2(a) and the boundary element regions were modelled using three noded isoparametric line elements as shown in Fig. 2(b). As the elements in each region are compatible, merging of the two regions can be achieved. The source boundary is placed outside the actual boundary of the component being analysed when the RIDM is used and previous studies, 6 on the application of the regular boundary element method for stress analysis, have indicated that there is a range of possible positions of the source boundary which yield good, stable results. These studies have revealed that when the source and component boundaries are too close or too far apart unsatisfactory results are obtained. The position of the source boundary is defined by a nondimensional factor A, which is the ratio of the distance of the source boundary from the boundary of the component to the length of the boundary element. This means that for boundary element models with elements of different lengths, the source boundary does not require to be at a fixed distance from the boundary of the component being analysed. For the case studies considered in this paper good stable results were obtained for values of ), between 1.0 and 10"0. In all analyses the value chosen for ), was 4.0. Due to the symmetrical nature of the problems discussed in the paper it was only necessary to analyse a quarter of the plate in each of the first three case studies and in the four point bending specimen the analysis was confined to half of the beam. All analyses were plane strain with a Von Mises yielding criterion and the first
Elasto-plastic analysis using a coupled boundary element finite element technique
43
--,--D.-
(b)
(a)
t
0.4~ 00
o.'s:4~
O.561
0/56 (d)
(c)
Fig. 3. Plate with central circular hole. Poisson's ratio = Uniaxial yield stress = Work hardening parameter =
three case studies had the following material properties: Young's modulus Poisson's ratio Uniaxial yield stress Work hardening parameter
= = = ---
7000 N/ram 2 0.33 24-3 N / m m 2 0-0 (perfectly plastic)
The material properties for the four point bending specimen were: Young's modulus = 200 000 N / m m 2
With respect to the fracture analyses an L E F M and E P F M analysis was undertaken for each case study. Quarter point elements were used in the finite element mesh to model the singularity at the crack tip. The results for the fracture analysis, using the combined method, were used in conjunction with a number
v,,k
.--.-I)-
,
(b)
(a)
O. 7 2 5 \
o.6
!: ....
0"28 400 N/ram 2 0"0 (perfectly plastic)
-"~-~
!. (d)
(c)
Fig. 4. Notched plate.
44
J. L. Wearing, M.C. Burstow 10mm
--
r
:
"
'
-/
/:
))
', Z . . . .
\-- ...........
-/ . . . .
I I I I t
I I 1 I I
15ram
/
/
Crack
J - Path
J - Path 8
h u'-
) . I Ill
Fig. 7. J integral paths for centre crack plate. 10ram v I
II I l l
Fig. 5. Centre cracked plate. of approaches to determine the stress intensity factors, the results of which are given in Tables 1 and 2 (later). Rectangular plate with central hole The plate has a length of 36 mm, a width of 20 mm and
_
)!
Y
l l l l l l l
the diameter of the central hole is 10ram. Details of the loading and the boundary conditions for the quarter plate model, used in the analysis, are given in Fig. 3(a). An initial tensile load which ensured an elastic stress distribution throughout the plate was applied initially as shown in Fig. 3(a). This load was scaled to a value which caused the initiation of yield at the edge of the hole and was followed by a number of load increments up to a value of 13.78 N/mm 2, which is 0"56 times the value of the yield stress of the material. The results showing the growth of the plastic zone, adjacent to the hole, are shown in Fig. 3(d) and indicate the extent of the plastic zone for each load increment, which is indicated as a percentage of the yield stress used as the applied load. The region adjacent to the hole containing the plastic region was modelled using 20 finite elements, having 130 nodes and the remainder of the plate was modelled using 16 boundary elements with a total of 32 nodes. Results were also obtained from a finite element model with 42 elements and 160 nodes. Details of the meshes are shown in Figs 3(b) and
3(c). Notched plate Details of the loading and boundary conditions for the quarter plate model are shown in Fig. 4(a). The length of the plate is 36 mm, its width is 18 mm and the depth of the 90 ° notch is 5 mm. The problem was analysed using two boundary element regions with a total of 34 Table 1. Stress intensity factors for centre cracked plate Method of calculating Ki (1)
.......
i i i i i
Fig. 6. Combined RIDM-FEM model for centre cracked plate.
(2)
LEFM analysis (a) Combined method Displacement extrapolation Strain energy release rate J-integrals (b) Rooke and Cartwrightll (c) J integrals using FEM EPFM analysis (a) J integrals using combined method (b) Theoretical plastic z o n e c o r r e c t i o n 12
Ki
2.31 2.33 2.32 2.31 2.28 34.85 34"26
Elasto-plastic analysis using a coupled boundary element finite element technique 140
45
loads being indicated as a percentage of the yield stress used as the applied load.
-
130 •
Centre cracked plate II0
9O 8O
~
50
Sl~-ess(N/ram2) ~0 lied
30
4.471
6.901
20
io 0
3.136 "~T 0.00
0.05
I
I
I
I
I
I
I
0.10
O. I S
0.20
0.~
0.30
0.35
0.40
0.486
Distance Across Ligament From Crack Tip (nma)
Fig. 8. Stress distribution across ligament of centre cracked plate. elements and a finite element region with 28 elements giving a total of 139 nodes as shown in Fig. 4(b). It was also analysed using 52 finite elements with a total of 185 nodes as shown in Fig. 4(c). The plate was initially loaded along its top edge with an elastic tensile load as shown in Fig. l(a). The load was then scaled to induce yielding at the notch and increased in six increments until the applied load was 0"725 times the yield stress of the material. Figure 4(d) shows the increase in the yield surface for each load increment with the
The dimensions and loading of the plate are shown in Fig. 5 with the combined R I D M - F E M mesh, having a total of 299 nodes, which was used to analyse the quarter plate model, being indicated in Fig. 6. Details of the finite element mesh and the J integral paths used in the LEFM and EPFM analysis are shown in Fig. 7. The problem was also analysed using a finite element mesh having 1000 nodes. For the LEFM analysis a tensile load was applied to the plate and this load was sufficiently small to ensure that the plastic zone, which must exist in the vicinity of the crack tip, was not large enough to encroach on any of the Gauss points close to the crack tip, ensuring a purely elastic analysis. The KI stress intensity factors obtained by a number of methods for the LEFM analysis are given in Table 1. For the EPFM analysis the KI values were obtained for a number of tensile loads ranging from 0-486 to 7.205N/mm 2 with the results for a tensile load of 7.205N/mm 2 being given in Table 1. The effect of increasing the load on the stress distribution along the plate ligament is shown in Fig. 8 and the growth of yield surface with increasing load is shown in Fig. 9.
Four point bending specimen As in the previous example, this case study presents the problem of a crack opening under mode I deformation.
Applied Stress = 4,471 N / r a m 2
Applied S t r e ~ = 5.686 N / r a m 2
Applied Stress = 6.598 N / r a m 2
Applied S t r ¢ ~ = 7.205 N/ram2
Fig. 9. Plastic zone growth for centre cracked plate. • FEM; - -
RIDM-FEM.
46
J. L. Wearing, M.C. Burstow 450
mm Crack
R = 0.254 mm ~ l / - -
19.05 mm 4
--I12.7
mm
d 31.75 mm -I
Fig. 10. Four point bending specimen. = :
= :
IIIFig. 11. Combined RIDM-FEM model for four point bending specimen. In this case however a bending load is applied and the crack is emanating from a stress concentration. The geometry of the specimen is shown in Fig. 10. The position of the supports ensures that the bending moment between the supports is constant. The crack extends 0"5 mm along the vertical centre line from the root of the notch. Figure 11 gives details of the combined R I D M - F E M mesh, which comprised four boundary element regions, with a total of 48 elements, and a finite element region for the plastic zone, with a total of 90 elements, giving a total of 361 nodes for the analysis. A finite element model comprising 180 elements with 850 nodes was also used to analyse the problem. The finite element mesh for the combined method and the J integral paths used in the LEFM and EPFM analyses are shown in Fig. 12. A deflection of 0.1 × 10-3mm was applied at the load points of the specimen for the LEFM analysis. This was small enough to ensure that the yield surface
adjacent to the crack tip was not large enough to be detected by the Gauss points. The KI stress intensity factors for the LEFM analysis, calculated from various methods are shown in Table 2. The elasto-plastic solution to the problem was obtained by scaling the elastic solution until the yield stress of the most highly stressed Gauss point was reached and the load was then applied incrementally until the deflection at the load points was 0.07ram. Figure 13 shows the growth of the plastic zone with increasing deflection and the distribution of the effect Von Mises' stress through the plate below the crack tip is shown in Fig. 14. DISCUSSION OF RESULTS The effectiveness of combining the finite element method and the regular indirect discrete method for the analysis of two dimensional elasto-plastic problems has been investigated using four case studies. Figures 3(d) and 4(d) show the results for the growth of the
/ J - Path 7
~__~J-Path 1
Table 2. Stress intensity factors for four point bending problem Method of calculating KI (1)
(2)
LEFM analysis (a) Combinedmethod Displacement extrapolation J integrals (b) Rooke and Cartwrightll (c) J integrals using FEM EPFM analysis (a) J integrals using combined method (b) Theoretical plastic zone correction12
KI
1' 15 1.17 1.16 1.17 821.5 816.7
/ ......5 ........L .......i Fig. 12. J integral paths for four point bending spcomen.
Elasto-plastic analysis using a coupled boundary element finite element technique
J/ /
/
/
AppliedDisplacement= 0.035 mm
/
47
/
AppliedDisplacement= 0.045 mm
AppliedDisplacement= 0.07 mm
Fig. 13. Plastic zone growth rate for four point bending specimen. • FEM; - -
plastic zone with increasing tensile load for the plane strain analysis of a rectangular plate with a central circular hole and a rectangular plate with symmetrical vee notches on opposite edges. The plates have been analysed using both the combined method and the FEM and the results from both analyses are in close correspondence. Due to the reduction in the number of degrees of freedom in the combined method compared to the FEM the analysis time using the R I D M - F E M approach was 30% less than the time taken for the FEM analysis. The results for the stress intensity factors from the LEFM and EPFM analyses are given in Tables 1 and 2. The results from the R I D M - F E M analyses have been obtained using a number of approaches as indicated in Tables 1 and 2. These results are compared with
RIDM-FEM.
the results from finite element analyses and from other sources and close agreement has been achieved. Figures 9 and 13 show the growth of the plastic zones, with increasing load, for the centre cracked plate and the four point bend specimen respectively and in each case good agreement with the results obtained from finite element analyses has been achieved. Considerable savings in model size were achieved when using the combined method compared to the FEM resulting in improved efficiency at the calculational phase. The results shown in Fig. 8 for the variation in stress across the ligament of the centre cracked plate show a smooth transition in the results between the boundary element region and the finite element region showing that the FEM and the RIDM can be successfully coupled.
48
J. L. Wearing, M.C. Burstow
400350
Applied displacements (ram)
300~
/
0.07
oIlll -~ 200 1so
0.035
'~ 100
o
0"3 0
I
2
3
4
5
6
7
8
Depth Through Beam Below Crack Tip (ram)
Fig. 14. Von Mises' stress distribution below crack tip for four point bending specimen.
CONCLUSIONS The results which have been obtained for the analysis of the case studies discussed in the paper clearly indicate that the modification to the R I D M to allow it to be combined with the FEM has proved successful for the analysis of two dimensional elasto-plastic problems and overcome the major disadvantages of the BEM and the FEM when the methods are used on their own for elasto-plastic analysis. The results from the combined analyses are in excellent agreement with finite element results and with results from other sources. H'12 The models for the combined approach require considerably fewer nodes than for the finite element models, which leads to a significant reduction in the size of the stiffness matrix for the R I D M - F E M analysis, compared to the stiffness matrix for the FEM. Further savings in analysis times are achieved from the combined method for elasto-plastic analysis since the boundary element region remains elastic and the stiffness matrix for this region retains its initial values. Only the stiffness matrix for the fmite element region needs to be recalculated since the plastic zone is confined to this region. Although the case studies discussed in this paper are all two-dimensional, there is clear evidence that the combined R I D M - F E M would be particularly advantageous for three-dimensional elasto-plastic analysis, since considerable savings in model size will be achieved for the combined method compared to the FEM and removes the problems associated with internal discretization when the BEM is used for such analyses.
REFERENCES 1. Aliabadi, M. H., Cartwright, D. J., Naehring, D. W. & Daney, W. A. Stress intensity factor weight functions for cracks at holes and half plane. Proceedings of the llth International Conference on Boundary Element Methods, Vol. 3, eds C. A. Brebbia & J. J. Connor. Springer-Verlag, Berlin, 1989, pp. 83-98. 2. Stok, B. & Bukovec, B. Stress and stress intensity factor analysis in cracked elastic bars under torsion by the boundary element method. Proceedings of the llth International Conference on Boundary Element Methods, Vol. 3, ed. C. A. Brebbia & J. J. Connor. Springer-Verlag, Berlin, 1989, pp. 111-22. 3. Kelly, D. W., Mustoe, G. G. W. & Zienkiewicz, O. C. Coupling boundary elements with other numerical methods. In Developments in Boundary Element Methods, Vol. 1, eds P. K. Baunerjee & R. Butterfield. Applied Science Publishers, London, 1979. 4. Scboliield, R. P. Development of the indirect discrete boundary method and its application to three dimensional design analysis. PhD thesis, The University of Sheffield, 1986. 5. Bannerjee, P. K. & Butterfield, R. Boundary Element Methods in Engineering Science. McGraw-Hill, Chichester, 1981. 6. Wearing, J. L., Sheikh, M. A., Patterson, C. & Abdul Rahman, A. G. A regular indirect boundary element method for stress analysis. Proceedings of the 9th International Conference on Boundary Element Methods, ¢ds C. A. Brebbia, W. L. Wendland & G. Khun. Computational Mechanics Publications, Southampton and Boston, 1987, pp. 183-98. 7. Wearing, J. L., Sheikh, M. A. & Hickson, A. J. On the application of a combined finite element boundary element method for design analysis. Proceedings of the International Conference on Boundary Element Technology - BETECH 89, eds C. A. Brebbia & N. G. Zamani.
Elasto-plastic analysis using a coupled boundary element finite element technique Computational Mechanics Publications, Southampton, 1989, pp. 99-110. 8. Wearing, J. L., Sheikh, M. A. & Rahmani, O. A combined finite element boundary element approach for three dimensional stress analysis., Proceedings of the International Symposium on Boundary Elements in Mechanical and Electrical Engineering, ed. C. A. Brebbia et al. SpringerVerlag, Ikdin, 1990, pp. 31-41. 9. Hinton, E. & Owen, D. R. J. Finite Element Programming. Academic Press, London, 1977.
49
10. Owen, D. R. J. & Hinton, E. Finite Elements in Plasticity: Theory and Practice. Pineridge Press, Swansea, 1980. 1I. Rooke, D. P. & Cartwright, D. J. A Compendium of Stress Intensity Factors. Her Majesty's Stationery Ofl~c¢, London, 1976. 12. Brock, D. Elementary Engineering Fracture Mechanics. Sijhoff and Noordhoff, Holland, 1978.