INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING, VOL. 39, 2837-2862 (1996) SENSITIVITY ANALYSIS FORMULATION FOR THREE-DIMENSIONAL CONDUCTION HEAT TRANSFER WITH COMPLEX GEOMETRIES USING A BOUNDARY ELEMENT METHOD SEONG JIN PARK AND TAI HUN KWON Department of Mechanical Engineering, Pohang University of Science and Technology, Pohang 790-784, South Korea SUMMARY A special boundary integral formulation had been proposed to analyse many engineering problems of conduction heat transfer in complex three-dimensional geometries (closely spaced surface and circular hole in infinite domain or simple modification of it) by Rezayat and Burton. One example of such geometries is the mold sets in the injection molding process. In this paper, an efficient and accurate approach for the design sensitivity analysis (DSA) is presented for these kinds of problems in the similar complex geometries using the direct differentiation approach (DDA) based on the above special boundary integral formulation. The present approach utilizes the implicit differentiation of the boundary integral equations with respect to the design variables (radii and locations of circular holes) to yield the sensitivity equations. A sample problem (heat transfer of injection molding cooling system) is solved to demonstrate the accuracy of the present sensitivity analysis formulation. Although the techniques introduced here are applied to a particular problem in heat transfer of injection molding cooling system, their potential application is quite broad. KEY WORDS: special boundary integral formulation; design sensitivity analysis (DSA); direct differentiation approach (DDA); three-dimensional conduction heat transfer 1. INTRODUCTION Many engineering problems involve heat transfer, stress analysis, fluid mechanics, electronics, acoustics, etc., in regions with complex three-dimensional geometries (closely spaced surface and circular hole in infinite domain or simple modification of it). Some examples of such geometries are the mold sets for the process of injection molding of plastics and the region of fracture mechanics problem for crack tip stress analysis, etc. The present study is concerned about a particular heat transfer in the injection molding process. In this process, molten plastic is injected into the cavity of the closed mold, solidified in place, and then mechanically ejected. Heat is removed from molds mostly by water flowing in channels drilled in the mold.' Rezayat and Burton2 proposed the special boundary integral formulation to analyse the heat transfer in injection molds of these complex geometries. The proposed formulation consists of the boundary formula, its gradient and the exterior formula to calculate the heat transfer between a thermoplastic part and the corresponding mold. Design sensitivity analysis (DSA) will ultimately play an important role in the design optimization process after the above analysis has been achieved. Design sensitivities represent the gradients of variables of interest with respect to design variables (DVs). This information is required by non-linear optimization algorithms with first-order methods to determine the search direction in CCC 0029-598 11961162837-26 0 1996 by John Wiley & Sons, Ltd. Received 10 April 1995 Revised 22 September 1995 2838 S. J. PARK AND T. H. KWON which to proceed in altering the existing design of an object as it evolves towards an optimal configuration. These algorithms change the configuration of the object incrementally and require a number of such increments or iterations to obtain the optimal configuration. Inaccurate evaluations of design sensitivities may lead to an increased number of iterations and may even render the effort unsuccessful. The determination of accurate design sensitivities is thus of crucial significance and has attracted the attention of numerous researchers. For example, there are some research works on this subject of optimization problems in manufacturing processes of plastic materials: Barone and Caulk3 and Matsumoto and c o - w o r k e r ~for ~ ~thermal ~ design of injection molds, Forcucci and Kwon6 and Matsumoto and co-workers7~*for thermal design of compression molds, etc. The Boundary Element Method (BEM) has special advantages when applied to design sensitivity analysis as pointed out by, among others, Barone and Yang' and Saigal and co-workers.1° Two methods, namely the direct differentiation approach (DDA) and the adjoint structure approach (ASA), emerge as the most significant ones for the evaluation of design sensitivities. Both the approaches have recently been developed in the context of boundary element method for elastic systems: for examples, see the references by Barone and Yang' for two-dimensional systems, Saigal and co-workers" for axisymmetric systems, etc. These approaches have more recently been applied to several thermal system such as steady-state diffusion problems: for examples, Reference 11 for two-dimensional and axisymmetric heat diffusion problems, Reference 12 for three-dimensional heat conduction problems, etc. In this paper, an efficient and accurate approach for the design sensitivity analysis is proposed for the heat transfer of injection molding process using the direct differentiation approach based on the special boundary integral formulation proposed by Rezayat and Burton.2 The major advantage of this approach lies in the fact that the accurate sensitivity analysis results are obtained simultaneously for all design variables. We can generally consider radii and locations of circular holes (cooling channels in an injection mold) as design variables in many engineering problem^.'^.'^ Using an example, the numerical results of the sensitivity analysis simulation developed in the present study are amply discussed in terms of accuracy and computational efficiency. The accuracy of the present formulation is demonstrated through favourable comparisons between the sensitivity analysis results obtained from the present direct differentiation approach and those obtained by a numerical derivative approach using the forward finite difference method (FDM) in this example. 2. BOUNDARY INTEGRAL FORMULATION FOR HEAT TRANSFER ANALYSIS 2.1. Governing equation The configuration of a very simple injection mold is schematically shown in Figures 1 and 2. 0, Sp, I?, Sc, SE,and ii denote the mold region, the part surface, the mid-surface of part, the cooling channel surfaces, the mold exterior surfaces, and the outward unit normal vector in part, respectively. We can describe the part surface by only the mid-surface of the part I? instead of the whole surface (positive and negative surface as indicated by Sp' and Si- in Figure 2.) because plastic parts are usually thin (sheet-like). It is interesting to note that in most injection molding applications, heat loss through exterior surface is very small (typically less than 5 per cent). In such cases, even though the normal treatment of the mold exterior surface can be applied, the mold exterior surface may be treated as an infinite adiabatic sphere without significant change of re~u1ts.l~ It may be noted that this approximation does not require modelling of the mold exterior surface and thus it leads to substantial savings in the memory requirement and CPU time. Thus, the given geometry consists of closely spaced surfaces and circular holes in an infinite domain. The governing differential equation for the heat transfer of injection mold cooling stage under the 2839 SENSITIVITY ANALYSIS FORMULATION Figure 1. A simple example of mold bounded by r, SC and SE ~moldaurfaas, region n Figure 2. Schematic diagram of a cross-sectional mold for notations so-called cycle-averaged concept’3% l 4 may be written as (steady-state heat conduction equation) V 2 T = 0 in 0 (1) + + where T is the temperature. The boundary conditions on the boundary S = Sp SC SE are given as -k- aT = h ( T - Tb) on SC an (2 ) aT on SE an where h, k , h, and Tb denote the outward unit normal vector, the thermal conductivity of the mold material, the heat transfer coefficient in the cooling channels and the coolant bulk temperature, respectively, and (aT/an)o is assumed to be specified varying with the location in the part. (As a matter of fact, (dT/an)o is to be determined as a part of solution of the system. However, the -=0 2840 S. J. PARK AND T. H. KWON following method for the sensitivity analysis is valid regardless of this assumption concerning (aT/an)O.) 2.2. Boundary integral formulation for analysis A standard boundary element formulation for the three-dimensional Laplace equation given by equation (1 ) governing the conduction heat transfer in a mold leads to” Here x and 5 are position vectors in space, r = 15 - X I , and a denotes a solid angle formed by the boundary surface at the position x . Note that ci = 211 on the smooth boundary surface, and a = 471 or 0 at internal points or external points of the domain, respectively. For any two closely spaced surfaces such as the mold cavity surfaces, equation (3) leads to redundancy in the final system of linear algebra equations; in such cases, a modified procedure as described in Reference 2 needs to be used. According to this modification, a mid-surface is considered rather than two closely spaced surfaces. For example, in the case of mold cavity surfaces S p which consist of two closely spaced surfaces, Sp’ and S ,; it is represented by the cavity midsurface I? as shown in Figure 2. It might be noted that each mid-surface element has two degrees of freedom (e.g. temperature on Sp’ and S;) in order to derive the extra equation corresponding to the additional degree of freedom, we also used an approach similar to that used for crack analysis in fracture mechanics.2 For this purpose, a derivative of equation (3) with respect to the normal direction, 3 ( i + in Figure 2), at the mid-surface element is taken. See Reference 2 for further details of this modified approach for cavity surfaces. For cooling channel (circular hole) surfaces (SC as shown in Figures 1 and 2), we use a special formulation based on the line-sink approximation?. 3*16- I8 This formulation avoids discretization of the circular channels along the circumferential direction and thus saves large amount of computer memory and time. When we approximate the exterior surface by an adiabatic sphere at infinite, the integral parts of equation (3) and its gradient over SE do not need to be evaluated because the following integrals hold:2,l9 [!g - $(!)TI dS =41cT, (4) In equation ( 4 ) , T, denotes the temperature on the exterior surface (at infinite point). The followings are the final BEM formulae for these modifications: for a point x on the mid-surface of the part (I‘) we have the following pair of integral equations: + 4xT, (6) SENSITIVITY ANALYSIS FORMULATION 284 1 and for a point on the axis of the cylindrical segment of the cooling channels we have the following integral equation: In equations (6)-(8), N is the total number of cooling channels, uk is the radius of kth cooling channel, 1 is the local co-ordinate in the axial direction of each cooling channel, and 0 is the local circumferential co-ordinate of each cooling channel as shown in Figure 3. It is noted again that ( A + in Figure 2) in equation (7) is the outward unit normal vector in the plus plane at the source point (x). Figure 3. Co-ordinate transformation rule in one cooling channel element for sensitivity analysis with respect to location of cooling channel 2842 S. J. PARK AND T. H. KWON 3. BOUNDARY INTEGRAL FORMULATION FOR SENSITIVITY ANALYSIS 3.1. Design variables The boundary element method is also a strong tool in design sensitivity analyses when the problem can be treated without the domain discretization and the design variables can be defined on the boundary. In the mold cooling system design, the size and the arrangement of the channels can be treated as the design variables defined only on the boundary. Therefore, the boundary element method can be effectively applied to the sensitivity analyses of this kind of heat transfer problems. The sensitivity values of the temperature and the heat flux with respect to these design variables give us valuable information for the design optimization. Therefore, in this paper, we consider the radii and the locations of the cooling channels as design variables. (The coolant flow rate and the coolant bulk temperature can be also considered as design variables in injection mold problems, but will not be discussed here.) In other engineering problems, the size and the arrangement of the circular holes may be very interesting design variables for designers of those systems. 3.2. Boundary integral formulation for sensitiviry analysis The design sensitivity analysis consists of determination of the variation of an objective hnction with respect to a variation of each design variable. In the present study, only the sizes and the locations of the cooling channels are considered as design variables. To determine the sensitivity of a heat transfer system with respect to each design variable, the implicit differentiation of equations (6)--(8) with respect to each design variable is performed. This leads to the following design sensitivity equations for each design variable. First, we consider the radius of jth cooling channel, u j , as a design variable: for the part surface, the implicit derivatives of equations (6) and (7) with respect to the radius of jth cooling channel, a,, give us the following pair of integral equations for the sensitivity analysis: - aT+(x) tl- aaj + tl+- aT-(x) aa, SENSITIVITY ANALYSIS FORMULATION 2843 (10) for the ith cooling channel surface, the implicit derivative of equation (8) with respect to the radius of jth cooling channel, u j , gives us the following integral equation for the sensitivity analysis: + 4K- a r m dUj (11) In equations (9)-(11), the third integral terms appear because r is the function of the design variable uj and the domain of integral (the surface of jth cooling channel) also changes. Next, we consider the location of jth cooling channel as a design variable (a vector): for the part surface, the implicit gradient of equations (6) and (7) with respect to the location of the jth cooling channel gives us a pair of integral equations for the sensitivity analysis as follows in a gradient form: cr-VT+(x) + a+VT-(x) + 4aVTm 2844 S. J. PARK AND T. H. KWON - 2 (a(:)) av an+ 1 (VT' - V T - ) dS(Q for the ith cooling channel surface, the implicit gradient of equation (8) with respect to the location of the jth cooling channel gives us the following integral equation for the sensitivity analysis, depending on the value of i: when i # j , 0= A [:{ (g) + g)} & V + 4nV T , when i = j, V( - 1 ( ! ) ( V T + - V T - ) dS(5) SENSITIVITY ANALYSIS FORMULATION 2845 + 4nVT, (15) In the above equations (12)-( 15), V t and 0, specifically denote gradient with respect to 5 and x, respectively, depending on which is the design variable, where r = 15 - XI. Because the location is a vector, the co-ordinate transformation rule should be introduced for the sensitivity analysis with respect to the location of the cooling channel. The co-ordinate transformation rule between global co-ordinate system and local co-ordinate system is illustrated in Figure 3 and briefly described in Appendix I. 3.3. Solution procedure Once the integrals on each element are calculated, the discretized boundary element formulae for the analysis can be manipulated to obtain the following matrix equation: where [Hi,] and [ Cij] are hnctions of geometry of boundary surfaces. Some of [Hi,] and [ Gi,] are singular integral, which can be evaluated analytically as described in Appendix 11. Also, the integrals over 8 in equations (6)-(8) can be evaluated in closed forms using the complete elliptic integrals, See Appendix I1 for details. Other integrals can be evaluated by the Gaussian quadrature integration method. Next, boundary conditions can be introduced into equation (16) to obtain a system of linear algebraic equations: [ ~ i j ] ={fiI (17) where [Aii] and { f i } reflect boundary conditions. For this particular problem, temperature T is taken to be an unknown on each element with dT/dn being eliminated with the help of a mixed boundary condition. This system of equations can be solved either by LU-decomposition through Gaussian elimination or by an under-relaxation iterative method. In similar manner, the discretized boundary element formulae for the sensitivity analysis can be manipulated to obtain the following form of equation: [ ~ i j ]{Tj,x) = {fi.x) - = {f;,x) [ ~ i j , {~ T] j ) (18) where X is a design variable and, 'X' stands for a partial derivative with respect to X. The integrals in above equations can be evaluated in the similar manner to the above analysis case. See Appendix I11 for the detailed integral formulae for various terms appearing in equations (9)(15). It should be noted that [Ai,] in equation (18) is exactly the same as that in equation (17) for the analysis, and that only the forcing terms { f ; , x } are to be evaluated for each design variable. Therefore, once { f i , x } are evaluated for all design variables, the sensitivity analysis results, i.e. { T ~ J } can , be determined simultaneously with multiple forcing matrices for all design variables, thus saving the computational time. This fact is the most fascinating point of the present sensitivity analysis approach. In this regard, the major advantage of LU-decomposition solver lies in the fact that the LU-factorization of the matrix [ A i j ] , formed during the analysis, can be multiply reused in this sensitivity analysis. On the other hand, an iterative method 2846 S. J. PARK AND T. H. KWON which takes less computer memory than a direct solver can also efficiently evaluate {Tj,x} for all design variables, thus would be very useful when one needs large number of boundary elements. 4. RESULTS AND DISCUSSION A simple injection mold cooling problem as shown in Figure 4 is presented here as an example to demonstrate the accuracy and the computational efficiency of the sensitivity analysis method formulated in the present study. The mold has two cooling channels, and the dimension of the part is 200 x 200 x 2 (unit: mm). The global co-ordinate system was chosen such that the x-co-ordinate is in the axial direction of the cooling channels, the y-co-ordinate is in the thickness direction of the part, and the z-co-ordinate is in the spanwise direction of the part, as shown in Figure 4. The sizes and the locations of cooling channel #1 and cooling channel #2 are made intentionally unsymmetrically to show their effects on the cooling analysis. The discretization for the boundary element analysis is shown in Figure 5. In Figure 5 , each asterisk represents a node. We used triangular elements for the part surface and line elements for the cooling channel surface. Number of elements in the part and two cooling channels are 558 and 30, respectively. It d i n g line H (rsdlur: 0.6 cm) 3.5 cm + plane t h i i w 0.2 an oMin COOliiline82 (radius: 0.4 c m I ) d b 20 cm Figure 4. Diagram for geometry of an example problem SENSITIVITY ANALYSIS FORMULATION 2847 Figure 5. Boundary element mesh for the example might be mentioned that each element has a constant temperature and a constant heat flux in the boundary element analysis in the present study. The thermal conductivity of mold material used in this example is 100 W/m .K.For this example, we used the following boundary conditions: = 10 OC/m, (a) for part surface: uniform temperature gradient (&'/an), (b) for both cooling channels: Tb = 20 "C, h = lo4 W/m2 . K. Figure 6 illustrates analysis results for the part surface temperature distribution in each plane. The elapsed CPU time for this analysis is 0:28:19 (h:min:s) in the machine SUN SPARC 10 (22.9 MFLOPS/102 mips). (All computations were run on the same machine with a UNIX operating system using the same Fortran 77 compiler and complier options.) We have used four significant digits for each evaluated number (possible in single precision calculation). What we mean by the number of significant digits, n, is that we used the following the maximum norm convergent criteria in using an iterative method: ""I TFew - T.Old cold 1 'lo-" where Tiis an unknown, i.e. temperature, in each element. The maximum surface temperature in each plane is located in the position which is opposite to two cooling channels in the sense of both y-co-ordinate and z-co-ordinate, respectively as indicated in Figure 6. Moreover, both the minimum temperature and maximum temperature in the plus plane are higher than those of the minus plane, thus cooling channel #1 has more dominant cooling effect than cooling channel #2, as expected. As far as the sensitivity analysis is concerned, we consider the radius, the y-co-ordinate and the z-co-ordinate of each cooling channel as design variables because the effect of the x-co-ordinate is negligible. Therefore, six design variables exist (three design variables for each cooling channel) in this example. The elapsed CPU time for the sensitivity analysis using the direct differentiation approach is 2:27:00 in the same machine, thus the total elapsed CPU time is 2:55:19 (one analysis and one sensitivity analysis). In the sensitivity analysis, we have used three significant digits for each evaluated number (possible in single precision calculation). In this kind of geometries, the exact solution does not exist. Therefore, in order to check the accuracy of the sensitivity analysis results, we should directly compare sensitivity analysis results obtained from the direct differentiation approach with those from another numerical method, namely a forward finite difference method. Concerning the accuracy of the forward finite difference method, we have investigated the effects of the number of significant digits used in the computation and the magnitude of difference 2848 S. J. PARK AND T. H. KWON max: 146.46 15: 14: 13: 12: 11: 10: 9: 8: 7: 6: 5: 143.39 140.33 137.27 134.21 131.15 128.08 125.02 121.96 118.90 115.84 112.78 4: 109.71 3: 106.65 2: 103.59 1: 100.53 min: 97.47 max: 138.75 15: 136.06 14: 13: 12: 11: 10 9: 8: 7: 6: 5: 4: 3: 2: 1: 133.36 130.67 127.97 125.27 122.58 119.88 117.18 114.48 111.79 109.09 108.40 103.70 101.01 98.31 min: 95.61 Figure 6. Surface temperature distribution in part (unit: "C): (a) plus plane, (b) minus plane in design variables (AX) on the sensitivity analysis results. Figure 7 illustrates those effects on the sensitivity analysis with respect to a design variable, the radius of the cooling channel #1, R1. In Figure 7, the ordinates of both plots (a) and (b) represent the root-mean-square value of SENSITIVITY ANALYSIS FORMULATION 2849 Figure 7. Some criteria for finite difference method: (a) criterion on the number of the significant digits, (b) criterion on the magnitude of design variable difference d ( T + - T-)/dRI, i.e. 1 with the abscissa of plot (a) and (b) representing the number of significant digits and the magnitude of difference, respectively. The magnitude of difference is set to be 0.01 percent (0.00006cm / 0.6cm) to obtain results in Figure 7(a), while the number of significant digits is set to be 10 for Figure 7(b). According to Figure 7, it was found that at least six-digit significant figure and 0.01 per cent difference are required for accurate computations for the sensitivity analysis. Therefore, 2850 S. J. PARK AND T. H. KWON max: -30.41 15: -31.54 14: -32.68 13: -33.81 12: -34.94 11: -36.08 10: -37.21 9: -38.35 8: -39.48 7: -40.62 6: -41.75 5: -42.89 4: -44.02 3: -45.15 2: -46.29 1: 47.42 min: -48.56 max: -38.70 15: 14: 13: 12: 11: 10: 9: 8: 7: 6: 5: 4: 3: -40.29 -41.89 -43.48 -45.08 46.67 -48.27 49.86 -51.46 -53.05 -54.66 -56.24 -57.04 -59.43 2: -61.03 1: -62.62 min: -64.22 Figure 8. Sensitivity results of part surface temperature by the direct differentiation approach with respect to radius of cooling channel #1 (unit: OC/cm): (a) plus plane, (b) minus plane we have used a double-precision calculation in order to satisfy the above requirements when the forward finite difference method was used for the sensitivity analysis. The sensitivity analysis by the forward finite difference method required the total elapsed CPU time of 10:14:23 (the elapsed CPU time is 1:17:46 for each run of the analysis program and total 7 runs of analysis program 285 1 SENSITIVITY ANALYSIS FORMULATION max: 30.63 15: -31.81 14: -32.98 13: -34.15 12: -35.33 11: -36.50 10: -37.68 9: 8: 7: 6: 38.85 -40.02 -41.20 -42.37 5: -43.55 4: 44.72 3: 45.90 2: 47.07 1: -48.24 min: -40.42 max: -39.20 15: -40.85 14: -42.50 13: -44.16 12: -45.81 11: -47.46 10: -40.11 9: -50.76 8: -52.41 7: -54.06 6: -55.71 5: 37.36 4: 4 . 0 1 3: -60.66 2: -62.31 1: -63.96 min: -65.61 Figure 9. Sensitivity results of part surface temperature by the forward finite difference method with respect to radius of cooling channel #1 (unit: OCkm): (a) plus plane, (b) minus plane are needed for 6 design variables), thus requiring about 3 times as much as that for sensitivity analysis by the direct differentiation approach in this particular example. It should be noted that the ratio of the total elapsed CPU time by the forward finite difference method to that by the direct differentiation approach approximately linearly increases with the number of design variables (i.e. cooling channels). As shown in Figure 7, sensitivity analysis results by the forward finite difference 2852 S. 1. PARK AND T. H. KWON 9 max: 3.49 L 1 15: 14: 13: 12: 11: 10: 9: 8: 7: 6: 5: 4: 3: 2: 1: 3.32 3.14 2.97 2.80 2.62 2.45 2.27 2.10 1.92 1.75 1.57 1.40 1.23 1.05 0.88 min: 0.70 max: 15: 14: 13: 12: 11: 10: 9: 8: 7: 6: 5: 4: 3: 2: 1: min: 10.93 10.30 9.66 9.02 8.38 7.74 7.10 6.46 5.83 5.19 4.55 3.91 3.27 2.63 2.00 1.36 0.72 Figure 10. Sensitivity results of part surface temperature by the direct differentiation approach with respect to y-co-ordinate of cooling channel #I (unit: OC/cm): (a) plus plane, (b) minus plane method converge to those by the direct differentiation approach as the number of significant digits increases and the magnitude of difference decreases, which in turn proves the accuracy of both approaches. Thus, we can conclude that the direct differentiation approach is more accurate and faster than the forward finite difference method. 2853 SENSITIVITY ANALYSIS FORMULATION I I - I max: 3.68 15: 14: 13: 12: 11: 10: 9: 8: 7: 6: 5: 4: 3: 2: 1: mln: max: 15: 14: 13: 12: 11: 10: 9: 8: 7: 6: 5: 4: 3: 2: 3.50 3.31 3.13 2.94 2.76 2.58 2.30 2.21 2.03 1.84 1.66 1.48 1.29 1.11 0.92 0.74 11.51 10.84 10.17 9.50 8.82 8.15 7.48 6.81 6.14 5.47 4.79 4.12 3.45 2.78 2.11 1: 1.43 min: 0.76 Figure 1 1. Sensitivity results of part surface temperature by the forward finite difference method with respect to y-co-ordinate of cooling channel #I (unit: OCkm): (a) plus plane, (b) minus plane For each design variable of cooling channel #I, sensitivities of part surface temperature with respect to the radius, the y-co-ordinate, the z-co-ordinate both by the direct differentiation approach and by the forward finite difference method are shown in Figures 8 and 9, Figures 10 and 11, and Figures 12 and 13, respectively. (Six-digit significant figure and 0.01 per cent difference in design 2854 S. J. PARK AND T. H. KWON max: 15: 14: 13: 12: 11: 10: 9: 8: 7: 6: 5: v I 4: 3: 2: 1: 1.26 1.08 0.91 0.73 0.56 0.38 0.21 0.03 -0.15 -0.32 -0.50 -0.67 -0.85 -1.02 -1.20 -1.37 min: -1.55 I max: 15: 14: 13: 12: 11: 10: 9: 8: 7: 6: 5: 4: 3: 2: 1: 5.50 4.92 4.33 3I75 3.17 2.59 2.00 1.42 0.84 0.26 -0.32 -0.91 -1.49 -2.07 -2.65 -3.24 min: -3.82 Figure 12. Sensitivity results of part surface temperature by the direct differentiation approach with respect to z-co-ordinate of cooling channel # I (unit: OC/cm) (a) plus plane. (b) minus plane variables are adopted in the finite difference calculations.) These figures reveal that the computed sensitivity results are in an excellent agreement with each other. The maximum relative error found in this example was of about 6-59 per cent. (The relative error on each element is defined as a ratio of the absolute value of the difference between sensitivity result by the forward finite difference 2855 SENSITIVITY ANALYSIS FORMULATION 1 max: 15: 14: 13: 12: 11: 10: 9: 8: 7: 6: 5: 4: 3: 2: 1: 1.32 1.13 0.95 0.76 0.58 0.40 0.21 0.03 -0.15 -0.34 -0.52 -0.70 -0.89 -1.07 -1.26 -1.44 min: -1.62 max: 5.76 15: 5.15 14: 4.54 13: 3.93 12: 3.32 11: 2.71 10: 2.10 9: 1.49 8: 0.88 7: 0.27 6: -0.34 5: -0.95 4: -1.56 3: -2.17 2: -2.78 1: -3.39 min: -4.00 Figure 13. Sensitivity results of part surface temperature by the forward finite di..-rence method with respect to z-coordinate of cooling channel # 1 (unit: OC/cm) (a) plus plane, (b) minus plane method and that by the direct differentiation approach to that by the direct differentiationapproach.) The maximum error may be reduced by increasing the number of significant digits and reducing the magnitude of difference in the forward finite difference method. From the above results, we have found the following characteristics of each design variable: 2856 s. J. PARK AND T. n. KWON (a) An increase of radius of cooling channel #1 makes the part surface temperature decrease mainly due to the increase of surface area of the cooling channel. (b) An increase of y-co-ordinate of cooling channel #1 make the part surface temperature increase due to the increase of the distance between the part and cooling channel #l. (c) As z-co-ordinate of cooling channel #1 increases (i.e. moving in the downward direction in Figure 4), the part surface temperature increases in the portion of the part where the distance between the part position and cooling channel #1 increases while it decreases in the other portion. To no surprise, the line of no change of part surface temperature, that is, zero sensitivity coincides exactly with the location of cooling channel #l. It may just be mentioned that we have also compared sensitivity analysis results in the case of each design variable of cooling channel #2 in the similar manner to the case of cooling channel #1, and obtained the similar accuracy and results to the case of cooling channel #l. 5 . CONCLUDING REMARKS An efficient and accurate approach for the design sensitivity analysis of steady-state heat conduction-type problem in complex three-dimensional geometries using boundary element method was presented. The presented formulae are based on the implicit differentiation of the special boundary integral equations with respect to the design variables (radii and locations of circular holes) to obtain the sensitivity equations. Using an illustrative example problem, the numerical results of the sensitivity analysis developed here is amply discussed in terms of accuracy and computational efficiency. The accuracy of the present formulation is demonstrated through comparisons between the present direct differentiation approach and the forward finite difference method. It was concluded that the former is faster and more accurate than the latter. These programs will be very useful for designers to obtain optimal configuration of radii and locations of circular holes (cooling channels in the injection mold) to minimize an objective function in many engineering problems. This optimization problem may be automated and solved numerically using any standard optimization software, for example, CONMINZowhich is based on Davidon-Fletcher-Powell method2' treating the constraints via an Augmented Lagrangian Multiplier Method. Finally, the procedures and techniques proposed in the present paper can be applied to similar three-dimensional heat conduction problem. The extension of these techniques to engineering problems which could be represented by the Laplace operator should be straightforward, because the kernels of the integral equations for heat transfer and those engineering problems are fundamentally similar. ACKNOWLEDGEMENTS Portions of the research discussed herein have been supported by grants from R & D Center in Samsung Electronics (Consumer Electronics Business) and Korea Science and Engineering Foundation (KOSEF). The authors gratefully acknowledge these financial supports. APPENDIX I Co-ordinate transformation rule The following orthogonal matix represents the co-ordinate transformation rule used for the sensitivity analysis with respect to the location of the cooling channel. The basis {i,j,i } represents 2857 SENSITIVITY ANALYSIS FORMULATION the local co-ordinate system in each line element of the cooling channel and the basis { f , j , k } represents the global co-ordinate system as shown in Figure 3. All symbols used in the following orthogonal matrix are shown in Figure 3: where P I = P, - x , Pz = P, - y , and P3 = P, - z In the special case that P is on the axis of the same element as Q, the above co-ordinate transformation rule can be simplified to permit the following form (in representing v, unprimed system { v,, vy,vz} for the global co-ordinate system and primed system { vi, v:, v:} for the local co-ordinate system): vii+ v:j = (vx - v:lx)i+ ( v y - v i l y ) j + (vz - V;I,)R v; L = v:(lxi + lJ + 1,R) where L = i = 1,i + l y j + lzK 1 - v, = 1 - 3 = lxvx + lYVY + l,v, APPENDIX I1 Integral formulae for heat transfer analysis Singular integrals in triangular element.2 The singular integrals over a flat area A in triangular element can be evaluated analytically. The symbols used in the following closed forms are shown in Figure 14: 1 + -1nl r3 I tan ((02 -k a3)/2) Wa3/2) I 1 12 + r tan ((03 -k ai )/2} tan ( aI /2) In the above equations, A denotes the area of the triangular element. I] 2858 S. 1. PARK AND T. H. KWON Figure 14. Notations for geometric information of triangular element Closed form of integrals in line element. The integrals over 6 can be evaluated in closed forms. In the special case that P is on the axis of the same element as Q, the integrals over 0 and 1 can be also evaluated in closed forms. The symbols used in the following definitions and closed forms are shown in Figure 3. We first define some fkctions used in evaluating integrals over 8. To avoid conflicts with other symbols in the paper, we shall use non-standard, calligraphic symbols for the following functions. X and & denote the complete elliptic integrals of the first and second kind, respectively, defined as X ( m )= de JiYza Lnj2 &(m) = J 0x ’ 2 d C Z & d ~ 9, 9, 9 and 9’denote simple modifications of X and & and can be evaluated using the derivatives of X and 8: where + m = B ~ / ( A B~ ~ ) SENSITIVITY ANALYSIS FORMULATION W , X,CY and 9’ are more general forms of tl = L t3 = ~ p , t2 = P,22, 9 and 9, respectively, d a 2 +z; Jq, + A’ = S: s;, B2 = 4aRp m = B ~ / ( A ~B ~ ) + In the analysis formulation, the following integral formulae are needed: a d 0 = 4 a W ( m ; 1) 12 = ~ 2 n - a( - ) 1a d 0 = 4 a X ( m ; s ~ , 2 R p ) an r 2859 2860 S. J. PARK AND T. H. KWON I3 = 14 = 2na I (--) ad0 1'" (% a a = 4a%(m;s3, -2avi) i (--))ud0= 12aY(m;~1~3,-2~~,-B~v~)+4av~%(m;-1,2) In the special case that P is on the axis of the same element as Q, the above formulae can be simplified to obtain closed forms of integrations over 1. In such a case, the following integral formulae are needed: JI = J2 = /dL 1'"1 1"12" (-) lLL2"av (--) a d o d l = 27ta In r a 1 a d o d l = 2n an r a 1 a dB d l = 27tav; J3 = J4 = l L 1 2 ' - ( - (a- ) ) a d 0 id l av an r (2+ 2) (i i) - =27ta2vi APPENDIX 111 Integral formulae for sensitivity analysis In the sensitivity analysis formulation with respect to the radius of the cooling channel, the following integral formulae are needed: all - 1 1 _ - - + 4a%(m; -SI, -2Rp) aa u -= '2 aa a 313 13 aa a ar4 14 - 12aq(m;s~,4Rps1,4R~) + 4a%(m; 1,0) - = - - 12aY(m;sls3, -2s4, -B2v:) -= aa + 4a v:%(m; 1, -2) - - 6 0 a ~ ( m ; s ~ s 3 , 2 s ~ s ~ , 4 -2B' R pRspv:) g, a +12a9(m;2slv: + s3, -2(3sl - Rp)v:, -8Rpv:) In the sensitivity analysis formulation with respect to the location of the cooling channel, the following integral formulae are needed: ~ t ( 1 1 ) = -Vx(Ii) = 4 a ~ ( m ; - s 1 , 2 a ) 2 ^ - 4 a s ~ ~ (Im , O;) ~ t ( 1 2 )= -v,(I~) = { 1 2 a ~ ( m-s:,~s:,B') ; i -I- 4 a ~ ( m1,-2)} ; 2^ - 12s2Y(m; s,, 2Rp, 0) + vt(I3) = -Vx(13) = - 1 2 a Y ( m ; s ~ s ~ , - 2 a ( s ~ v s3),4a2vi)zI ~ + 4 8 a 3 v ~ ~ ( m ; 0 , - 1 , 1 ) j -12as2Y(m;s3,-2av~,O) i + 4 a O%(m; 1,O) QU4) = -VXI4) = { -6Oab(m;sfs3, -2s:(s3 + av:),4as7,2aB2v:) + 12a?Y(m;slv: + s3, -2Fgr8av:)} + {-240a3v~Z(m;0,s~,3Rp- a , - 2 R p ) + 9 6 a 2 v ~ ~ ( m ; 01,, - l ) } j 2^ SENSITIVITY ANALYSIS FORMULATION 2861 Next, in the sensitivity analysis formulation with respect to the radius of the cooling channel for the special case that P is on the axis of the same element as Q,the following integral formulae are needed aa a Finally, in the sensitivity analysis formulation with respect to the location of the cooling channel for the same special case as above, the following integral formulae are also needed: V c ( J l )= -V,(JI) = -2za Vt(J2)= -V,(J2) = -272a2 ($ - $) REFERENCES I . K. J. Singh, ‘Design of mold cooling system’, in A. 1. Isayev (ed.), Injection and Compression Molding Fundamentals. Marcel Dekker, New York, 1987, pp. 567-605. 2. M. Rezayat and T. Burton, ‘A boundary-integral formulation for complex three-dimensional geometries’, Inr. j . numer. methods eng., 29, 263-273 (1990). 3. M. R. Barone and D. A. Caulk, ‘Optimal thermal design of injection molds for filled thermosets’, Polym. Eng. Sci., 25, 608617 (1985). 4. T. Matsumoto and M.Tanaka, ‘Optimum design of cooling lines in injection moulds by using boundary element design sensitivity analysis’, Finite Elements Anal. and Des., 14, 177-185 (1993). 5. T. Matsumoto, M. Tanaka and T. Takahashi, ‘Determination of cooling tempemture history in injection molding’, in H. D. Bui, M. Tanaka, M. Bonnet, H. Maigre, and M. Reynier (eds), Inuerse Problems in Engineering Mechanics, Bakema, Rotterdam, 1994, pp. 295-300. 6. S. J. Forcucci and T. H. Kwon. ‘A computer aided design system for three-dimensional compression mold heating’, ASME J. Eng. Indusrry, 111, 361-368 (1989). 2862 S . J. PARK AND T. H. KWON 7. T. Matsumoto, M. Tanaka and M. Miyagawa, ‘Boundary element system for mold cooling/heating design’, in C. A. Brebbia and J. J. Rencis (eds), Boundary Element XV-Vol. 2: Stress Analysis, Computational Mechanics Publications, Boston 1993, pp. 461-475. 8. T. Matsumoto, M. Tanaka and N. Ishii, ‘Boundary element system for design sensitivity analysis of compression mold heating’, in B. M. Kwak and M. Tanaka (eds), Computational Engineering, Elsevier Science Publishers, Amsterdam 1993, pp. 465-470. 9. M. Barone and R. J. Yang, ‘Boundary integral equations for recovery of design sensitivities in shape optimization’, AIAA J., 26, 589-594 (1988). 10. S. Saigal, J. T. Borggaard and J. H. Kane, ‘Boundary element implicit differentiation equations for design sensitivities of axisymmetric structures’, Inr. J. Solids Struct., 25, 527-538 (1989). 1 1. S. Saigal and A. Chandra, ‘Shape sensitivities and optimal configurations for heat diffusion problems: a BEM approach’, J. Heat Transfer, 113, 287-295 (1991). 12. K. G. Prasad and J. H. Kane, ‘Three-dimensional boundary element thermal shape sensitivity analysis’, Inr. J. Heat Mass Transfer, 35, 1427-1 439 (1 992). 13. K. Himasekhar, J. Lottey and K. K. Wang, ‘CAD of mold cooling in injection molding using a three-dimensional numerical simulation’, ASME J. Eng. Industry, 144, 213-221 (1992). 14. T. H. Kwon, ‘Mold cooling system design using boundary element method’, A S M E J. Eng. Industry, 110, 384-394 (1988). 15. C. A. Brebbia, J. C. F. Telles and L. C. Wrobel, Boundary Element Techniques-Theory and Application in Engineering, Springer, Berlin, 1984. 16. M. R. Barone and D. A. Caulk, ‘Special boundary integral equations for approximate solution of Laplace’s equation in two-dimensional regions with circular holes’, Q. J. Mech. Appl. Marh., 34, 265-286 (1981). 17. J. C. F. Telles, L. C. Wrobel, W. J. Mansur and J. P. S. Azevedo, ‘Boundary elements for cathodic protection problems’, in C. A. Brebbia and G. Maier (eds), Boundary Element VII-Proc. 7th Int. Con$, Springer, Berlin, 1985. 18. J. C. F. Telles, W. J. Mansur, L. C. M‘robel and L. C. Marinho, ‘Numerical simulation of a cathodically protected semi-submersible platform using the PROCAT system’, Corrosion J., 46, 5 13-5 18 ( 1990). 19. J. C. F. Telles, W. J. Mansur and L. C. Wrobel, ‘On boundary elements for external potential problems’, Mech. Res. Commun., 11, 373-377 (1984). 20. P. C. Haarhoff and J. D. Buys, ‘A new method for the optimization of a nonlinear function subject to nonlinear constraint’, Comput. J . , 13, 178-184 (1970). 21. R. Fletcher and M. J. D. Powell, ‘A rapidly convergent descent method for minimization’, Comput. J., 6, 163-168 (1963).

1/--страниц