- Email: [email protected]
Improved sand surface element for residual stress determination
Applied Mathematical Modelling 28 (2004) 533–546 www.elsevier.com/locate/apm
Improved sand surface element for residual stress determination Anthony Chang *, Jonathan Dantzig University of Illinois at Urbana-Champaign, 351 meb, mc 244, 1206 West Green Street, Urbana, IL 61801, USA Received 15 July 2002; received in revised form 2 September 2003; accepted 24 October 2003
Abstract A method has been developed to efficiently predict residual stresses in castings whose molds include steel chills and other high stiffness mold components. Mold elements are replaced by surface elements that apply an equivalent normal force to the surface of the part. A special contact condition, consisting of two separate functions, was developed to handle a large range of mold stiffnesses. Solidification of several geometric shapes was simulated and the solutions using the new element were compared with solutions computed using a mesh for the full mold. The surface element was found to provide a good estimate for the displacement and equivalent plastic strain for the casting models. The time necessary to compute the solution was also significantly decreased when using this new technique. Ó 2003 Elsevier Inc. All rights reserved.
1. Introduction Casting is an important industrial process used to make near net shaped products such as engine blocks, crankshafts and turbines. The determination of residual stresses in these complex parts is important for preventing problems such as cold cracking, hot tearing and warping. Test castings are commonly used to ensure that mold patterns produce parts that are free of defects. However, with the advancement of computer technology, it has become feasible to simulate this process to determine problem areas, and improve designs. Numerical simulations using the finite element method (FEM) have been successful in predicting residual stress patterns in complex parts [1–3]. With this method, it is also possible to evaluate changes to the geometry of the part without resorting to costly test castings. Some *
Corresponding author. Tel.: +1-217-244-6000. E-mail address: [email protected] (A. Chang).
0307-904X/$ - see front matter Ó 2003 Elsevier Inc. All rights reserved. doi:10.1016/j.apm.2003.10.016
534
A. Chang, J. Dantzig / Appl. Math. Modelling 28 (2004) 533–546
experience and skill is needed to produce the complex meshes that capture the geometry of intricate parts. While advanced mesh generation programs are making it easier to construct such complex meshes, solution time remains a significant problem for complicated simulations. Although the mold provides important thermal and mechanical boundary conditions, the detailed information produced from the interior is seldom of interest to the analyst. However, a large fraction of the total nodes in the mesh are required to model the mold. If the mold nodes were to be eliminated, and replaced by an equivalent boundary condition, then the time needed to obtain a solution would be significantly reduced. The mold and casting both deform under non-uniform thermal loads, creating areas of contact, where forces are transmitted, and gaps, where the mold and casting are free to move. Capturing the interactions between the mold and the part is essential to obtaining an accurate solution for residual stresses in casting simulations. Several different methods have been developed to model the contact between two solid bodies. Chaudhary and Bathe [4] used a Lagrange multiplier method to enforce compatibility of displacements for contacting surfaces. Oden and Kikuchi [5] used a reduced integration-penalty method to simulate rigid body contact with an elastic medium. Additional forces can also be directly added to the right-hand side of the stiffness equations [6]. Bellet et al. [1] used a rigid mold condition that allowed for gap formation. However, in this model the mold stiffness is altered to maintain a zero displacement condition and it is possible that the stiffness of the mold may be overestimated. Later, Bellet et al. [2] used a penalty method to restrain the penetration of the metal into the mold. Although these methods can be used to model the interaction between the mold and the part, they all require the discretization of all of the contacting bodies. Some simpler methods have also been used to model the forces applied by the mold surface. An elastic foundation element can be used to replace the mold nodes, essentially, attaching a series of springs to the surface of the part [7]. This model does not allow gap formation since the element produces mechanical resistance to motion in either direction. In some problems, it may be feasible to assume that no interactions with the mold take place and therefore, the mold nodes can simply be removed from the calculations [8]. However, this is often not a realistic assumption. More elaborate methods must be formulated to deal with the varying conditions of contact found in casting simulations. The previous work of Dantzig and coworkers [9–11] on the boundary curvature method was used as a model for the present development of a surface contact element. In the boundary curvature method, the thermal resistance of the mold is simulated by a surface element that has similar local heat transfer characteristics. This significantly reduces the computation time necessary to obtain a transient thermal solution for casting simulations. By assuming that the thermal contact resistance between the mold and the part is independent of the displacement, it is possible to partially decouple the mechanical and thermal solution. The thermal solution can then be used as an input to the mechanical problem. For this work, we assume the thermal resistance is negligible since we are primarily concerned with the modeling of the mechanical constraints provided by the mold. Taking the contact resistance to be a known function of surface temperature retains the independence of the thermal and mechanical solutions. A subsequent stress analysis can be performed using the sand surface element developed by Dantzig and coworkers [3]. This element applies an equivalent normal force to the surface of the casting when compaction occurs, and gradually releases when a gap appears. This is accomplished
A. Chang, J. Dantzig / Appl. Math. Modelling 28 (2004) 533–546
Iron
535
Mold
Ratchet
Iron
Spring
Fig. 1. Conceptual sketch of the surface contact element. The top right shows a metal part surrounded by a fully meshed mold. The bottom image shows the metal part with the surface element applied. [3]
by creating a non-linear spring and ratchet element. The method was implemented via the user element function (UEL) in ABAQUSTM . The element is illustrated in Fig. 1. The spring provides the necessary mechanical resistance to compaction and the ratchet mechanism accounts for permanent deformation of the mold surface. The work reported here expands the capabilities of the surface element to accommodate the presence of higher stiffness mold materials. It will be shown that the combination of the boundary curvature method and the sand surface element to simulate residual stresses in castings can greatly decrease the amount of time necessary to obtain solutions to these types of problems, while retaining acceptable accuracy. 2. Theory 2.1. Surface element The contact element is formulated as a 4-noded bilinear quadrilateral surface element that is to be placed on the external surface of an 8-noded trilinear volume element. The element is implemented using the user element (UEL) function in ABAQUSTM . ABAQUSTM uses NewtonÕs method to solve for the nodal displacements. It requires a trial solution vector of displacements U, oF . The essence of our method and returns the nodal force vector F and the mold tangent stiffness, oU is the precalculation of mold stiffness at each location on the surface, and we describe that now. To determine the local mold stiffness we discretize the mold without the part, as illustrated in Fig. 2. Consider the shaded region of the mold shown in Fig. 2. If a unit displacement is applied normal to the casting surface at one end, and the other end is fixed, the effective stiffness b can be computed as simply b¼
E ; L
ð1Þ
536
A. Chang, J. Dantzig / Appl. Math. Modelling 28 (2004) 533–546
Fig. 2. Unit normal displacement at each interface to determine local stiffness.
where E is the elastic modulus of the sand and L is the length of the beam. In more complex shapes, the local stiffness of each element can be found by generalizing this procedure as follows: 1. Impose a unit displacement normal to the surface of the mold, as illustrated in Fig. 2. 2. Perform a static mechanical analysis on the mold, treating it as a linear elastic continuum to obtain the stresses due to the imposed displacement. Since this is a simple static analysis, the computation time necessary to determine the stresses is relatively short compared to the transient solution to solve for the residual stresses. 3. The average stress in the direction normal to the surface of the mold is then determined and the local mold stiffness is assigned as: en ; bm ¼ r
ð2Þ
en is the average stress in the direction normal to the where bm is the local mold stiffness and r surface. The mold mesh is no longer necessary and it is removed from the analysis and replaced by the surface elements. oF ¼ be . However, thermal shrinkage, may cause the When the part is in contact with the mold, oU part to lose contact with the mold, and the stiffness must be reduced to represent this phenomenon. We discuss the method in the context of its implementation in ABAQUSTM using the UEL, but the method does not require that ABAQUSTM be used. oF proceeds in several steps. Some of the details have been described The calculation of F and oU elsewhere [12] and we give here only enough for the present development. First, the outward normal vector is determined for each element on the surface of the part. A local coordinate system can be calculated for each element based on the node numbering convention in ABAQUSTM . Then a rotation matrix is used to convert the displacement vector associated with each node from the global to the local coordinate system. The stiffness be is modified based on the direction and
A. Chang, J. Dantzig / Appl. Math. Modelling 28 (2004) 533–546
537
magnitude of the displacement vector. The local stiffness matrix is then formulated from these modified stiffnesses. Finally, the local stiffness matrix is transformed back into the global cooroF . The local displacement is multiplied by the local stiffness matrix to dinate and is returned as oU produce the element contribution to the global force vector, F. Gap formation is neglected in the thermal analysis, so the temperature history is determined before the mechanical analysis. We use FIDAPTM for this purpose, but any program could be used. The thermal solution is imported into ABAQUSTM using FIDAQUS, a conversion program developed by our lab. The thermal history is then used to determine material properties during the stress analysis. The entire process is illustrated in Fig. 3.
Thermal Solution from FIDAP
Transfer solution to ABAQUS input file
ABAQUS Input File
Unit displacement calculation
Stress from Gaussian quadrature points
Remove Sand Nodes
Calculate rotation matrix
Compute local displacement vector
Compute local stiffness matrix
Compute global stiffness matrix ABAQUS
TM
Calculate residual force
Fig. 3. Flowchart for simulation procedure.
UEL
Assign stiffness
538
A. Chang, J. Dantzig / Appl. Math. Modelling 28 (2004) 533–546
2.2. Contact models The most difficult aspect of the problem is to handle the contact/gap condition. We describe several possible approaches here to motivate our new work. The simplest contact function depends on the position of the part relative to the mold. If thermal shrinkage causes a gap to form, then no mechanical resistance is applied. However, if the part displaces into the mold, the element provides a resistance equivalent to the stiffness, bm , determined from the unit normal displacement. In functional form: oF ¼ oU
bm 0
: d < 0; : d P 0;
ð3Þ
where d is the gap thickness and bm is the local stiffness obtained from the unit normal displacement. Although this algorithm provides a realistic model of the interaction between the mold and the part, it causes numerical difficulties. The sharp transition from zero stiffness to bm can cause ‘‘contact bounce’’ where the stiffness oscillates between the two conditions. This, in turn, causes a loss of convergence and an inability to determine a solution for the residual stress. The zero stiffness condition can lead to ill-conditioning of the stiffness matrix, also causing loss of convergence. Metzger et al. [3] developed a method to model the stiffness of the mold while avoiding the numerical difficulties associated with the simple contact model. The model introduces a gradual release of contact using two parameters, c and a, oF d < 0; bm ð4Þ ¼ cbm þ ð1 cÞbm ead d P 0; oU c represents a minimum stiffness that allows the function to avoid problems with ill conditioning. The ‘‘soft contact’’ parameter, a, defines how quickly the stiffness decreases as a function of gap thickness. To avoid contact bounce and provide an adequate rate of decay for the stiffness, a value of a ¼ 10 mm1 was used. A value of c ¼ 0:1 was used to specify a minimum local stiffness of approximately 20 MPa. Since gravity is not included in these simulations, it is conceivable that a convex shape could lose contact completely from the mold. The minimum local stiffness prevents ill conditioning of the stiffness matrix when such a case arises. This function was intended to model low stiffness mold materials, and so the position of the mold surface is stored after compaction. This value is then used in subsequent calculations to evaluate d. This function works quite well in modeling the interactions between the part and soft mold materials [3]. Fig. 4 compares the low stiffness contact condition with the simple contact model. Although both conditions seem similar there is a problem. The real mold surface is able to immediately release the casting if it contracts away from the mold. In the soft contact model, the surface element gives mechanical resistance to gap formation that is equal to bm . For relatively low stiffness materials, this resistance does not provide enough force to significantly affect gap formation. However, if the mold stiffness is relatively high compared to the elastic modulus of the casting, then this artificial mechanical restraint can cause the part to stick in the high stiffness
A. Chang, J. Dantzig / Appl. Math. Modelling 28 (2004) 533–546
Simple Contact Soft Contact Method
6000
Stiffness (MPa)
539
β m 4000
2000
0 –1.0
–0.5
0.0
0.5
1.0
Gap Thickness (mm)
Fig. 4. Comparison between the simple contact and soft contact stiffness conditions.
regions. Another method must be used to model contact interactions between the casting and these high stiffness molds. 2.3. High stiffness contact model Chills made of high stiffness materials are commonly used to provide directional solidification in sand casting. These high stiffness mold materials are modeled using a new function to modify the surface element stiffness, introducing the new parameters, x, B, and C: 8 bm dþx > > 1 þ c þ ð1 cÞ tanh < oF 2 C ¼ b dþx oU > m > :Bþ 1 þ c þ ð1 cÞ tanh C 2
d > 0; ð5Þ d 6 0;
x is a constant factor that shifts the contact point to be slightly inside the mold, and C controls the width of the hyperbolic tangent function used to smoothly increase the stiffness. A small discontinuity parameter, B, is also used to increase the initial mechanical restraint to compaction. To preserve a low resistance to gap formation, the discontinuity parameter does not appear when d is positive. For steel chills, B ¼ 150 MPa was found to provide suitable resistance to penetration of the part into the mold. B was chosen to provide the maximum possible resistance to initial compaction without causing contact bounce. A minimum stiffness was defined by setting c ¼ 0:0001, which provided a minimum stiffness on the same order of magnitude as the low stiffness contact model. After some experimentation a shift factor of x ¼ 2:6 mm and width parameter C ¼ 0:8 mm were found to provide a separation force similar to the forces computed using the low stiffness contact model. This contact model is sensitive to the chosen value of B which tends to dominate the initial contact stiffness of the mold. C and x are used to control the stability of the algorithm, but they do not tend to change the results of the simulations significantly. Using these parameters, the model lowers the interface stiffness to 0.7% of the value determined from the unit normal displacement when d > 0. The behavior of the new contact function is illustrated in Fig. 5.
540
A. Chang, J. Dantzig / Appl. Math. Modelling 28 (2004) 533–546
Simple Contact Modified
6000.0
Stiffness (MPa)
βm 4000.0
ω C 2000.0
0.0 –5.0
B –3.0
–1.0
1.0
Gap Thickness (mm)
Fig. 5. Comparison between the simple contact and high stiffness contact conditions.
Using the new high stiffness mold contact function, it is possible to lower the separation forces at the surface of the part to a level that permits gap formation. The interface position is not changed in this model, since it is considered unlikely that the thermal contraction of the part will cause significant permanent deformation to the high stiffness chills. This function allows some artificial compaction to develop since substantial penetration must occur before the maximum stiffness is encountered. However, the value of B is chosen such that it is significantly larger than the stiffness of sand. This generally prevents compaction at any point on the surface of the steel mold. A combination of the high stiffness and low stiffness contact conditions can be used to span a large range of mold materials. In this implementation, if the mold stiffness is larger than 200 MPa, then the high stiffness contact method is used. When the stiffness of the mold is less than 200 MPa, the low stiffness contact method is implemented. The method is illustrated by example problems in Section 3.
3. Results 3.1. L-shaped bar To validate the algorithm used to model high stiffness mold materials, we constructed a simple two-dimensional, plane strain, L-shaped model, illustrated in Fig. 6. A plane of symmetry is imposed on the right-hand side. The casting is made of grey iron, and mold is a silica-based green sand. The thermal and mechanical properties are tabulated by Metzger et al. [3]. The steel chill cools the end blocks, which constrain the thin center section from contracting. This leads to the development of large plastic strains in this region of the casting. A good approximation of the plastic strain can be determined using ABAQUSTM gap elements. These elements are used to model gap formation while the mold elements are used to provide the appropriate mechanical restraint to compaction. Although this method models the surface conditions at the interface of the part more exactly, it requires the sand mold to be included in the analysis. Gap elements are generally non-linear and require large numbers of iterations to obtain convergence, making them
A. Chang, J. Dantzig / Appl. Math. Modelling 28 (2004) 533–546
541
Fixed
Steel Chill Sand Mold 40 mm Symmetry
90 mm (a)
30 mm
50 mm
160 mm
Symmetry
Fixed
180 mm
End block
CL
Sand Mold Fixed (b)
Fig. 6. 2-D L-shaped model and the displacements for the full mold solution: (a) L-shaped 2-D plane strain model, (b) deformed mesh after 6000 s (displacements magnified by a factor of 10).
computationally expensive. Large rotations and displacements may cause problems if gap orientation changes significantly during the analysis. However, the elements provide an excellent means of calculation contact forces. The thermal solution for the full mold model was determined using FIDAPTM . The left side of the chill was set to a constant temperature of 25 °C and the initial temperature of the part was specified as 1350 °C. Perfect contact between the chill, sand mold and part was assumed in the thermal solution. The thermal history was then sent to ABAQUSTM to perform the stress analysis. Symmetry conditions were applied on the right face of the model while the other outer mold surfaces were held fixed. Symmetry conditions were also applied to the front and back of the model, simulating plane strain. The full mold solution contained 3720 nodes while the surface element model required only 1052 nodes. In both of these simulations, the mold surface friction and gravity are neglected. There is no fundamental reason to prevent the extension of this element to include the effects of these phenomena. Both the thermal and stress analysis were simulated to 6000 s. The maximum time step used in the thermal solution was 50 s, and the largest time step taken in the stress solution was 200 s. The computed deformed shape after 6000 s is shown in Fig. 6. The centerline displacements of the casting computed using the two methods are compared in Fig. 7. The penetration of the end block into the steel chill is overestimated by about 0.1 mm when using the UEL, due to the rampup part of the high stiffness contact algorithm underestimating the resistance to deformation. However, the general characteristics of the deformed shape are captured well. To obtain a more quantitative measure of the difference between the solutions using the two methods, an error analysis was done using the equivalent plastic strain at the quadrature points of each element. The RMS error was calculated as sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi P 2 ðDÞ ; ð6Þ RMS % ¼ 100 n where D is the difference in plastic strain between the two solutions and n is the number of elements. The maximum percent difference was also computed:
542
A. Chang, J. Dantzig / Appl. Math. Modelling 28 (2004) 533–546
Displacement (mm)
1.0
0.5
Full Mold UEL
0.0
0.5
0.0
50.0
100.0
Position (mm)
150.0
200.0
Fig. 7. Vertical displacements along the centerline of the L-shaped bar for full mold and user element models.
Max Disp % ¼ 100
Dmax ;
ð7Þ
where Dmax is the maximum difference in equivalent plastic strain at a point and is the plastic strain in the full mold solution at the same position and time as Dmax . Since the majority of the high plastic strains occur in the test section, only the nodes in this area are used in the error analysis. Otherwise, the RMS error would be artificially low. The RMS error for this model was determined to be 0.21% and the maximum percent error was 24%. Fig. 8 shows the contour plots of equivalent plastic strain obtained from the surface element and full mold methods. Fig. 9 illustrates the absolute difference in equivalent plastic strain between the two solution methods. The maximum difference occurs in the high strain test region. Some large errors are also calculated in the compressive region of the beam and several of the corners. However, the strain differences in the corners are generally isolated to a single element. These errors probably stem from the difficulty in calculating an unambiguous representation of the mold stiffness in corners. The time necessary to complete these simulations is significantly reduced by implementing the new surface element. To obtain the stress solution, the full mold method required 7298 s for
Fig. 8. Equivalent plastic strain for L-shaped model at 6000 s obtained using the full mold and UEL methods: (a) equivalent plastic strain in the full mold solution, (b) equivalent plastic strain using the UEL.
A. Chang, J. Dantzig / Appl. Math. Modelling 28 (2004) 533–546
543
PEEQ Diff 0.00587705 0.00545835 0.00503965 0.00462095 0.00420225 0.00378355 0.00336485 0.00294615 0.00252745 0.00210875 0.00169005 0.00127135 0.00085265 0.00043395
Fig. 9. Plastic strain difference in the L-shaped bar for the two solution methods after 6000 s.
112 Newton iterations while the UEL method completed its solution in 1368 s and 177 Newton iterations. These times were obtained using a SunTM Ultra 10 Workstation running at 333 MHz. This reduction in time is directly related to removal of a large number of mold elements. 3.2. Shear bar Fig. 10 illustrates another 2-D plane strain model used for comparison between the full mold and surface element methods. This geometry introduces significant bending and rotation. The end blocks tend to restrict contraction of the part, causing shear strains to develop in the test section. The shear bar has symmetry conditions on the front and back faces to simulate plane strain. The full mold analysis used 7102 nodes, while the corresponding surface element model required only 1058 nodes. Each analysis simulated 6000 s of time. The maximum time step for the thermal analysis was again set to 50 s while the stress analysis had time steps no greater than 100 s. The thermal solution for this model was obtained using the boundary curvature method. This method also assumes that perfect contact is maintained between the chills, sand mold, and cast iron part. The backs of the steel chills were set to a constant 25 °C and the sand mold was also assumed to initially be at 25 °C. The metal part had a constant initial temperature of 1350 °C. This thermal solution was used for both analysis methods.
100 mm End block Sand Mold
260 mm
50 mm
270 mm
130 mm
50 mm
300 mm Sand Mold 60 mm Steel Chill Blocks
Fig. 10. Shearing 2-D, plane strain model used in verification.
544
A. Chang, J. Dantzig / Appl. Math. Modelling 28 (2004) 533–546
Fig. 11. 2-D, plane strain mesh configuration and deformed full mold solution: (a) mesh configuration for shear model, (b) deformed mesh after 6000 s using gap elements. Displacement have been magnified by a factor of 10.
The deformed shape of the shear bar is shown in Fig. 11. The general shape of the deformation is similar for both solutions and the magnitudes compare reasonably well. The deformation along the centerline of the upper left section for the UEL and full mold analysis is shown in Fig. 12. Some discrepancies arise because the end block compaction of the steel chill is again overestimated, due to the lack of initial mechanical restraint provided by the high stiffness contact algorithm. Similar error analysis to the previous case compared the values of the plastic strains computed using the two methods. Again, the largest values of the plastic strain appear in the test section and therefore, only the values from this region are used in the error analysis. The RMS error was determined to be 0.28% for this model and a maximum error of 20% was obtained. In both the Lshaped bar and the shear bar, the maximum difference in plastic strain is a bit high. However, the bending and rotation of this part is a worst case scenario for the surface element. Even under these conditions, the RMS error is reasonable. The equivalent plastic strain contour plots are shown in Fig. 13. Fig. 14 shows the absolute difference in equivalent plastic strain for the two solution methods. Although both solutions appear symmetric in Fig. 13, the fact that we use an absolute difference of the two solutions produces the slight loss of symmetry in Fig. 14. The maximum difference is
Displacement (mm)
1.0
0.5
0.0
–0.5 Full Mold UEL
–1.0 0.0
100.0
200.0
300.0
Position (mm)
Fig. 12. Vertical displacement along the centerline of the top shear beam member.
A. Chang, J. Dantzig / Appl. Math. Modelling 28 (2004) 533–546
PEEQ
545
VALUE +8.53E-03 +1.45E-02 +1.95E-02 +2.44E-02 +2.94E-02 +3.43E-02 +3.93E-02 +4.42E-02 +4.92E-02 +5.41E-02 +5.91E-02 +6.40E-02 +6.90E-02 +7.44E-02
(a)
(b)
Fig. 13. Equivalent plastic strain for shear bar model at 6000 s obtained using the full mold and UEL methods: (a) equivalent plastic strain in the full mold solution, (b) equivalent plastic strain using UEL.
PEEQ Diff 0.0106411 0.00988194 0.00912281 0.00836368 0.00760455 0.00684542 0.00608629 0.00532716 0.00456803 0.0038089 0.00304977 0.00229064 0.00153151 0.00077238
Fig. 14. Plastic strain difference in the shear bar for the two solution methods after 6000 s.
again found in the areas of the highest strain. The compressive side of the bending shear bars are also associated with some fairly high errors. This seems to indicate that the model has some difficulties in predicting bending and the interaction of corners. However, these errors are again relegated to just one or two elements in the high strain regions. The time necessary to complete the stress analysis using the full mold method was 18163 s. The UEL method required only 1427 s to obtain a solution. Although some accuracy is compromised, the time necessary to complete this analysis is significantly reduced with the use of the surface element. Since a larger fraction of the total nodes in the shear bar model are located in the mold, a greater reduction in time is obtained, when compared with the L-shaped bar analysis.
546
A. Chang, J. Dantzig / Appl. Math. Modelling 28 (2004) 533–546
4. Conclusions Simulations to predict residual stresses in casting are usually performed including the mold in the analysis. An overwhelming amount of computation time can be required to obtain both the thermal and mechanical solution in the mold. However, the mold elements are really useful only for providing accurate boundary conditions on the surface of the cast part. The development of a surface element by Metzger et al. [3] provided a method to eliminate the mold nodes from the stress analysis and replace them with surface elements that provide an equivalent normal force to the part. The contact algorithm used in that work failed for higher stiffness materials, such as those commonly used for chills. This problem was alleviated by adding a new contact condition that adequately describes higher stiffness materials. The new contact algorithm was implemented using the user element function in ABAQUSTM . Two cases were examined to test the ability of the user element to accurately predict the deformation and plastic strain determined by using a full mold solution. The general deformed shape can be predicted quite well and the overall plastic strain can be calculated accurately. Local errors at a few points can be large. This is accomplished while reducing the time necessary to obtain these solutions significantly. It appears that the surface element approach has some difficulties predicting bending and rotation of the cast part. Some penetration is also overestimated due to the significant amount of compaction necessary for the initiation of the maximum stiffness in the high stiffness contact model.
References [1] M. Bellet, F. Decultieux, M. Menai, F. Bay, C. Levaillant, J. Chenot, P. Schmidt, I. Svensson, Metall. Trans. B 27B (1996) 81. [2] M. Bellet, C. Aliaga, O. Jaouen, in: P.R. Sahm, P.N. Hansen, J.G. Conley (Eds.), Modeling of Casting, Welding and Advanced Solidification Processes––IX, Shaker verlag, Aachen, 2001. [3] D. Metzger, K.J. New, J. Dantzig, Appl. Math. Model. 25 (2001) 825. [4] A.B. Chaudhary, K. Bathe, Comp. Struct. 24 (6) (1986) 855. [5] J.T. Oden, N. Kikuchi, Int. J. Numer. Meth. Eng. 18 (1982) 701. [6] S.K. Arya, G.A. Hegemier, J. Struct. Div. ASCE 108 (1982) 327. [7] D.L. Metzger, Development of a Sand Surface Element for Improved Residual Stress Determination, MasterÕs Thesis, University of Illinois at Urbana-Champaign, 1999. [8] A. Jacot, S. Cockcroft, Metall. Mater. Trans. A 31A (2000) 1201. [9] J.A. Dantzig, S.C. Lu, Metall. Trans. B 16B (2) (1985) 195. [10] J.A. Dantzig, J.W. Wiese, Metall. Trans. B 16B (2) (1985) 203. [11] J.W. Wiese, J.A. Dantzig, Appl. Math. Model. 12 (1988) 213. [12] K.J. New, An Improved Method for Simulating the Contact Resistance of Sand Molds on Gray Iron Castings, MasterÕs Thesis, University of Illinois at Urbana-Champaign, 1997.
Improved sand surface element for residual stress determination Anthony Chang *, Jonathan Dantzig University of Illinois at Urbana-Champaign, 351 meb, mc 244, 1206 West Green Street, Urbana, IL 61801, USA Received 15 July 2002; received in revised form 2 September 2003; accepted 24 October 2003
Abstract A method has been developed to efficiently predict residual stresses in castings whose molds include steel chills and other high stiffness mold components. Mold elements are replaced by surface elements that apply an equivalent normal force to the surface of the part. A special contact condition, consisting of two separate functions, was developed to handle a large range of mold stiffnesses. Solidification of several geometric shapes was simulated and the solutions using the new element were compared with solutions computed using a mesh for the full mold. The surface element was found to provide a good estimate for the displacement and equivalent plastic strain for the casting models. The time necessary to compute the solution was also significantly decreased when using this new technique. Ó 2003 Elsevier Inc. All rights reserved.
1. Introduction Casting is an important industrial process used to make near net shaped products such as engine blocks, crankshafts and turbines. The determination of residual stresses in these complex parts is important for preventing problems such as cold cracking, hot tearing and warping. Test castings are commonly used to ensure that mold patterns produce parts that are free of defects. However, with the advancement of computer technology, it has become feasible to simulate this process to determine problem areas, and improve designs. Numerical simulations using the finite element method (FEM) have been successful in predicting residual stress patterns in complex parts [1–3]. With this method, it is also possible to evaluate changes to the geometry of the part without resorting to costly test castings. Some *
Corresponding author. Tel.: +1-217-244-6000. E-mail address: [email protected] (A. Chang).
0307-904X/$ - see front matter Ó 2003 Elsevier Inc. All rights reserved. doi:10.1016/j.apm.2003.10.016
534
A. Chang, J. Dantzig / Appl. Math. Modelling 28 (2004) 533–546
experience and skill is needed to produce the complex meshes that capture the geometry of intricate parts. While advanced mesh generation programs are making it easier to construct such complex meshes, solution time remains a significant problem for complicated simulations. Although the mold provides important thermal and mechanical boundary conditions, the detailed information produced from the interior is seldom of interest to the analyst. However, a large fraction of the total nodes in the mesh are required to model the mold. If the mold nodes were to be eliminated, and replaced by an equivalent boundary condition, then the time needed to obtain a solution would be significantly reduced. The mold and casting both deform under non-uniform thermal loads, creating areas of contact, where forces are transmitted, and gaps, where the mold and casting are free to move. Capturing the interactions between the mold and the part is essential to obtaining an accurate solution for residual stresses in casting simulations. Several different methods have been developed to model the contact between two solid bodies. Chaudhary and Bathe [4] used a Lagrange multiplier method to enforce compatibility of displacements for contacting surfaces. Oden and Kikuchi [5] used a reduced integration-penalty method to simulate rigid body contact with an elastic medium. Additional forces can also be directly added to the right-hand side of the stiffness equations [6]. Bellet et al. [1] used a rigid mold condition that allowed for gap formation. However, in this model the mold stiffness is altered to maintain a zero displacement condition and it is possible that the stiffness of the mold may be overestimated. Later, Bellet et al. [2] used a penalty method to restrain the penetration of the metal into the mold. Although these methods can be used to model the interaction between the mold and the part, they all require the discretization of all of the contacting bodies. Some simpler methods have also been used to model the forces applied by the mold surface. An elastic foundation element can be used to replace the mold nodes, essentially, attaching a series of springs to the surface of the part [7]. This model does not allow gap formation since the element produces mechanical resistance to motion in either direction. In some problems, it may be feasible to assume that no interactions with the mold take place and therefore, the mold nodes can simply be removed from the calculations [8]. However, this is often not a realistic assumption. More elaborate methods must be formulated to deal with the varying conditions of contact found in casting simulations. The previous work of Dantzig and coworkers [9–11] on the boundary curvature method was used as a model for the present development of a surface contact element. In the boundary curvature method, the thermal resistance of the mold is simulated by a surface element that has similar local heat transfer characteristics. This significantly reduces the computation time necessary to obtain a transient thermal solution for casting simulations. By assuming that the thermal contact resistance between the mold and the part is independent of the displacement, it is possible to partially decouple the mechanical and thermal solution. The thermal solution can then be used as an input to the mechanical problem. For this work, we assume the thermal resistance is negligible since we are primarily concerned with the modeling of the mechanical constraints provided by the mold. Taking the contact resistance to be a known function of surface temperature retains the independence of the thermal and mechanical solutions. A subsequent stress analysis can be performed using the sand surface element developed by Dantzig and coworkers [3]. This element applies an equivalent normal force to the surface of the casting when compaction occurs, and gradually releases when a gap appears. This is accomplished
A. Chang, J. Dantzig / Appl. Math. Modelling 28 (2004) 533–546
Iron
535
Mold
Ratchet
Iron
Spring
Fig. 1. Conceptual sketch of the surface contact element. The top right shows a metal part surrounded by a fully meshed mold. The bottom image shows the metal part with the surface element applied. [3]
by creating a non-linear spring and ratchet element. The method was implemented via the user element function (UEL) in ABAQUSTM . The element is illustrated in Fig. 1. The spring provides the necessary mechanical resistance to compaction and the ratchet mechanism accounts for permanent deformation of the mold surface. The work reported here expands the capabilities of the surface element to accommodate the presence of higher stiffness mold materials. It will be shown that the combination of the boundary curvature method and the sand surface element to simulate residual stresses in castings can greatly decrease the amount of time necessary to obtain solutions to these types of problems, while retaining acceptable accuracy. 2. Theory 2.1. Surface element The contact element is formulated as a 4-noded bilinear quadrilateral surface element that is to be placed on the external surface of an 8-noded trilinear volume element. The element is implemented using the user element (UEL) function in ABAQUSTM . ABAQUSTM uses NewtonÕs method to solve for the nodal displacements. It requires a trial solution vector of displacements U, oF . The essence of our method and returns the nodal force vector F and the mold tangent stiffness, oU is the precalculation of mold stiffness at each location on the surface, and we describe that now. To determine the local mold stiffness we discretize the mold without the part, as illustrated in Fig. 2. Consider the shaded region of the mold shown in Fig. 2. If a unit displacement is applied normal to the casting surface at one end, and the other end is fixed, the effective stiffness b can be computed as simply b¼
E ; L
ð1Þ
536
A. Chang, J. Dantzig / Appl. Math. Modelling 28 (2004) 533–546
Fig. 2. Unit normal displacement at each interface to determine local stiffness.
where E is the elastic modulus of the sand and L is the length of the beam. In more complex shapes, the local stiffness of each element can be found by generalizing this procedure as follows: 1. Impose a unit displacement normal to the surface of the mold, as illustrated in Fig. 2. 2. Perform a static mechanical analysis on the mold, treating it as a linear elastic continuum to obtain the stresses due to the imposed displacement. Since this is a simple static analysis, the computation time necessary to determine the stresses is relatively short compared to the transient solution to solve for the residual stresses. 3. The average stress in the direction normal to the surface of the mold is then determined and the local mold stiffness is assigned as: en ; bm ¼ r
ð2Þ
en is the average stress in the direction normal to the where bm is the local mold stiffness and r surface. The mold mesh is no longer necessary and it is removed from the analysis and replaced by the surface elements. oF ¼ be . However, thermal shrinkage, may cause the When the part is in contact with the mold, oU part to lose contact with the mold, and the stiffness must be reduced to represent this phenomenon. We discuss the method in the context of its implementation in ABAQUSTM using the UEL, but the method does not require that ABAQUSTM be used. oF proceeds in several steps. Some of the details have been described The calculation of F and oU elsewhere [12] and we give here only enough for the present development. First, the outward normal vector is determined for each element on the surface of the part. A local coordinate system can be calculated for each element based on the node numbering convention in ABAQUSTM . Then a rotation matrix is used to convert the displacement vector associated with each node from the global to the local coordinate system. The stiffness be is modified based on the direction and
A. Chang, J. Dantzig / Appl. Math. Modelling 28 (2004) 533–546
537
magnitude of the displacement vector. The local stiffness matrix is then formulated from these modified stiffnesses. Finally, the local stiffness matrix is transformed back into the global cooroF . The local displacement is multiplied by the local stiffness matrix to dinate and is returned as oU produce the element contribution to the global force vector, F. Gap formation is neglected in the thermal analysis, so the temperature history is determined before the mechanical analysis. We use FIDAPTM for this purpose, but any program could be used. The thermal solution is imported into ABAQUSTM using FIDAQUS, a conversion program developed by our lab. The thermal history is then used to determine material properties during the stress analysis. The entire process is illustrated in Fig. 3.
Thermal Solution from FIDAP
Transfer solution to ABAQUS input file
ABAQUS Input File
Unit displacement calculation
Stress from Gaussian quadrature points
Remove Sand Nodes
Calculate rotation matrix
Compute local displacement vector
Compute local stiffness matrix
Compute global stiffness matrix ABAQUS
TM
Calculate residual force
Fig. 3. Flowchart for simulation procedure.
UEL
Assign stiffness
538
A. Chang, J. Dantzig / Appl. Math. Modelling 28 (2004) 533–546
2.2. Contact models The most difficult aspect of the problem is to handle the contact/gap condition. We describe several possible approaches here to motivate our new work. The simplest contact function depends on the position of the part relative to the mold. If thermal shrinkage causes a gap to form, then no mechanical resistance is applied. However, if the part displaces into the mold, the element provides a resistance equivalent to the stiffness, bm , determined from the unit normal displacement. In functional form: oF ¼ oU
bm 0
: d < 0; : d P 0;
ð3Þ
where d is the gap thickness and bm is the local stiffness obtained from the unit normal displacement. Although this algorithm provides a realistic model of the interaction between the mold and the part, it causes numerical difficulties. The sharp transition from zero stiffness to bm can cause ‘‘contact bounce’’ where the stiffness oscillates between the two conditions. This, in turn, causes a loss of convergence and an inability to determine a solution for the residual stress. The zero stiffness condition can lead to ill-conditioning of the stiffness matrix, also causing loss of convergence. Metzger et al. [3] developed a method to model the stiffness of the mold while avoiding the numerical difficulties associated with the simple contact model. The model introduces a gradual release of contact using two parameters, c and a, oF d < 0; bm ð4Þ ¼ cbm þ ð1 cÞbm ead d P 0; oU c represents a minimum stiffness that allows the function to avoid problems with ill conditioning. The ‘‘soft contact’’ parameter, a, defines how quickly the stiffness decreases as a function of gap thickness. To avoid contact bounce and provide an adequate rate of decay for the stiffness, a value of a ¼ 10 mm1 was used. A value of c ¼ 0:1 was used to specify a minimum local stiffness of approximately 20 MPa. Since gravity is not included in these simulations, it is conceivable that a convex shape could lose contact completely from the mold. The minimum local stiffness prevents ill conditioning of the stiffness matrix when such a case arises. This function was intended to model low stiffness mold materials, and so the position of the mold surface is stored after compaction. This value is then used in subsequent calculations to evaluate d. This function works quite well in modeling the interactions between the part and soft mold materials [3]. Fig. 4 compares the low stiffness contact condition with the simple contact model. Although both conditions seem similar there is a problem. The real mold surface is able to immediately release the casting if it contracts away from the mold. In the soft contact model, the surface element gives mechanical resistance to gap formation that is equal to bm . For relatively low stiffness materials, this resistance does not provide enough force to significantly affect gap formation. However, if the mold stiffness is relatively high compared to the elastic modulus of the casting, then this artificial mechanical restraint can cause the part to stick in the high stiffness
A. Chang, J. Dantzig / Appl. Math. Modelling 28 (2004) 533–546
Simple Contact Soft Contact Method
6000
Stiffness (MPa)
539
β m 4000
2000
0 –1.0
–0.5
0.0
0.5
1.0
Gap Thickness (mm)
Fig. 4. Comparison between the simple contact and soft contact stiffness conditions.
regions. Another method must be used to model contact interactions between the casting and these high stiffness molds. 2.3. High stiffness contact model Chills made of high stiffness materials are commonly used to provide directional solidification in sand casting. These high stiffness mold materials are modeled using a new function to modify the surface element stiffness, introducing the new parameters, x, B, and C: 8 bm dþx > > 1 þ c þ ð1 cÞ tanh < oF 2 C ¼ b dþx oU > m > :Bþ 1 þ c þ ð1 cÞ tanh C 2
d > 0; ð5Þ d 6 0;
x is a constant factor that shifts the contact point to be slightly inside the mold, and C controls the width of the hyperbolic tangent function used to smoothly increase the stiffness. A small discontinuity parameter, B, is also used to increase the initial mechanical restraint to compaction. To preserve a low resistance to gap formation, the discontinuity parameter does not appear when d is positive. For steel chills, B ¼ 150 MPa was found to provide suitable resistance to penetration of the part into the mold. B was chosen to provide the maximum possible resistance to initial compaction without causing contact bounce. A minimum stiffness was defined by setting c ¼ 0:0001, which provided a minimum stiffness on the same order of magnitude as the low stiffness contact model. After some experimentation a shift factor of x ¼ 2:6 mm and width parameter C ¼ 0:8 mm were found to provide a separation force similar to the forces computed using the low stiffness contact model. This contact model is sensitive to the chosen value of B which tends to dominate the initial contact stiffness of the mold. C and x are used to control the stability of the algorithm, but they do not tend to change the results of the simulations significantly. Using these parameters, the model lowers the interface stiffness to 0.7% of the value determined from the unit normal displacement when d > 0. The behavior of the new contact function is illustrated in Fig. 5.
540
A. Chang, J. Dantzig / Appl. Math. Modelling 28 (2004) 533–546
Simple Contact Modified
6000.0
Stiffness (MPa)
βm 4000.0
ω C 2000.0
0.0 –5.0
B –3.0
–1.0
1.0
Gap Thickness (mm)
Fig. 5. Comparison between the simple contact and high stiffness contact conditions.
Using the new high stiffness mold contact function, it is possible to lower the separation forces at the surface of the part to a level that permits gap formation. The interface position is not changed in this model, since it is considered unlikely that the thermal contraction of the part will cause significant permanent deformation to the high stiffness chills. This function allows some artificial compaction to develop since substantial penetration must occur before the maximum stiffness is encountered. However, the value of B is chosen such that it is significantly larger than the stiffness of sand. This generally prevents compaction at any point on the surface of the steel mold. A combination of the high stiffness and low stiffness contact conditions can be used to span a large range of mold materials. In this implementation, if the mold stiffness is larger than 200 MPa, then the high stiffness contact method is used. When the stiffness of the mold is less than 200 MPa, the low stiffness contact method is implemented. The method is illustrated by example problems in Section 3.
3. Results 3.1. L-shaped bar To validate the algorithm used to model high stiffness mold materials, we constructed a simple two-dimensional, plane strain, L-shaped model, illustrated in Fig. 6. A plane of symmetry is imposed on the right-hand side. The casting is made of grey iron, and mold is a silica-based green sand. The thermal and mechanical properties are tabulated by Metzger et al. [3]. The steel chill cools the end blocks, which constrain the thin center section from contracting. This leads to the development of large plastic strains in this region of the casting. A good approximation of the plastic strain can be determined using ABAQUSTM gap elements. These elements are used to model gap formation while the mold elements are used to provide the appropriate mechanical restraint to compaction. Although this method models the surface conditions at the interface of the part more exactly, it requires the sand mold to be included in the analysis. Gap elements are generally non-linear and require large numbers of iterations to obtain convergence, making them
A. Chang, J. Dantzig / Appl. Math. Modelling 28 (2004) 533–546
541
Fixed
Steel Chill Sand Mold 40 mm Symmetry
90 mm (a)
30 mm
50 mm
160 mm
Symmetry
Fixed
180 mm
End block
CL
Sand Mold Fixed (b)
Fig. 6. 2-D L-shaped model and the displacements for the full mold solution: (a) L-shaped 2-D plane strain model, (b) deformed mesh after 6000 s (displacements magnified by a factor of 10).
computationally expensive. Large rotations and displacements may cause problems if gap orientation changes significantly during the analysis. However, the elements provide an excellent means of calculation contact forces. The thermal solution for the full mold model was determined using FIDAPTM . The left side of the chill was set to a constant temperature of 25 °C and the initial temperature of the part was specified as 1350 °C. Perfect contact between the chill, sand mold and part was assumed in the thermal solution. The thermal history was then sent to ABAQUSTM to perform the stress analysis. Symmetry conditions were applied on the right face of the model while the other outer mold surfaces were held fixed. Symmetry conditions were also applied to the front and back of the model, simulating plane strain. The full mold solution contained 3720 nodes while the surface element model required only 1052 nodes. In both of these simulations, the mold surface friction and gravity are neglected. There is no fundamental reason to prevent the extension of this element to include the effects of these phenomena. Both the thermal and stress analysis were simulated to 6000 s. The maximum time step used in the thermal solution was 50 s, and the largest time step taken in the stress solution was 200 s. The computed deformed shape after 6000 s is shown in Fig. 6. The centerline displacements of the casting computed using the two methods are compared in Fig. 7. The penetration of the end block into the steel chill is overestimated by about 0.1 mm when using the UEL, due to the rampup part of the high stiffness contact algorithm underestimating the resistance to deformation. However, the general characteristics of the deformed shape are captured well. To obtain a more quantitative measure of the difference between the solutions using the two methods, an error analysis was done using the equivalent plastic strain at the quadrature points of each element. The RMS error was calculated as sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi P 2 ðDÞ ; ð6Þ RMS % ¼ 100 n where D is the difference in plastic strain between the two solutions and n is the number of elements. The maximum percent difference was also computed:
542
A. Chang, J. Dantzig / Appl. Math. Modelling 28 (2004) 533–546
Displacement (mm)
1.0
0.5
Full Mold UEL
0.0
0.5
0.0
50.0
100.0
Position (mm)
150.0
200.0
Fig. 7. Vertical displacements along the centerline of the L-shaped bar for full mold and user element models.
Max Disp % ¼ 100
Dmax ;
ð7Þ
where Dmax is the maximum difference in equivalent plastic strain at a point and is the plastic strain in the full mold solution at the same position and time as Dmax . Since the majority of the high plastic strains occur in the test section, only the nodes in this area are used in the error analysis. Otherwise, the RMS error would be artificially low. The RMS error for this model was determined to be 0.21% and the maximum percent error was 24%. Fig. 8 shows the contour plots of equivalent plastic strain obtained from the surface element and full mold methods. Fig. 9 illustrates the absolute difference in equivalent plastic strain between the two solution methods. The maximum difference occurs in the high strain test region. Some large errors are also calculated in the compressive region of the beam and several of the corners. However, the strain differences in the corners are generally isolated to a single element. These errors probably stem from the difficulty in calculating an unambiguous representation of the mold stiffness in corners. The time necessary to complete these simulations is significantly reduced by implementing the new surface element. To obtain the stress solution, the full mold method required 7298 s for
Fig. 8. Equivalent plastic strain for L-shaped model at 6000 s obtained using the full mold and UEL methods: (a) equivalent plastic strain in the full mold solution, (b) equivalent plastic strain using the UEL.
A. Chang, J. Dantzig / Appl. Math. Modelling 28 (2004) 533–546
543
PEEQ Diff 0.00587705 0.00545835 0.00503965 0.00462095 0.00420225 0.00378355 0.00336485 0.00294615 0.00252745 0.00210875 0.00169005 0.00127135 0.00085265 0.00043395
Fig. 9. Plastic strain difference in the L-shaped bar for the two solution methods after 6000 s.
112 Newton iterations while the UEL method completed its solution in 1368 s and 177 Newton iterations. These times were obtained using a SunTM Ultra 10 Workstation running at 333 MHz. This reduction in time is directly related to removal of a large number of mold elements. 3.2. Shear bar Fig. 10 illustrates another 2-D plane strain model used for comparison between the full mold and surface element methods. This geometry introduces significant bending and rotation. The end blocks tend to restrict contraction of the part, causing shear strains to develop in the test section. The shear bar has symmetry conditions on the front and back faces to simulate plane strain. The full mold analysis used 7102 nodes, while the corresponding surface element model required only 1058 nodes. Each analysis simulated 6000 s of time. The maximum time step for the thermal analysis was again set to 50 s while the stress analysis had time steps no greater than 100 s. The thermal solution for this model was obtained using the boundary curvature method. This method also assumes that perfect contact is maintained between the chills, sand mold, and cast iron part. The backs of the steel chills were set to a constant 25 °C and the sand mold was also assumed to initially be at 25 °C. The metal part had a constant initial temperature of 1350 °C. This thermal solution was used for both analysis methods.
100 mm End block Sand Mold
260 mm
50 mm
270 mm
130 mm
50 mm
300 mm Sand Mold 60 mm Steel Chill Blocks
Fig. 10. Shearing 2-D, plane strain model used in verification.
544
A. Chang, J. Dantzig / Appl. Math. Modelling 28 (2004) 533–546
Fig. 11. 2-D, plane strain mesh configuration and deformed full mold solution: (a) mesh configuration for shear model, (b) deformed mesh after 6000 s using gap elements. Displacement have been magnified by a factor of 10.
The deformed shape of the shear bar is shown in Fig. 11. The general shape of the deformation is similar for both solutions and the magnitudes compare reasonably well. The deformation along the centerline of the upper left section for the UEL and full mold analysis is shown in Fig. 12. Some discrepancies arise because the end block compaction of the steel chill is again overestimated, due to the lack of initial mechanical restraint provided by the high stiffness contact algorithm. Similar error analysis to the previous case compared the values of the plastic strains computed using the two methods. Again, the largest values of the plastic strain appear in the test section and therefore, only the values from this region are used in the error analysis. The RMS error was determined to be 0.28% for this model and a maximum error of 20% was obtained. In both the Lshaped bar and the shear bar, the maximum difference in plastic strain is a bit high. However, the bending and rotation of this part is a worst case scenario for the surface element. Even under these conditions, the RMS error is reasonable. The equivalent plastic strain contour plots are shown in Fig. 13. Fig. 14 shows the absolute difference in equivalent plastic strain for the two solution methods. Although both solutions appear symmetric in Fig. 13, the fact that we use an absolute difference of the two solutions produces the slight loss of symmetry in Fig. 14. The maximum difference is
Displacement (mm)
1.0
0.5
0.0
–0.5 Full Mold UEL
–1.0 0.0
100.0
200.0
300.0
Position (mm)
Fig. 12. Vertical displacement along the centerline of the top shear beam member.
A. Chang, J. Dantzig / Appl. Math. Modelling 28 (2004) 533–546
PEEQ
545
VALUE +8.53E-03 +1.45E-02 +1.95E-02 +2.44E-02 +2.94E-02 +3.43E-02 +3.93E-02 +4.42E-02 +4.92E-02 +5.41E-02 +5.91E-02 +6.40E-02 +6.90E-02 +7.44E-02
(a)
(b)
Fig. 13. Equivalent plastic strain for shear bar model at 6000 s obtained using the full mold and UEL methods: (a) equivalent plastic strain in the full mold solution, (b) equivalent plastic strain using UEL.
PEEQ Diff 0.0106411 0.00988194 0.00912281 0.00836368 0.00760455 0.00684542 0.00608629 0.00532716 0.00456803 0.0038089 0.00304977 0.00229064 0.00153151 0.00077238
Fig. 14. Plastic strain difference in the shear bar for the two solution methods after 6000 s.
again found in the areas of the highest strain. The compressive side of the bending shear bars are also associated with some fairly high errors. This seems to indicate that the model has some difficulties in predicting bending and the interaction of corners. However, these errors are again relegated to just one or two elements in the high strain regions. The time necessary to complete the stress analysis using the full mold method was 18163 s. The UEL method required only 1427 s to obtain a solution. Although some accuracy is compromised, the time necessary to complete this analysis is significantly reduced with the use of the surface element. Since a larger fraction of the total nodes in the shear bar model are located in the mold, a greater reduction in time is obtained, when compared with the L-shaped bar analysis.
546
A. Chang, J. Dantzig / Appl. Math. Modelling 28 (2004) 533–546
4. Conclusions Simulations to predict residual stresses in casting are usually performed including the mold in the analysis. An overwhelming amount of computation time can be required to obtain both the thermal and mechanical solution in the mold. However, the mold elements are really useful only for providing accurate boundary conditions on the surface of the cast part. The development of a surface element by Metzger et al. [3] provided a method to eliminate the mold nodes from the stress analysis and replace them with surface elements that provide an equivalent normal force to the part. The contact algorithm used in that work failed for higher stiffness materials, such as those commonly used for chills. This problem was alleviated by adding a new contact condition that adequately describes higher stiffness materials. The new contact algorithm was implemented using the user element function in ABAQUSTM . Two cases were examined to test the ability of the user element to accurately predict the deformation and plastic strain determined by using a full mold solution. The general deformed shape can be predicted quite well and the overall plastic strain can be calculated accurately. Local errors at a few points can be large. This is accomplished while reducing the time necessary to obtain these solutions significantly. It appears that the surface element approach has some difficulties predicting bending and rotation of the cast part. Some penetration is also overestimated due to the significant amount of compaction necessary for the initiation of the maximum stiffness in the high stiffness contact model.
References [1] M. Bellet, F. Decultieux, M. Menai, F. Bay, C. Levaillant, J. Chenot, P. Schmidt, I. Svensson, Metall. Trans. B 27B (1996) 81. [2] M. Bellet, C. Aliaga, O. Jaouen, in: P.R. Sahm, P.N. Hansen, J.G. Conley (Eds.), Modeling of Casting, Welding and Advanced Solidification Processes––IX, Shaker verlag, Aachen, 2001. [3] D. Metzger, K.J. New, J. Dantzig, Appl. Math. Model. 25 (2001) 825. [4] A.B. Chaudhary, K. Bathe, Comp. Struct. 24 (6) (1986) 855. [5] J.T. Oden, N. Kikuchi, Int. J. Numer. Meth. Eng. 18 (1982) 701. [6] S.K. Arya, G.A. Hegemier, J. Struct. Div. ASCE 108 (1982) 327. [7] D.L. Metzger, Development of a Sand Surface Element for Improved Residual Stress Determination, MasterÕs Thesis, University of Illinois at Urbana-Champaign, 1999. [8] A. Jacot, S. Cockcroft, Metall. Mater. Trans. A 31A (2000) 1201. [9] J.A. Dantzig, S.C. Lu, Metall. Trans. B 16B (2) (1985) 195. [10] J.A. Dantzig, J.W. Wiese, Metall. Trans. B 16B (2) (1985) 203. [11] J.W. Wiese, J.A. Dantzig, Appl. Math. Model. 12 (1988) 213. [12] K.J. New, An Improved Method for Simulating the Contact Resistance of Sand Molds on Gray Iron Castings, MasterÕs Thesis, University of Illinois at Urbana-Champaign, 1997.