- Email: [email protected]
Automatic remeshing for three-dimensional finite element simulation of welding
Computer methods in applied mechanics and engineerlng ELSEVIER
Comput.
Automatic
Methods
Appl.
Mech.
Engrg.
147 (1997)
401-409
remeshing for three-dimensional simulation of welding
finite element
L.-E. L,indgren”,*, H.-A. Hfiggblad”, J.M.J. McDillb, A.S. Oddy” “LuleB University of Technology, 971 87 Luled, Sweden bCarleton University, Ottawa, KlS 5B6, Canada ‘OddylMcDill Numerical Investigation Sciences Inc., Ottawa, Kl G OW8, Canada Received
17 January
1997; revised
31 January
1997
Abstract Three-dimensional finite element simulation of electron beam welding of a large copper canister has been performed. The use of an automatic remeshing algorithm, based on a graded hexahedral element was found to be effective. With this algorithm the strongly nonlinear thermomechanical effects locally close to the moving heat source can accurately be modelled using a dense element mesh that follows the heat source.
1. Introduction
Finite element (FE) simulation of welding has been performed for more than twenty years [l-4]. The early models were only two-dimensional (2D). The first three-dimensional (3D) models were presented around 1990 [5-71. Fig. 1 shows problem size versus year of publication for [6-131 and current paper. Despite the increase in problem size during the later years the computational task can be overwhelming when a good spatial resolution is required in the neighbourhood of the weld. The thermal and mechanical fields from the welding process have large gradients near the arc. Therefore, it is natural to apply some remesh:ing techniques as the arc is moving. This paper presents an application where automatic remeshing of a graded hexahedral element is utilised. Two models, one with and one without remeshing, that give the same accuracy were used. The application is a large copper canister that is intended to be used for storage of nuclear waste. It is made of two halves tha; are-joined by electron beam welding. The aim of-the investigation was to find residual plastic strains and stresses in order to estimate subsequent creep-damage of the canister [14,15].
2. Spatial resolutiou
required
in simulation
of welding
The use of a moving, concentrated heat source is common for most welding processes. The characteristic length can vary from less than a millimetre, in the case of laser welding, up to a couple of centimetres, in the case of a gas flame. Rykalin [16] represented this by a diagram showing thermal flux
* Corresponding
author.
0045-7825/97/$17.00 0 1997 Elsevier PIZ SOO45-7825(97)00025-X
Science
S.A.
All rights
reserved
402
L.-E.
Lindgren
et al. I Comput. Methods Appl. Mech. Engrg.
147 (1997) 401-409
lE+7
p31
lE+5
1975
1990
1995
1990
1995
moo
YEW Fig. 1. Size, number of unknowns in mechanical analysis times the number of time steps, of computational of welding. The reference numbers are given in the diagram, * denotes this study.
models
for simulation
and source size. This localisation creates a large thermal gradient with corresponding gradients in the mechanical fields. However, this is not the sole factor determining the required minimum size of the finite elements. It is also determined by the scope of the investigation. If the effect of the residual stresses and deformations on the subsequent use of the component is the focus, then it is possible to use a coarser mesh than if hot cracking is studied. It may be enough to have a few elements in the heat-effected zone (HAZ) in the first case. However, it may turn out that many elements are needed, requiring element sizes as small as 50 microns as in some laser welding applications [17]. On the other hand, the computational model must be large enough to give correct boundary conditions for the thermal and mechanical fields. The latter is more demanding due to the elliptic nature of the quasistatic thermal stress problem. The combination of a model including small elements that are part of a finite element model of a large structure is further complicated by the motion of the heat source. The domain with large gradients is moving. The obvious solution to this is to have a mesh with a fine region that is moving with the heat source. This basic approach was utilised by Rosenthal [18] who used a coordinate system fixed relative to the heat source. A co-moving frame applied in a finite element analysis of a process with moving heat source was used by Bergheau et al. [19]. Brown et al. [20] used dynamic substructuring where a denser mesh was travelling with the heat source on an overall coarser mesh. Linear hexahedral elements are superior to linear tetrahedral elements [21,22]. They are also better than quadratic tetrahedron elements when plastic deformation occurs [22]. However, most analyses involving remeshing use tetrahedron elements. This is due to the good availability of procedures for generating arbitrary graded meshes using tetrahedral elements. The models in this work are based on a graded hexahedron, described in the next section, that alleviates the creation of a mesh with a refined region. The element is combined with a remeshing strategy that makes it possible to move the refined region with the heat source.
3. Remeshing of graded hexahedron It is well known that h-adaptive techniques are preferable wherever the solution is rough; i.e., the gradients are not continuous [23]. The h-adaptive strategy used here follows the techniques developed by McDill and co-workers [24-271. The method has three components, the basic grading hexahedral
L.-E.
Lindgren
et al. I Comput. Methodc Appl. Mech. Engrg.
147 (1997) 401-409
403
elements [24], the adaptive remeshing technique [25,26], and the data transfer between meshes [27]. The h-adaptive strategy was selected because of the nature of the temperature, stress and microstructure fields involved in welding. Although continuous, these fields may not always vary smoothly because of the temperature dependence of material properties as well as internal boundaries caused by phase transformations and plasticity. These internal boundaries are physically manifested as HAZ boundaries and plastic zone boundaries. In addition, singularities such as crack tips, re-entrant corners and point loads are better captured with h-adaptive schemes unless special singularity elements are used. 3.1. Grading hexahedron
elements
The technique is built around the 8-26 node grading hexahedron originally described in [24]. This is a variation of the familiar linear element in that temperatures or displacements vary linearly with position within each element. It can be extended to a quadratic, non-conforming grading hexahedron as in [28] for use in welding analyses [13] but this form has not been used in this investigation. The structure of the element requires eight mandatory corner nodes (see Fig. 2). Nodes 9 to 20 are optional, located at the mid-edges and nodes 21 to 26 are optional mid-face nodes. These optional nodes may be included in any combination. The presence of a mid-face node modifies the basis functions of existing mid-edge nodes and corner nodes present on that face. Mid-edge nodes modify the basis functions of the corner nodes on that edge. This preserves interelement continuity. The basis functions are shown in Table 1. The table should be evaluated, in decreasing order, from Nz6 to N,. If a node is not present, its basis function and all subsequent references to it should be removed from the table. The element is divided into, at most eight, linear subdomains, depending on the arrangement of mid-edge and mid-face nodes. Numerical integration occurs over each subdomain. The net result is very similar to the assembly of a collection of standard linear elements with some degrees-of-freedom condensed out. However, by embedding the constraints in the basis functions, all nodes are active. Global constraint equations, matrix condensation operations and additional data structures to identify redundant nodes are not required. These elements have proven to be computationally very effective, especially around point loads and weak singularities such as welding heat sources. As described in [24], a point thermal load in an infinite half space has a singularity proportional to the inverse of the radial distance. For a cubic domain with the point load at one corner, if 12is the ratio of the size of the cube to the size of the smallest element then the number of elements in a standard, uniform mesh is n3. In a recursively graded mesh only 7 *log,(n), that is 7 times logarithm base 2 of IZ, elements are required to obtain elements of a similar size near the singularity. This grading has been the key to successful 3D thermal-mechanical finite element analysis of welds since it allows an easy concentration of elements and nodes around the welding source where gradients change rapidly. 3.2. Adaptive remeshing technique The structure of the grading hexahedrons led to an adaptive remeshing strategy as described in [25,26]. It allows easy modifications to be made in the mesh by adding extra nodes and subdividing elements only in regions where it is required. Element division into octants preserves shape regularity
Fig. 2. The eight
to twenty-six
node,
grading
hexahedron
element.
404
L.-E.
Lindgren
et al. I Comput.
Methods
Appl.
Mech.
Engrg.
147 (1997) 401-409
Table 1 Basis functions
26 25 24 23 22 21 20 19 18 17 16 15 14 13 12 11 10 9 8 7 6 4 3 2
0.5(1- kl)U - lbm- 0 O-5(1- kl>(1- /d)(1 + 5) w1 - 121)(1-17)(1Irl1
031 - ml +7fx1- Irl> 0.5(1 - 5x1 - lM1 - 14% 0.5(1+ S)(llb7lN- lil) o.q1+ 5X1-71)(1- IiD - 0.5(&
+ K,)
0.25(1-t)(l-q)(l-Iii)-0.5(&z+&) 0.25(1- <)(I +v)(l - kl) - 0.5@‘,, + k) 0.25(1 + E>(l +rl)(l - I&/) - 0.5(&, + f%) 0.25(1 + 5)(1 - b11)(1- 5) - 0.5(N,, + N,,) 0.25( - IS[)(I - v)(l - J-) - 0.5(%, + W 0.25(1- [)(IIvl)(l - %I- O.W, + N,,) 0.25(I - kl)(l+~)(l0 -0.5(&x + N,,) Q.25(1 + 5)(1b11>(1+ l) - 0.5(&, + KS) 0.25(- [[I)(1 -77)(1+5)-0.5(N,,+N,,) 0.25(1- 5)(1- hl)(l + 5) - 0.5(N,, + N,,) 0.25(f - ISl)(l +d(l + 5) - 0.5(N,, + K,) 0.125(1+ 5)(1 -v)(l -i) - 0.25(N,, + Nz4 + N,,) 0.125( 1 - .$)(l - v)(l - J) - 0.25(N,, + A’,, •t N,,) 0.125(1 - 5)(1 + v)(l - 5) - 0.25(N,, + Nz3 + N,,) 0.125( 1+ t)( 1 + v)( 1 - 5) - 0.25@‘,, + h’z3 + N,,) 0.125(1+5)(1 -v)(l + 6) -0.25@‘,, + Nu +N& 0.125( 1 - .$)( 1 - q)( 1 + i) - 0.25(&, + Nz4 + N,,) 0.125(1 - t)(l + TJ)(~+ 5) - 0.25(N,, + h’>x+ N,,) 0.125(1 + [)(l +r,)(l + r) -0.25(&, + N,, f IV,,)
-0.5(N,, - 0-5(N,, - 0.5@‘,, - O-5@‘,, -0.5(N,, - 0.5(N,, - 0.5(N, -O-5@&
+N,6 + Nm) + N,, + N,,) + N,, + N,d + N,, + N,,) +N,, +Nx,) + N,, + N,,) + N,, + N18) + N,, + N,,)
but is not an essential requirement. This retains the basic mesh and changes are localized. Simple array type data structures identify the elements attached to corner, mid-edge and mid-face nodes are all that are required for automatic refinement. Coarsening; i.e. the replacement of many small elements with fewer large elements is also possible. The element assembly technique in [26] eliminates the need to store or retrieve the former topology. Either fully arbitrary [26] or tree-type coarsening strategies are possible. In this application coarsening returns the mesh to a former configuration by retracing the refinement path. A tree structure which identifies the central nodes at each refinement level is used to retrace the path. The rules and operations followed in refinement and coarsening are described in [25,26]. Briefly, bisection refinement and binary grading preserve shape regularity and avoid abrupt transitions from coarse to fine regions. Refinement propagates to neighbouring elements. Coarsening is not permitted to propagate to neighbours and only occurs if all members of a local group allow it and if it results in a valid mesh. Only redundant, i.e. unnecessary, nodes are eliminated. The creation of new meshes with refinement in some areas and coarsening in others changes the topology of the mesh. This changes the number and locations of nodes and Gaussian integration points. In time dependent problems information associated with nodes and Gauss points must be transferred from the old mesh to the new mesh. If this is not done with some care it is quite possible that the transfer will introduce errors in the temperature, stress and plastic strain fields. This changes the total energy in the system which may be thermal energy in the form of temperatures, elastic energy as stresses, or energy dissipated as plastic work. The data transfer from the old to the new mesh depends on whether there is coarsening or refining of the mesh and if it is nodal data or Gauss point data. Interpolation is used for nodal data and Gauss point data in refinement. The coarsening procedure is somewhat more complex than the refining procedure as the coarser mesh can only approximately represent the old, finer mesh. It was decided to use a least-square method based on a local domain for transfer of nodal data. Gauss point data were transferred by interpolation. See [27] for details.
L.-E.
3.3.
Lindgren
Measures for determining
et al. I Comput. Methods Appl. Mech. Engrg.
147 (1997) 401-409
405
element sizes
In this investigation a simple geometric measure; h/r where h is the size of the element and r is the shortest distance from the element centre to the path of the heat source until next remeshing is performed, was used to drive the refinement and coarsening. This makes it possible to refine in advance of the motion of the heat source. The user can then choose an optimal combination of the amount of refinement /coarsening each time and how often this should be performed. Refinement or coarsening decisions are ideally driven by error measures based on the solution and the thermal, stress and microstructural fields present. In the future, generic error measures such as [29] which can evaluate the error in the temperature, stress, plastic strain or microstructure fields will be used to determine which elements to refine or coarsen. A posteriori measures of this type can only change the mesh once a trial solution is available. With a moving heat source, as is found in welding analysis, it would be a distinct advantage to alter the mesh in response to the source motion before generating a solution at each step.
4. Longitudinal
butt welding of a copper canister
4.1. Background In the Swedish nuclear waste program it has been proposed that spent nuclear fuel shall be placed in composite copper-steel canisters. These canisters will be placed in holes in tunnels located some 500 m underground in a rock storage. The canister consists of two cylinders of 5 m length, one inner cylinder made of steel and one outer cylinder made of copper. The outer diameter of the canister is 500 mm and the wall thickness is 50 mm. For storage, the steel cylinder, which contains the spent nuclear fuel, is placed inside the copper cylinder. Thereafter, a copper end is butt welded to the copper cylinder using electron beam weld:ing. The primary objective in the mechanical design of the canister is to ensure no leakage of radioactive particles to the surroundings. One of the alternative fabrication schemes for the outer copper cylinder is based on cold forming (rolling) and subsequent longitudinal seam welding of the two copper half cylinders. 4.2. Computational models The potential of the proposed algorithm is demonstrated by studying the longitudinal butt welding of two copper half cylinders previously studied in [15]. The models used do not include residual stress from the rolling process. The transient and residual temperature, stress and strain fields present in the canister during longitudinal seam welding and after cooling to ambient temperature are calculated numerically by use of the coupled thermo-mechanical FE-code SIMPLE [30]. Due to symmetry the FE-model consists of one half of a copper half cylinder. The models extend 500 mm in the axial direction including the rear end where the welding finishes. The cylinder is restrained axially at the front end, the rear lend is free. In order to have a satisfying accuracy in the simulations, the elements close to the weld-path must be around 5 mm in size, transverse to the welding direction. Two models are used: one with and the other without remeshing. The meshes were chosen in order to achieve about the same accuracy. The thermal material data for the copper alloy were taken from handbooks (see [8] for the data used). A thermo-elastic-plastic model with linear isotropic hardening is used in the mechanical analysis. The material data for the mechanical analysis were determined from uniaxial tensile tests. See [15] for the elastic and plastic material data used. The load is applied by defining a moving box along the cylinder. In each loading-step the nodes in the box and those passed by the box are subjected to a uniform input of energy. The net heat power is 8 kW and the welding speed is 5 mm/s. The loading is applied during the first 100 s of the simulation which is performed in a total of 1000 s. Both calculations are performed in 324 loading-steps. The FE-model without remeshing consists of 5783 nodes and 4220 eight-node brick elements. To reduce the number of elements in the model without reducing the accuracy a dense mesh is required in the vicinity of the moving heat source but a
406
L.-E.
Lindgren
et al. I Comput. Methods Appl. Mech. Engrg.
147 (1997) 401-409
coarser mesh is used elsewhere. See Fig. 3(b) for the FE-mesh of the non-remeshing model. The FE-model used for remeshing consists initially of 1000 eight-node brick elements (1683 nodes) but increases to maximum 2715 elements (441 nodes) at 70 s. The remeshing is performed nine times. See Fig. 3(a) for the FE-mesh at different stages of the calculation. 4.3. Computed temperatures and stresses The two FE-calculations, without and with remeshing, were performed on an IBM RS6000, model 590. The calculation without remeshing required 34 h and 20 min. of CPU-time. The CPU-time for the model with remeshing was 14 h and 5 min. The size of the stiffness and the conductivity matrix of the non-remeshing model totals 9.7 million terms. The corresponding size for the remeshing model varies from 1.6 to 5.3 million. The evolution of the calculated axial stress is shown in Fig. 3. The plastic strain fields at the final stage are shown in Fig. 4. Comparison of the residual stress-states is shown in Fig. 5. The temperature at 50 s (when the beam is half-way along the cylinder) is shown in Fig. 6. The
Fig. 3. Comparison of the axial stress field (z-stress). (a) With remeshing at 10, 50, 100 and 200 s; (b) without remeshing at 50 and 100 s.
L:,E.
Lindgren
mill
I
et al. I Comput. Methods Appl. Mech. Engrg.
0.0
nnx=o.o3e
147 (1997) 401-409
407
mh=o.o mLx=O.OU
Y
a)
W
Fig. 4. Comparison of the residual effective plastic strain. (a) With remeshing; (b) without remeshing. Fringe Iev6l6 [MPal
60 75 69 64 59 53 48
43 37 32 27 21 16
min.
1.95 MPrI
llWXE61.8MPB
mm * 1.65 f4Pe mar Ia%5 MP.¶
11 5 0 15
a)
W
Fig. 5. Comparison of the residual effective stress (von Mises). (a) With remeshing; (b) without remeshing.
maximum temperature at the nodes in the centre of the weld pool reach approximatively temperature, 1083°C. However, the used postprocessor plots element average temperatures a maximum temperature of 520°C in Fig. 6.
4.4. Comparisons
the melting resulting in
of computational models
As can be seen in Figs. 3 to 5, the simulation using remeshing gives similar results to the simulation without remeshing. The reduction in size of the models decreases the required CPU-time. It also reduces the data-storage requirement substantially, this can in some applications be the key factor determining whether a successful FE-simulation can be completed or not.
L.-E.
408
Lindgren
et al. I Comput. Methods Appl. Mech. Engrg.
Fig. 6. Temperature
distribution
147 (1997) 401-409
at 50 s.
5. Discussion and conclusions The automatic remeshing procedure is found to be effective. The required CPU-time was reduced by 60% for the studied application without any loss in accuracy. Larger gains will be expected for application with longer welds as the zone associated with large gradients will be smaller relative to the total length of the weld. The current logic for deciding the size of the elements in relation to the distance to the arc will be combined with a posteriori error measure, in order to have a versatile tool that includes a predictive capability.
Acknowledgements This work was funded by the Swedish Nuclear Fuel and Waste Management Sweden and by Lulei University of Technology, Lulea, Sweden.
Co (SKB) in Stockholm,
References [II L. Karisson,
PI
[31 [41 PI PI [71 PI
Thermal stresses in welding, in: R. Hetnarski, ed., Thermal Stresses, Ch. 5 (Elsevier, Amsterdam, 1986). J. Goldak, A. Oddy, M. Gu, W. Ma, A. Mashaie and E. Hughes, Coupling heat transfer, microstructure evolution and thermal stress analysis in weld mechanics, in: L. Karlsson, M. Johnson and L.-E. Lindgren, eds., Proc. of IUTAM Symp. Mechanical Effects of Welding (Springer-Verlag, Berlin, 1992) l-30. D. Radaj, Heat Effects of Welding (Springer-Verlag, Berlin, 1992). Y. Ueda, M. Murakawa, K. Nakacho and N.-X. Ma, Establishment of computational weld mechanics, Trans. JWRI 24 (1995) 73-86. A.S. Oddy, J.A. Goldak and J.M.J. McDill, Transformation effects in the 3D finite element analysis of welds, in: Proc. of 2nd Inter. Conf. on Trends in Welding Research, Gatlinburg, USA, 1989). R.I. Karlsson and L. Josefson, Three-dimensional finite element analysis of temperatures and stresses in single-pass butt-welded pipe, Trans. ASME, J. Press. Vessel & Tech. 112 (1990) 76-84. L.-E. Lindgren and L. Karlsson, Deformations and stresses in welding of shell structures, Int. J. Numer. Methods Engrg. 25 (1988) 635-655. E. Friedman, Thermomechanical analysis of the welding process using the finite element method, Trans. ASME, J. Pressure Vessel Tech. (1975) 206-213.
L.-E.
Lindgren
et al. I Comput. Methods Appl. Mech. Engrg.
147 (1997) 401-409
409
[9] B. Andersson and I,. Karlsson, Thermal stresses in large butt-welded plates, J. Thermal Stresses 4 (1981) 491-500. [lo] M. Jonsson, L. Karlsson and L.-E. Lindgren, Deformations and stresses in butt-welding of large plates, in: R.W. Lewis, ed., Numerical Methods in Heat Transfer, III (1985) 35-38. [ll] A. Oddy, J. Goldak and M. McDilI, Transformation plasticity and residual stresses in single-pass repair welds, Trans. ASME, J. Pressure Vessel Tech. 11 (1992) 33-38. [12] J. Goldak, A. Oddy and D. Dorling, Finite element analysis of welding on fluid filled pressurized pipelines, in: Proc. 3rd Int. Conf. on Trends in Welding Research, Gathnburg, USA, 1992, 45-50. [13] A. Oddy and J.M. E&Dill, Residual stresses in weaved repair welds, in: Proc. 1996 ASME Pressure Vessel & Piping Conf., Montreal, Canada, 1966, 147-152. [14] B.L. Josefson, L. Karlsson, L.-E. Lindgren and M. Jonsson, Thermo-mechanical FE-analysis of butt-welding of Cu-Fe canister for spent fu’el, in: K.F. Kussmaul, ed., Trans. 12th Int. Conf. on Structural Mechanics in Reactor Technology, paper NO115 in Vol. N (Stuttgart, Germany, 1993). [15] B.L. Josefson, L.-E.. Lindgren, H.-A. Haggblad and L. Karlsson, Thermo-mechanical FE-analysis of the fabrication of a Cu-Fe canister for spent nuclear fuel, in: S.V. Mbller, ed., Trans. 13th Int. Conf. on Structural Mechanics in Reactor Technology (SMiRT-13), Universidade Federal do Rio Grande do Sul, Brazil, 1995, 275-280. [16] N.N. Rykalin, Energy sources used in welding, Soudage et Techniques Connexes l(12) (1974) 471-485. [17] H. Runnemalm, L...E. Lindgren, M.O. Nisstriim and C. Lampa, Accuracy in thermal analysis of laser welding, in: G.M. Carlomagno and C.A. Brebbia, eds., Proc. Comput. Methods and Experimental Measurements VII (Computational Mechanics Publications, Southampton, 1995) 85-92. [18] D. Rosenthal, Mathematical theory of heat distribution during welding and cutting, Welding Journal, welding research supplement (1941) 220-234. [19] J.M. Bergheau, D. Pont and J.B. Leblond, Three-dimensional simulation of a laser surface treatment through steady-state computation in the heat source’s comoving frame, in: L. Karlsson, M. Johnson and L.-E. Lindgren, eds., Proc. IUTAM Symp. Mech. Effects of Welding (Springer-Verlag, Berlin, 1992) 85-92. [20] S. Brown, K. Christie and H. Song, Three dimensional finite element modeling of welded structures, in: L. Karlsson, M. Johnson and L.-E. Lindgren, eds., Proc. IUTAM Symp. Mech. Effects of Welding (Springer-Verlag, Berlin, 1992) 181-188. [21] A.O. Cifuentes and A. Kalbag, A performance study of tetrahedral and hexahedral elements in 3-D finite element structural analysis, Fin. Elem. Anal. Des. 12 (1992) 313-318. [22] S.E. Benzley, E. Perry, K. Merkely, B. Clark and G.D. Sjaardama, A comparison of all hexagonal and all tetrahedral finite element meshes for elastic and elasto-plastic analysis, in: Proc. 14th Ann. Int. Meshing Roundtable, Albuquerque, USA, 1995. [23] B.A. Szabo, Estimation and control of error based on P-convergence, in: Proc. Int. Conf. on Accuracy Estimates and Adaptive Refinements in Finite Element Computations (ARFEC), Lisbon, Portugal, 1984. [24] J.M. McDill, J.A. Goldak, A.S. Oddy and M.J. Bibby, Isoparametric quadrilaterals and hexahedrons for mesh grading algorithms, Comm. Appl. Numer. Methods 3 (1987) 155-163. [25] J.M.J. McDill, AS. Oddy and J.A. Goldak, An adaptive mesh-management algorithm for three-dimensional automatic finite element analysis, Trans. CSME 15 (1991) 57-70. [26] J.M.J. McDill and .4.S. Oddy, Arbitrary coarsening for adaptive mesh management in automatic finite element analysis, J. Math. Modelling Sci. Comput. 2B (1993) 1072-1077. [27] J.M.J. McDill, ME. Klein and A.S. Oddy, Data transfer for 3D h-adaptive thermal-elasto-plastic finite element analysis, simulations of materials processing: theory, methods and applications, in: Proc. Numiform 95, Ithaca, New York, 1995, 463-468. [28] J.M.J. McDill and A.S. Oddy, A non-conforming eight to 26-node hexahedron for three-dimensional thermal-elastoplastic finite element analysis, Comput. Struct. 54 (1995) 183-189. [29] J.M.J. McDill, M.A. Johnston and A.S. Oddy, A generic error estimate for three-dimensional h-adaptive finite element analysis, submitted to the Eleventh Int. Conf. on Mathematical Modelling and Scientific Computing, Sept. 1996. [3O] L.-E. Lindgren, User’s manual for SIMPLE-Finite Element Program for Simulation of Material Processing, Lulei University of Technology, Sweden, 1996.
Comput.
Automatic
Methods
Appl.
Mech.
Engrg.
147 (1997)
401-409
remeshing for three-dimensional simulation of welding
finite element
L.-E. L,indgren”,*, H.-A. Hfiggblad”, J.M.J. McDillb, A.S. Oddy” “LuleB University of Technology, 971 87 Luled, Sweden bCarleton University, Ottawa, KlS 5B6, Canada ‘OddylMcDill Numerical Investigation Sciences Inc., Ottawa, Kl G OW8, Canada Received
17 January
1997; revised
31 January
1997
Abstract Three-dimensional finite element simulation of electron beam welding of a large copper canister has been performed. The use of an automatic remeshing algorithm, based on a graded hexahedral element was found to be effective. With this algorithm the strongly nonlinear thermomechanical effects locally close to the moving heat source can accurately be modelled using a dense element mesh that follows the heat source.
1. Introduction
Finite element (FE) simulation of welding has been performed for more than twenty years [l-4]. The early models were only two-dimensional (2D). The first three-dimensional (3D) models were presented around 1990 [5-71. Fig. 1 shows problem size versus year of publication for [6-131 and current paper. Despite the increase in problem size during the later years the computational task can be overwhelming when a good spatial resolution is required in the neighbourhood of the weld. The thermal and mechanical fields from the welding process have large gradients near the arc. Therefore, it is natural to apply some remesh:ing techniques as the arc is moving. This paper presents an application where automatic remeshing of a graded hexahedral element is utilised. Two models, one with and one without remeshing, that give the same accuracy were used. The application is a large copper canister that is intended to be used for storage of nuclear waste. It is made of two halves tha; are-joined by electron beam welding. The aim of-the investigation was to find residual plastic strains and stresses in order to estimate subsequent creep-damage of the canister [14,15].
2. Spatial resolutiou
required
in simulation
of welding
The use of a moving, concentrated heat source is common for most welding processes. The characteristic length can vary from less than a millimetre, in the case of laser welding, up to a couple of centimetres, in the case of a gas flame. Rykalin [16] represented this by a diagram showing thermal flux
* Corresponding
author.
0045-7825/97/$17.00 0 1997 Elsevier PIZ SOO45-7825(97)00025-X
Science
S.A.
All rights
reserved
402
L.-E.
Lindgren
et al. I Comput. Methods Appl. Mech. Engrg.
147 (1997) 401-409
lE+7
p31
lE+5
1975
1990
1995
1990
1995
moo
YEW Fig. 1. Size, number of unknowns in mechanical analysis times the number of time steps, of computational of welding. The reference numbers are given in the diagram, * denotes this study.
models
for simulation
and source size. This localisation creates a large thermal gradient with corresponding gradients in the mechanical fields. However, this is not the sole factor determining the required minimum size of the finite elements. It is also determined by the scope of the investigation. If the effect of the residual stresses and deformations on the subsequent use of the component is the focus, then it is possible to use a coarser mesh than if hot cracking is studied. It may be enough to have a few elements in the heat-effected zone (HAZ) in the first case. However, it may turn out that many elements are needed, requiring element sizes as small as 50 microns as in some laser welding applications [17]. On the other hand, the computational model must be large enough to give correct boundary conditions for the thermal and mechanical fields. The latter is more demanding due to the elliptic nature of the quasistatic thermal stress problem. The combination of a model including small elements that are part of a finite element model of a large structure is further complicated by the motion of the heat source. The domain with large gradients is moving. The obvious solution to this is to have a mesh with a fine region that is moving with the heat source. This basic approach was utilised by Rosenthal [18] who used a coordinate system fixed relative to the heat source. A co-moving frame applied in a finite element analysis of a process with moving heat source was used by Bergheau et al. [19]. Brown et al. [20] used dynamic substructuring where a denser mesh was travelling with the heat source on an overall coarser mesh. Linear hexahedral elements are superior to linear tetrahedral elements [21,22]. They are also better than quadratic tetrahedron elements when plastic deformation occurs [22]. However, most analyses involving remeshing use tetrahedron elements. This is due to the good availability of procedures for generating arbitrary graded meshes using tetrahedral elements. The models in this work are based on a graded hexahedron, described in the next section, that alleviates the creation of a mesh with a refined region. The element is combined with a remeshing strategy that makes it possible to move the refined region with the heat source.
3. Remeshing of graded hexahedron It is well known that h-adaptive techniques are preferable wherever the solution is rough; i.e., the gradients are not continuous [23]. The h-adaptive strategy used here follows the techniques developed by McDill and co-workers [24-271. The method has three components, the basic grading hexahedral
L.-E.
Lindgren
et al. I Comput. Methodc Appl. Mech. Engrg.
147 (1997) 401-409
403
elements [24], the adaptive remeshing technique [25,26], and the data transfer between meshes [27]. The h-adaptive strategy was selected because of the nature of the temperature, stress and microstructure fields involved in welding. Although continuous, these fields may not always vary smoothly because of the temperature dependence of material properties as well as internal boundaries caused by phase transformations and plasticity. These internal boundaries are physically manifested as HAZ boundaries and plastic zone boundaries. In addition, singularities such as crack tips, re-entrant corners and point loads are better captured with h-adaptive schemes unless special singularity elements are used. 3.1. Grading hexahedron
elements
The technique is built around the 8-26 node grading hexahedron originally described in [24]. This is a variation of the familiar linear element in that temperatures or displacements vary linearly with position within each element. It can be extended to a quadratic, non-conforming grading hexahedron as in [28] for use in welding analyses [13] but this form has not been used in this investigation. The structure of the element requires eight mandatory corner nodes (see Fig. 2). Nodes 9 to 20 are optional, located at the mid-edges and nodes 21 to 26 are optional mid-face nodes. These optional nodes may be included in any combination. The presence of a mid-face node modifies the basis functions of existing mid-edge nodes and corner nodes present on that face. Mid-edge nodes modify the basis functions of the corner nodes on that edge. This preserves interelement continuity. The basis functions are shown in Table 1. The table should be evaluated, in decreasing order, from Nz6 to N,. If a node is not present, its basis function and all subsequent references to it should be removed from the table. The element is divided into, at most eight, linear subdomains, depending on the arrangement of mid-edge and mid-face nodes. Numerical integration occurs over each subdomain. The net result is very similar to the assembly of a collection of standard linear elements with some degrees-of-freedom condensed out. However, by embedding the constraints in the basis functions, all nodes are active. Global constraint equations, matrix condensation operations and additional data structures to identify redundant nodes are not required. These elements have proven to be computationally very effective, especially around point loads and weak singularities such as welding heat sources. As described in [24], a point thermal load in an infinite half space has a singularity proportional to the inverse of the radial distance. For a cubic domain with the point load at one corner, if 12is the ratio of the size of the cube to the size of the smallest element then the number of elements in a standard, uniform mesh is n3. In a recursively graded mesh only 7 *log,(n), that is 7 times logarithm base 2 of IZ, elements are required to obtain elements of a similar size near the singularity. This grading has been the key to successful 3D thermal-mechanical finite element analysis of welds since it allows an easy concentration of elements and nodes around the welding source where gradients change rapidly. 3.2. Adaptive remeshing technique The structure of the grading hexahedrons led to an adaptive remeshing strategy as described in [25,26]. It allows easy modifications to be made in the mesh by adding extra nodes and subdividing elements only in regions where it is required. Element division into octants preserves shape regularity
Fig. 2. The eight
to twenty-six
node,
grading
hexahedron
element.
404
L.-E.
Lindgren
et al. I Comput.
Methods
Appl.
Mech.
Engrg.
147 (1997) 401-409
Table 1 Basis functions
26 25 24 23 22 21 20 19 18 17 16 15 14 13 12 11 10 9 8 7 6 4 3 2
0.5(1- kl)U - lbm- 0 O-5(1- kl>(1- /d)(1 + 5) w1 - 121)(1-17)(1Irl1
031 - ml +7fx1- Irl> 0.5(1 - 5x1 - lM1 - 14% 0.5(1+ S)(llb7lN- lil) o.q1+ 5X1-71)(1- IiD - 0.5(&
+ K,)
0.25(1-t)(l-q)(l-Iii)-0.5(&z+&) 0.25(1- <)(I +v)(l - kl) - 0.5@‘,, + k) 0.25(1 + E>(l +rl)(l - I&/) - 0.5(&, + f%) 0.25(1 + 5)(1 - b11)(1- 5) - 0.5(N,, + N,,) 0.25( - IS[)(I - v)(l - J-) - 0.5(%, + W 0.25(1- [)(IIvl)(l - %I- O.W, + N,,) 0.25(I - kl)(l+~)(l0 -0.5(&x + N,,) Q.25(1 + 5)(1b11>(1+ l) - 0.5(&, + KS) 0.25(- [[I)(1 -77)(1+5)-0.5(N,,+N,,) 0.25(1- 5)(1- hl)(l + 5) - 0.5(N,, + N,,) 0.25(f - ISl)(l +d(l + 5) - 0.5(N,, + K,) 0.125(1+ 5)(1 -v)(l -i) - 0.25(N,, + Nz4 + N,,) 0.125( 1 - .$)(l - v)(l - J) - 0.25(N,, + A’,, •t N,,) 0.125(1 - 5)(1 + v)(l - 5) - 0.25(N,, + Nz3 + N,,) 0.125( 1+ t)( 1 + v)( 1 - 5) - 0.25@‘,, + h’z3 + N,,) 0.125(1+5)(1 -v)(l + 6) -0.25@‘,, + Nu +N& 0.125( 1 - .$)( 1 - q)( 1 + i) - 0.25(&, + Nz4 + N,,) 0.125(1 - t)(l + TJ)(~+ 5) - 0.25(N,, + h’>x+ N,,) 0.125(1 + [)(l +r,)(l + r) -0.25(&, + N,, f IV,,)
-0.5(N,, - 0-5(N,, - 0.5@‘,, - O-5@‘,, -0.5(N,, - 0.5(N,, - 0.5(N, -O-5@&
+N,6 + Nm) + N,, + N,,) + N,, + N,d + N,, + N,,) +N,, +Nx,) + N,, + N,,) + N,, + N18) + N,, + N,,)
but is not an essential requirement. This retains the basic mesh and changes are localized. Simple array type data structures identify the elements attached to corner, mid-edge and mid-face nodes are all that are required for automatic refinement. Coarsening; i.e. the replacement of many small elements with fewer large elements is also possible. The element assembly technique in [26] eliminates the need to store or retrieve the former topology. Either fully arbitrary [26] or tree-type coarsening strategies are possible. In this application coarsening returns the mesh to a former configuration by retracing the refinement path. A tree structure which identifies the central nodes at each refinement level is used to retrace the path. The rules and operations followed in refinement and coarsening are described in [25,26]. Briefly, bisection refinement and binary grading preserve shape regularity and avoid abrupt transitions from coarse to fine regions. Refinement propagates to neighbouring elements. Coarsening is not permitted to propagate to neighbours and only occurs if all members of a local group allow it and if it results in a valid mesh. Only redundant, i.e. unnecessary, nodes are eliminated. The creation of new meshes with refinement in some areas and coarsening in others changes the topology of the mesh. This changes the number and locations of nodes and Gaussian integration points. In time dependent problems information associated with nodes and Gauss points must be transferred from the old mesh to the new mesh. If this is not done with some care it is quite possible that the transfer will introduce errors in the temperature, stress and plastic strain fields. This changes the total energy in the system which may be thermal energy in the form of temperatures, elastic energy as stresses, or energy dissipated as plastic work. The data transfer from the old to the new mesh depends on whether there is coarsening or refining of the mesh and if it is nodal data or Gauss point data. Interpolation is used for nodal data and Gauss point data in refinement. The coarsening procedure is somewhat more complex than the refining procedure as the coarser mesh can only approximately represent the old, finer mesh. It was decided to use a least-square method based on a local domain for transfer of nodal data. Gauss point data were transferred by interpolation. See [27] for details.
L.-E.
3.3.
Lindgren
Measures for determining
et al. I Comput. Methods Appl. Mech. Engrg.
147 (1997) 401-409
405
element sizes
In this investigation a simple geometric measure; h/r where h is the size of the element and r is the shortest distance from the element centre to the path of the heat source until next remeshing is performed, was used to drive the refinement and coarsening. This makes it possible to refine in advance of the motion of the heat source. The user can then choose an optimal combination of the amount of refinement /coarsening each time and how often this should be performed. Refinement or coarsening decisions are ideally driven by error measures based on the solution and the thermal, stress and microstructural fields present. In the future, generic error measures such as [29] which can evaluate the error in the temperature, stress, plastic strain or microstructure fields will be used to determine which elements to refine or coarsen. A posteriori measures of this type can only change the mesh once a trial solution is available. With a moving heat source, as is found in welding analysis, it would be a distinct advantage to alter the mesh in response to the source motion before generating a solution at each step.
4. Longitudinal
butt welding of a copper canister
4.1. Background In the Swedish nuclear waste program it has been proposed that spent nuclear fuel shall be placed in composite copper-steel canisters. These canisters will be placed in holes in tunnels located some 500 m underground in a rock storage. The canister consists of two cylinders of 5 m length, one inner cylinder made of steel and one outer cylinder made of copper. The outer diameter of the canister is 500 mm and the wall thickness is 50 mm. For storage, the steel cylinder, which contains the spent nuclear fuel, is placed inside the copper cylinder. Thereafter, a copper end is butt welded to the copper cylinder using electron beam weld:ing. The primary objective in the mechanical design of the canister is to ensure no leakage of radioactive particles to the surroundings. One of the alternative fabrication schemes for the outer copper cylinder is based on cold forming (rolling) and subsequent longitudinal seam welding of the two copper half cylinders. 4.2. Computational models The potential of the proposed algorithm is demonstrated by studying the longitudinal butt welding of two copper half cylinders previously studied in [15]. The models used do not include residual stress from the rolling process. The transient and residual temperature, stress and strain fields present in the canister during longitudinal seam welding and after cooling to ambient temperature are calculated numerically by use of the coupled thermo-mechanical FE-code SIMPLE [30]. Due to symmetry the FE-model consists of one half of a copper half cylinder. The models extend 500 mm in the axial direction including the rear end where the welding finishes. The cylinder is restrained axially at the front end, the rear lend is free. In order to have a satisfying accuracy in the simulations, the elements close to the weld-path must be around 5 mm in size, transverse to the welding direction. Two models are used: one with and the other without remeshing. The meshes were chosen in order to achieve about the same accuracy. The thermal material data for the copper alloy were taken from handbooks (see [8] for the data used). A thermo-elastic-plastic model with linear isotropic hardening is used in the mechanical analysis. The material data for the mechanical analysis were determined from uniaxial tensile tests. See [15] for the elastic and plastic material data used. The load is applied by defining a moving box along the cylinder. In each loading-step the nodes in the box and those passed by the box are subjected to a uniform input of energy. The net heat power is 8 kW and the welding speed is 5 mm/s. The loading is applied during the first 100 s of the simulation which is performed in a total of 1000 s. Both calculations are performed in 324 loading-steps. The FE-model without remeshing consists of 5783 nodes and 4220 eight-node brick elements. To reduce the number of elements in the model without reducing the accuracy a dense mesh is required in the vicinity of the moving heat source but a
406
L.-E.
Lindgren
et al. I Comput. Methods Appl. Mech. Engrg.
147 (1997) 401-409
coarser mesh is used elsewhere. See Fig. 3(b) for the FE-mesh of the non-remeshing model. The FE-model used for remeshing consists initially of 1000 eight-node brick elements (1683 nodes) but increases to maximum 2715 elements (441 nodes) at 70 s. The remeshing is performed nine times. See Fig. 3(a) for the FE-mesh at different stages of the calculation. 4.3. Computed temperatures and stresses The two FE-calculations, without and with remeshing, were performed on an IBM RS6000, model 590. The calculation without remeshing required 34 h and 20 min. of CPU-time. The CPU-time for the model with remeshing was 14 h and 5 min. The size of the stiffness and the conductivity matrix of the non-remeshing model totals 9.7 million terms. The corresponding size for the remeshing model varies from 1.6 to 5.3 million. The evolution of the calculated axial stress is shown in Fig. 3. The plastic strain fields at the final stage are shown in Fig. 4. Comparison of the residual stress-states is shown in Fig. 5. The temperature at 50 s (when the beam is half-way along the cylinder) is shown in Fig. 6. The
Fig. 3. Comparison of the axial stress field (z-stress). (a) With remeshing at 10, 50, 100 and 200 s; (b) without remeshing at 50 and 100 s.
L:,E.
Lindgren
mill
I
et al. I Comput. Methods Appl. Mech. Engrg.
0.0
nnx=o.o3e
147 (1997) 401-409
407
mh=o.o mLx=O.OU
Y
a)
W
Fig. 4. Comparison of the residual effective plastic strain. (a) With remeshing; (b) without remeshing. Fringe Iev6l6 [MPal
60 75 69 64 59 53 48
43 37 32 27 21 16
min.
1.95 MPrI
llWXE61.8MPB
mm * 1.65 f4Pe mar Ia%5 MP.¶
11 5 0 15
a)
W
Fig. 5. Comparison of the residual effective stress (von Mises). (a) With remeshing; (b) without remeshing.
maximum temperature at the nodes in the centre of the weld pool reach approximatively temperature, 1083°C. However, the used postprocessor plots element average temperatures a maximum temperature of 520°C in Fig. 6.
4.4. Comparisons
the melting resulting in
of computational models
As can be seen in Figs. 3 to 5, the simulation using remeshing gives similar results to the simulation without remeshing. The reduction in size of the models decreases the required CPU-time. It also reduces the data-storage requirement substantially, this can in some applications be the key factor determining whether a successful FE-simulation can be completed or not.
L.-E.
408
Lindgren
et al. I Comput. Methods Appl. Mech. Engrg.
Fig. 6. Temperature
distribution
147 (1997) 401-409
at 50 s.
5. Discussion and conclusions The automatic remeshing procedure is found to be effective. The required CPU-time was reduced by 60% for the studied application without any loss in accuracy. Larger gains will be expected for application with longer welds as the zone associated with large gradients will be smaller relative to the total length of the weld. The current logic for deciding the size of the elements in relation to the distance to the arc will be combined with a posteriori error measure, in order to have a versatile tool that includes a predictive capability.
Acknowledgements This work was funded by the Swedish Nuclear Fuel and Waste Management Sweden and by Lulei University of Technology, Lulea, Sweden.
Co (SKB) in Stockholm,
References [II L. Karisson,
PI
[31 [41 PI PI [71 PI
Thermal stresses in welding, in: R. Hetnarski, ed., Thermal Stresses, Ch. 5 (Elsevier, Amsterdam, 1986). J. Goldak, A. Oddy, M. Gu, W. Ma, A. Mashaie and E. Hughes, Coupling heat transfer, microstructure evolution and thermal stress analysis in weld mechanics, in: L. Karlsson, M. Johnson and L.-E. Lindgren, eds., Proc. of IUTAM Symp. Mechanical Effects of Welding (Springer-Verlag, Berlin, 1992) l-30. D. Radaj, Heat Effects of Welding (Springer-Verlag, Berlin, 1992). Y. Ueda, M. Murakawa, K. Nakacho and N.-X. Ma, Establishment of computational weld mechanics, Trans. JWRI 24 (1995) 73-86. A.S. Oddy, J.A. Goldak and J.M.J. McDill, Transformation effects in the 3D finite element analysis of welds, in: Proc. of 2nd Inter. Conf. on Trends in Welding Research, Gatlinburg, USA, 1989). R.I. Karlsson and L. Josefson, Three-dimensional finite element analysis of temperatures and stresses in single-pass butt-welded pipe, Trans. ASME, J. Press. Vessel & Tech. 112 (1990) 76-84. L.-E. Lindgren and L. Karlsson, Deformations and stresses in welding of shell structures, Int. J. Numer. Methods Engrg. 25 (1988) 635-655. E. Friedman, Thermomechanical analysis of the welding process using the finite element method, Trans. ASME, J. Pressure Vessel Tech. (1975) 206-213.
L.-E.
Lindgren
et al. I Comput. Methods Appl. Mech. Engrg.
147 (1997) 401-409
409
[9] B. Andersson and I,. Karlsson, Thermal stresses in large butt-welded plates, J. Thermal Stresses 4 (1981) 491-500. [lo] M. Jonsson, L. Karlsson and L.-E. Lindgren, Deformations and stresses in butt-welding of large plates, in: R.W. Lewis, ed., Numerical Methods in Heat Transfer, III (1985) 35-38. [ll] A. Oddy, J. Goldak and M. McDilI, Transformation plasticity and residual stresses in single-pass repair welds, Trans. ASME, J. Pressure Vessel Tech. 11 (1992) 33-38. [12] J. Goldak, A. Oddy and D. Dorling, Finite element analysis of welding on fluid filled pressurized pipelines, in: Proc. 3rd Int. Conf. on Trends in Welding Research, Gathnburg, USA, 1992, 45-50. [13] A. Oddy and J.M. E&Dill, Residual stresses in weaved repair welds, in: Proc. 1996 ASME Pressure Vessel & Piping Conf., Montreal, Canada, 1966, 147-152. [14] B.L. Josefson, L. Karlsson, L.-E. Lindgren and M. Jonsson, Thermo-mechanical FE-analysis of butt-welding of Cu-Fe canister for spent fu’el, in: K.F. Kussmaul, ed., Trans. 12th Int. Conf. on Structural Mechanics in Reactor Technology, paper NO115 in Vol. N (Stuttgart, Germany, 1993). [15] B.L. Josefson, L.-E.. Lindgren, H.-A. Haggblad and L. Karlsson, Thermo-mechanical FE-analysis of the fabrication of a Cu-Fe canister for spent nuclear fuel, in: S.V. Mbller, ed., Trans. 13th Int. Conf. on Structural Mechanics in Reactor Technology (SMiRT-13), Universidade Federal do Rio Grande do Sul, Brazil, 1995, 275-280. [16] N.N. Rykalin, Energy sources used in welding, Soudage et Techniques Connexes l(12) (1974) 471-485. [17] H. Runnemalm, L...E. Lindgren, M.O. Nisstriim and C. Lampa, Accuracy in thermal analysis of laser welding, in: G.M. Carlomagno and C.A. Brebbia, eds., Proc. Comput. Methods and Experimental Measurements VII (Computational Mechanics Publications, Southampton, 1995) 85-92. [18] D. Rosenthal, Mathematical theory of heat distribution during welding and cutting, Welding Journal, welding research supplement (1941) 220-234. [19] J.M. Bergheau, D. Pont and J.B. Leblond, Three-dimensional simulation of a laser surface treatment through steady-state computation in the heat source’s comoving frame, in: L. Karlsson, M. Johnson and L.-E. Lindgren, eds., Proc. IUTAM Symp. Mech. Effects of Welding (Springer-Verlag, Berlin, 1992) 85-92. [20] S. Brown, K. Christie and H. Song, Three dimensional finite element modeling of welded structures, in: L. Karlsson, M. Johnson and L.-E. Lindgren, eds., Proc. IUTAM Symp. Mech. Effects of Welding (Springer-Verlag, Berlin, 1992) 181-188. [21] A.O. Cifuentes and A. Kalbag, A performance study of tetrahedral and hexahedral elements in 3-D finite element structural analysis, Fin. Elem. Anal. Des. 12 (1992) 313-318. [22] S.E. Benzley, E. Perry, K. Merkely, B. Clark and G.D. Sjaardama, A comparison of all hexagonal and all tetrahedral finite element meshes for elastic and elasto-plastic analysis, in: Proc. 14th Ann. Int. Meshing Roundtable, Albuquerque, USA, 1995. [23] B.A. Szabo, Estimation and control of error based on P-convergence, in: Proc. Int. Conf. on Accuracy Estimates and Adaptive Refinements in Finite Element Computations (ARFEC), Lisbon, Portugal, 1984. [24] J.M. McDill, J.A. Goldak, A.S. Oddy and M.J. Bibby, Isoparametric quadrilaterals and hexahedrons for mesh grading algorithms, Comm. Appl. Numer. Methods 3 (1987) 155-163. [25] J.M.J. McDill, AS. Oddy and J.A. Goldak, An adaptive mesh-management algorithm for three-dimensional automatic finite element analysis, Trans. CSME 15 (1991) 57-70. [26] J.M.J. McDill and .4.S. Oddy, Arbitrary coarsening for adaptive mesh management in automatic finite element analysis, J. Math. Modelling Sci. Comput. 2B (1993) 1072-1077. [27] J.M.J. McDill, ME. Klein and A.S. Oddy, Data transfer for 3D h-adaptive thermal-elasto-plastic finite element analysis, simulations of materials processing: theory, methods and applications, in: Proc. Numiform 95, Ithaca, New York, 1995, 463-468. [28] J.M.J. McDill and A.S. Oddy, A non-conforming eight to 26-node hexahedron for three-dimensional thermal-elastoplastic finite element analysis, Comput. Struct. 54 (1995) 183-189. [29] J.M.J. McDill, M.A. Johnston and A.S. Oddy, A generic error estimate for three-dimensional h-adaptive finite element analysis, submitted to the Eleventh Int. Conf. on Mathematical Modelling and Scientific Computing, Sept. 1996. [3O] L.-E. Lindgren, User’s manual for SIMPLE-Finite Element Program for Simulation of Material Processing, Lulei University of Technology, Sweden, 1996.