Acta Polytechnica doi:10.14311/AP.2016.56.0472 Acta Polytechnica 56(6):472–477, 2016 © Czech Technical University in Prague, 2016 available online at http://ojs.cvut.cz/ojs/index.php/ap SIMULATION OF SURFACE HEATING FOR ARBITRARY SHAPE’S MOVING BODIES/SOURCES BY USING R-FUNCTIONS Sergiy Plankovskyy, Olga Shypul, Yevgen Tsegelnyk∗, Oleg Tryfonov, Ivan Golovin National Aerospace University, “Kharkiv Aviation Institute” named after N. Ye. Zhukovsky, Faculty of Aircraft Engineering, Department of Aircraft Manufacturing Technologies, Chkalova 17, 61070, Kharkiv, Ukraine ∗ corresponding author: y.tsegelnyk@gmail.com Abstract. The purpose of this article is to propose an efficient algorithm for determining the place of an action of a heat source with a given motion law for a body of an arbitrary shape using methods of analytical geometry. The solution to this problem is an important part of a modeling of a laser, plasma, ion beam treatment. In addition, it can also be used for mass transfer problems, such as simulation of coating, sputtering, painting etc. The problem is solved by the method of R-functions to define the shape of the test body and the heat source and the analytical determination zone shadowing. As an example, we consider the problem of using the method of ion cleaning parameters optimization considering temperature limitations. Application of the R-functions can significantly reduce the amount of computation with usage of the ray tracing algorithm. The numerical realization of the proposed method requires an accurate creation of a numerical mesh. The best results in terms of accuracy of determination the scope of the source can be expected when applying adaptive tunable meshes. In case of integration of the R-functions into the CAD system, the use of the proposed method would be simple enough. The proposed method allows to determine the range of the source by the expression, which is constructed only once for the body and the source of arbitrary geometric shapes moving in any law. This distinguishes the proposed approach against all known algorithms for ray tracing. The proposed method can also be used for time-dependent multisource with arbitrary shapes, which move in different directions. Keywords: numerical methods; moving heat sources; ray tracing; R-functions method. 1. Introduction Moving body/source heating analysis has an applica- tion in several manufacturing processes [1, 2]. In most well-known studies, in this formulation, the problem is considered for a circular heat source, which moves in a straight line [3, 4]. However, there are papers that outline the temperature distribution in a half-space, because of the complexity of the shape of the moving heat source [5], and because of the fact that a heat source moves according to a more complex laws [6]. The problem of heating bodies of finite dimensions by moving heat sources is considered in fewer studies. The objects of the study are likely to be the bodies of simple shapes (cylinder or parallelepipeds). For example, a number of studies were carried out by [7], to investigate the temperature distribution in a rotat- ing cylinder heated with the laser heat source. This problem may be associated with the calculation of laser hardening regimes, laser-assisted machining, etc. At the same time, it is necessary to calculate the temperature of an arbitrary shaped body caused by heating of a moving heat source. The law of motion and the shape of the heat source can be quite complex. Such problem is typical for cases, when body heat happens because of an energy flux action (radiation, heating, laser, ion or electron beam processing). In this case, the definition of the heated area’s borders is a complicated problem. In this paper, we propose an analytical method for the solution of this problem by using the R-functions. The geometric shapes of the body, the source and the law and their relative motion can be arbitrary. After definition of the action zone of the heat source , the further temperature calculation was determined numerically by a finite element method. 2. Mathematical formulation The problem with determining the action zone of the moving heat source is very similar to the problem with determining the shadow zones location that is typical to computer graphics.. Beginning with the studies [8, 9] such problems are solved by different algorithms of ray tracing. It required creating a lot of rays from points of the body surface in the direction of the radiation source. These tasks require large amount of computations. At the same time, for complex shape bodies, the precise definition of the shaded area is associated with considerable difficulties in case of the moving sources, which change the intensity of the radiation. We assume that the geometry of the source and the body is known. It can be independent of time, 472 http://dx.doi.org/10.14311/AP.2016.56.0472 http://ojs.cvut.cz/ojs/index.php/ap vol. 56 no. 6/2016 Simulation of Surface Heating Figure 1. Description of geometric shapes with R- functions. or change according to the known law. The source intensity is able to change arbitrary on time and on its cross section. The law of the source movement is also known. Solution of the problem is reduced to the construc- tion of a switch-function, which would be equal to 1 on the heated surface and to 0 on the rest of the surface of the part. For this purpose, it is enough to create a function that has the following properties:  ω(x,y,z,t) > 0 inside the heated surface part, ω(x,y,z,t) = 0 on the border of the heated surface part, ω(x,y,z,t) < 0 ouside the heeated surface part. (1) Such relation can be described using the R- functions [10]. The R-function method was devel- oped as an improvement of Ritz methods for solving boundary-value problems. It is known for its utiliza- tion of solving the heat conduction problems in [11]. However, in this paper it will be used only as a tool of analytical geometry. The R-functions are functions of continuous real arguments. Their sign is determined by the sign of the argument. If, instead the sign of the “+” and “−”, we use the values 0 and 1, R-functions are equivalent to some Boolean logic functions. Boolean function equiv- alent to a particular R-function, called the companion function. Almost in every CAD system, a geometric shape of a complex object is created by using Boolean operations with geometric primitives. These primitives can be defined by algebraic equations or inequalities. For example, inequality ω1 = (R2 −x2 −y2 ≥ 0) in the plane X0Y defines a circle of radius R centered at the origin. It is proved by [12] that if the geometry of the complex object is created by using Boolean opera- tions (∨,∧,¬) with geometric primitives which are described by the inequalities ωi ≥ 0, the replacement of companion Boolean function with R-functions, al- Figure 2. Definition of shadow zones. lows to obtain the inequality ω = f(ω1, . . . ,ωn) ≥ 0, with the properties (1). The simplest complete system of these R-functions is the following one: • f ≡−f (logical negation); • f1 ∧f2 = f1 + f2 − √ f21 + f22 (logical conjunction); • f1 ∨f2 = f1 + f2 + √ f21 + f22 (logical disjunction). Let it be required to create an expression ω(x,y) ≥ 0 for the area which is shown in Figure 1. As geometric primitives, the following items are selected: ω1 =( 1 −|x/a|− |y/b| ≥ 0 ) – the part of the plane inside a diamond with vertices (±a, 0), (0,±b); ω2 = (x2 + y2 − r2 ≥ 0) – the part of the plane outside of the circle with radius r centered at the origin. The shape of the region is determined by ω = ω1 ∧ ω2, or after replacing the Boolean operations with R-functions: ω = ( 1 −|x/a|− |y/b| + x2 + y2 −r2 − √ (1 −|x/a|− |y/b|)2 + (x2 + y2 −r2)2 ≥ 0 ) . Further, for simplicity, the symbols ∧R, ∨R will be used instead of the expanded form of R-function. Method of analytical determination of the coverage of the source considers the example of two-dimensional problem. We propose dependences that describe the geometric shape of the heat source ωsource ≥ 0 and body ωbody = R(ωi) ≥ 0, where R(ωi)– system of R-functions; ωi, i = 1, . . . ,N – geometric primitives (Figure 2). 473 S. Plankovskyy, O. Shypui, Y. Tsegelnyk et al. Acta Polytechnica Figure 3. Design scheme for the test problem. We propose the expression gi = −(ωi ∧R ωsource)2, which is equal to zero in the surface regions Γi, primitives fall into the domain ωsource and less than zero at all other points. We also use the expression gbody = −(ωbody ∧R ωsource)2, which is equal to zero for all Γi. Expression Ψi = 4H(gi) × H(gbody), where H (Heaviside function) is equal to one for all points on the surface of the body, lying on the Γi, and zero for the remaining points. It means that ray tracing are needed only for points, where Ψi is equal to one. Without loss of generality, we assume that the ac- tion of the source is directed along the axis ysource. We use lines fi = xsource −xpoint for points on the surface lying on Γi. Expression H(−f2i ) × Ψi ×|ysource| de- termines the distance between the source and the point Γi, which lies on fi. Point of the line fi, which lies closest to the source, will have coordinate yi,minsource = minbody ( H(−f2i ) × Ψi ×|ysource| ) . Seeking function equal to 1 in areas which are af- fected by the source and 0 for the rest of the Γi can be written as: Φ = 2H ( yi,minsource − Ψi ×|ysource| ) . (2) With the known law of motion of the body, connect- ing the coordinates (xbody,ybody) and (xsource,ysource), functions H(gbody) and Φ clearly defines the scope of the source for the body of arbitrary shape given by the expression ωbody. The method of creating Equation (2) may seems cumbersome. However, this expression is created only once. This distinguishes the proposed approach against all known algorithms for ray tracing. The proposed method can also be used for time-dependent multisource with arbitrary shapes which move in dif- ferent directions. At solving such problems using the finite element method, quality of domain finite element mesh has a serious influence on the accuracy of the source lo- cale. If the finite element mesh of the object is not performed accurately enough, the mesh generator de- termines the coordinates of the surface nodes with an error. In this case, the elements heated by source on the surface of the body are defined incorrectly. You can see such influence in the following test problem. Determination of the moving heat source site of action for complex shape body was solved as a test problem. A rectangular shape source with a cut performs a reciprocating motion along the axis xbody, and rotates around the axis ysource (Figure 3). For this problem, the coordinate system of the source and body are connected by expressions:  xsource = (xbody + l− 2l sin ϕ̇1t) cos ϕ̇2t +zbody sin ϕ̇2t, ysource = ybody, zsource = −(xbody + l− 2l sin ϕ̇1t) sin ϕ̇2t +zbody cos ϕ̇2t. Figure 4 shows the results of numerical calculations for two cases: when using automatic generation of unstructured mesh without additional conditions and with the mesh obtained by the forced association of mesh nodes to the body surface. Area of the source is defined incorrectly on curved surfaces and planes parallel to it for case one. 3. Numerical simulation result and discussion Simulation of gas turbine engine blade heating dur- ing ion sputtering has been done using the proposed method. The process is performed on an ion-plasma device and is a preparatory stage before coating. Sput- tering is a slow process, so to improve performance, a set of parts was treated at the same time. So batch of blades is set on a turntable which rotates with an angular velocity ϕ̇1 (Figure 5). For more uniform sputtering and coating, the blade also rotates around its own axis with an angular velocity ϕ̇2. Energy of the ion beam, sputtering time, the rotation speed of the table and of the blades must be assigned so as to provide a predetermined quality of treatment and to prevent the blade material from overheating above the phase transition temperature. For the heat resistant alloys, which are most commonly used in a turbine blade manufacturing, this temperature is equal to 1270 K. The intersection of the ion beam and the blade surface causes the blade heating. It was assumed that the ion flux occupies a cylindri- cal shape area with radius r, located perpendicularly to the plane Z0X, displaced along the axis z from the plane of the table on a distance h. Movement of ions occurs along a positive direction of the axis y (Figure 5). The blade heating occurs near the ion source. There- fore, the shape of the ion beam was set by the expres- sion ωsource = (r2 −x2 −(z−h)2)∧R−y. Coordinate system of the body and the source are related by the 474 vol. 56 no. 6/2016 Simulation of Surface Heating Figure 4. The results of numerical simulations using unstructured mesh without conditions (left) and mesh obtained by the forced association of mesh nodes from the body surface (right). expressions:  zsource = zbody, xsource = (xbody + R cos ϕ̇1t) cos(−ϕ̇2)t +(ybody + R sin ϕ̇1t) sin(−ϕ̇2)t, ysource = −(xbody + R cos ϕ̇1t) sin(−ϕ̇2)t +(ybody + R sin ϕ̇1t) cos(−ϕ̇2)t. Blade surface geometry was defined by point cloud. The heat flux is determined on the basis of the energy balance at the surface between input heat Qin and heat loss Qout. For plasma processes Qin described by the expres- sion [13]: Qin = Jrad,in + Jch + Jn + Jads + Jreact,in + Jext,in. Here Jrad,in is the heat radiation towards the sur- face; Jch is the power transferred by charge carriers (electrons and ions); Jn is the contribution of neutral species of the background gas and the neutral par- ticles; Jads and Jreact,in are energies released by an absorption or a condensation and the reaction energy of exothermic processes including molecular surface recombination; Jext,in is an input heat flux by exter- nal sources that influences the thermal balance of the substrate. Clearing by ion sputtering is carried out at low pressure with using a high-purity neutral gases, with- out additional heating and cooling. In this case, the expression for determining Qin can be simplified [13]: Qin = Jrad,in + Jch = σ(εT 4rad −εsT 4 s ) + jiEion + 4kcMiMs (Mi + Ms)2 sin2 θ2 eiVbias, where ε is the spectral emittance of the radiation source at a temperature Trad; εs represents the spec- tral absorbance of the substrate surface at a temper- ature Ts; σ denotes the Stefan–Boltzmann constant; ji is the ion flux density of the surface; Eion is the ionization potential of the incident ion; kc is the en- ergy transfer coefficient; Mi, Ms are masses ratio of the colliding particles (ion and surface atom); θ is the angle of incidence; ei is the ion charge; Vbias is the sum of the plasma potential and the substrate potential. The heat loss Qout of the substrate during plasma processing consist of the following terms [13]: Qout = Jrad,out + Jparticle + Jdes + Jreact,out + Jext,out. Here Jrad,out is the energy radiated from the substrate at a temperature Ts; Jparticle is the energy transport 475 S. Plankovskyy, O. Shypui, Y. Tsegelnyk et al. Acta Polytechnica Figure 5. The equipment for sputtering and coating – turntable with blades and design scheme for the task of the blade heating under ion sputtering. from the substrate due to sputtering of surface atoms and the secondary electron emission; Jdes is the energy sink due to the desorption of particles into the gas phase and the diffusion into the solid bulk; Jreact,out is the reaction energy of exothermic processes, including molecular surface recombination; Jext,out is the heat loss caused by an external cooling. Subject to the terms of the ion cleaning, the expres- sion for determining Qout can be written as: Qout = Jrad,out + Jparticle ≈ Jrad,out + Jsputtering = σ(εsT 4s −εenvT 4 env) + ksputjiEcoh, where εenv is the emissivity of the environment; Tenv is the environmental temperature (reactor walls, etc.); ksput is the sputtering coefficient dependent on the ion energy and angle of incidence; Ecoh is the cohesive energy of sputtered material’s atoms. The temperature field in the clearing blade was calculated during the simulation. For this, the tran- sient heat conduction equation was solved. The heat flux through surface of the blade was specified by the expression: W = σ(εsT 4s −εenvT 4 env) + Φ ( σ(εT 4rad −εsT 4 s ) + jiEion + 4kcMiMs (Mi + Ms)2 sin2 θ2 eiVbias −ksputjiEcoh ) , where Φ is a switch-function, which is determined by the expression (2). Dependence of the heat flux value on the angle be- tween the direction of ions flux and the blade surface is a feature of the process. This feature required calcu- lating the directional cos of the normal to the surface. Using R-functions makes it simple to solve this prob- lem. If the equation ωi ≥ 0 for geometric primitives is known, and for a complex domain the equation ω(ωi) ≥ 0 is constructed by using R-functions, than Figure 6. Maximum blade surface temperature ver- sus the sputtering time. for the function ω∗i ≥ 0, where ω∗i = ωi√ ω2i + (grad ωi)2 , the expressions ∂ω(ω∗i )/∂x, ∂ω(ω ∗ i )/∂y, ∂ω(ω ∗ i )/∂z will determine the direction cosines of the interior normal to the corresponding coordinate axes [12]. Parameters of sputtering (such as value of the ions energy, rotation speed of the turntable and the blade) have been determined with respect to the temperature limit, so the maximum surface temperature should not exceed the temperature of the material phase transition. Additionally, the condition , in which the minimum value of the sputtered material layer should not be less than a predetermined value throughout surface of the blade, was applied. The graphs of the maximum temperature change in a checkpoint of the blade at a various ion beam energy is shown in Figure 6. Based on the simulation results, recommendations on the choice parameters of the sputtering of the blades have been made. Criterion of choice was the minimum time of the sputtering while respecting the temperature limitations. 476 vol. 56 no. 6/2016 Simulation of Surface Heating 4. Conclusions The use of R-functions can significantly reduce the amount of tracing computations. The above- mentioned effect is due to the analytical determination of points falling within the scope of the source. The point clouds produced by the 3D scanners can be used for this purpose. The numerical realization of the proposed method requires an accurate build numerical mesh. The best results in terms of accuracy of determination the scope of the source can be expected, when applying adaptive tunable meshes. The proposed approach can be applied to cases with an arbitrary number of heat sources considering that the law of body motion is known. It can also be used in mass transfer problems, such as simulation coating, painting or clearing bodies with complex geometric shapes. In case of integration of the R-functions in CAD system, the use of the proposed method would be simple enough. References [1] Komanduri, R., Hou, Z.B. Thermal modeling of the metal cutting process – Part II: temperature rise distribution due to frictional heat source at the tool–chip interface. International Journal of Mechanical Sciences 43:57-88, 2001. doi:10.1016/S0020-7403(99)00104-6 [2] Cline, H.E., Anthony, T.R. Heat treating and melting material with a scanning laser or electron beam. Journal of Applied Physics 48(9):3895-900, 1977. doi:10.1063/1.324261 [3] Rosenthal, D. The theory of moving sources of heat and its application to metal treatments. Transaction of the American Society of Mechanical Engineers 68:849-66, 1946. [4] Liu, S., Lannou, S., Wang, Q., Keer, L. Solutions for temperature rise in stationary/moving bodies caused by surface heating with surface convection. Journal of Heat Transfer 126:776-85, 2004. doi:10.1115/1.1795234 [5] Akbari, M., Sinton, D., Bahrami, M. Geometrical effects on the temperature distribution in a half-space due to a moving heat source. Journal of Heat Transfer 133:064502-1-10, 2011. doi:10.1115/1.4003155 [6] Zhou, H. Temperature rise induced by a rotating or dithering laser beam. Advanced Studies in Theoretical Physics 5(10):443-68, 2011. http://hdl.handle.net/10945/25571 [7] Jung, J.W., Lee, C.M. Cutting temperature and laser beam temperature effects on cutting tool deformation in laser-assisted machining. In Proceedings of the International MultiConference of Engineers and Computer Scientists (IMECS 2009) Vol. II, Hong Kong, March 18-20, 2009, pp. 1817-22. [8] Kay, D.S. Transparency, Refraction and Ray Tracing for Computer Synthesized Images. Master’s Thesis in Program of Computer Graphics, Cornell University, USA, 1979. [9] Whitted, T. An improved illumination model for shaded display. Communications of the ACM 23(6): 343-49, 1980. doi:10.1145/965103.807419 [10] Rvachev, V.L., Sheiko, T.I., Shapiro, V. The R-function method in boundary-value problems with geometric and physical symmetry. Journal of Mathematical Sciences 97(1):3888-99, 1999. doi:10.1007/BF02364929 [11] Voronyanskaya, M.E. Maksimenko-Sheiko, K.V., Sheiko, T.I. Mathematical modeling of heat conduction processes for structural elements of nuclear power plants by the method of R-functions. Journal of Mathematical Sciences 170(6):776-93, 2010. doi:10.1007/s10958-010-0120-x [12] Rvachev, V.L. On the analytical description of some geometric objects. Reports of Ukrainian Academy of Sciences 153(4):765-67, 1963. (in Russian) [13] Kersten, H., Deutch, H., Steffen, H., Kroesen, G.M.W., Hippler, R. The energy balance at substrate surfaces during plasma processing. Vacuum 63:385-431, 2001. doi:10.1016/S0042-207X(01)00350-5 477 http://dx.doi.org/10.1016/S0020-7403(99)00104-6 http://dx.doi.org/10.1063/1.324261 http://dx.doi.org/10.1115/1.1795234 http://dx.doi.org/10.1115/1.4003155 http://hdl.handle.net/10945/25571 http://dx.doi.org/10.1145/965103.807419 http://dx.doi.org/10.1007/BF02364929 http://dx.doi.org/10.1007/s10958-010-0120-x http://dx.doi.org/10.1016/S0042-207X(01)00350-5 Acta Polytechnica 56(6):472–477, 2016 1 Introduction 2 Mathematical formulation 3 Numerical simulation result and discussion 4 Conclusions References