**World Journal of Mechanics**

Vol.04 No.09(2014), Article ID:50202,12 pages

10.4236/wjm.2014.49029

Efficient Simulation of Nonlinear Heat Transfer during Thermal Spraying of Complex Workpieces

Rolf Berthelsen^{1*}, Thomas Wiederkehr^{2}, Ralf Denzer^{1}, Andreas Menzel^{1,3}, Heinrich Müller^{2}

^{1}Institute of Mechanics, TU Dortmund, Leonhard-Euler-Straße 5, Dortmund, Germany

^{2}Lehrstuhl Informatik VII, TU Dortmund, Otto-Hahn-Straße 16, Dortmund, Germany

^{3}Division of Solid Mechanics, Lund University, Ole Römers väg, Lund, Sweden

Email: ^{*}rolf.berthelsen@udo.edu

Copyright © 2014 by authors and Scientific Research Publishing Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY).

Received 11 July 2014; revised 9 August 2014; accepted 1 September 2014

ABSTRACT

The quality of coatings, produced by thermal spraying processes, considerably decreases with the occurrence of higher residual stresses, which are especially pronounced for complex workpiece geometries. To understand the occurring effects and to aid in the planning of coating processes, simulations of the highly transient energy flux of the HVOF spray gun into the substrate are of great value. In this article, a software framework for the simulation of nonlinear heat transfer during (HVOF) thermal spraying is presented. One part of this framework employs an efficient GPU-based simulation algorithm to compute the time-dependent input boundary conditions for a spray gun that moves along a complex workpiece of arbitrary shape. The other part employs a finite-element model for a rigid heat conductor adhering to the computed boundary conditions. The model is derived from the fundamental equations of continuum thermodynamics where nonlinear temperature-depending heat conduction is assumed.

**Keywords:**

Nonlinear Heat Conduction, Continuum Thermodynamics, GPU Computing

1. Introduction

Thermal spraying is a cost-efficient coating technique for the production of wear-resistant surfaces consisting of various materials tailored to particular applications. For the coating of forming tools, e.g., hard material coatings of tungsten carbide (WC) and cobalt (Co) are used, [1] , because of their superior wear-resistance compared to chrome (Cr) and nickel (Ni) coatings, [2] . As a disadvantage, the High Velocity Oxygen Fuel (HVOF) thermal spraying process induces a large amount of energy into the heterogeneous coating and the substrate which leads to a complex transient thermo-mechanical problem, as illustrated in Figure 1. For an overview about HVOF thermal spraying, the reader is referred to, e.g., [3] [4] . While thermal-spraying of planar work pieces delivers coatings of quite satisfying quality, the quality considerably decreases with the rising degree of complexity of the work pieces’ geometry, for instance due to radii or curvatures. Experimental analyses suggest significant temperature exaltations in the substrate in dependence of the component geometry, whereby the coating quality decreases due to residual stresses during the cooling process. Thus, for the computer aided planning of thermal spray processes, the ability to predict the temperature development for a given workpiece is of high value.

This article presents a simulation framework for the computation of the temperature development for a given workpiece during the thermal spraying process. Therein the “outer” part of the framework handles the time stepping process, the spray gun movement and computes the time-dependent input boundary conditions for the “inner” simulation module, which evaluates a thermodynamically consistent transient and nonlinear heat conduction formulation for a rigid heat conductor based on the governing equations of continuum thermodynamics, as introduced by Coleman and Noll [5] .

For complex workpieces the setup of the input boundary conditions is a computationally expensive procedure by itself, since the (non-ambient) load distribution on the workpiece surface due to thermal radiation is not locally restricted and may affect a large fraction of the workpiece. The computation of the affected parts constitutes a computationally expensive visibility problem. For this, a novel approach is presented, which makes use of an elegant GPU-acceleration technique and is in fact an adaption of a method developed earlier in the context of mass deposition simulation [6] . This two-level simulation approach enables the highly detailed thermodynamic model, which is specifically adjusted to the HVOF spraying process, to be efficiently applied even to larger workpieces such as forming tools in the automotive industry or complex geometries like turbine blades. Another contribution of this work is the “inner” finite element simulation module including a multi-threaded C++ based implementation which uses the solver library Eigen.

The following two sections describe the fundamental equations of continuum thermodynamics for a rigid heat conductor and the resulting finite element discretisation. Section 4 outlines the functionality of the robot guided spray simulation, [6] , and provides important details regarding the GPU-accelerated computation of the input boundary conditions. Section 5 presents a demonstration of the software tool for the simulation of a real work piece. The material parameters of the underlying constitutive relations―the heat capacity as well as the heat conduction coefficient―are represented by suitable functions which are fitted to experimental data.

Figure 1. Photograph of the HVOF spraying process, kindly provided by LWT, TU Dortmund.

2. Continuum Thermodynamics

The governing well-established equations of continuum thermodynamics for a rigid heat conductor are illustrated in the following section. The fundamental balance equation for a rigid heat conductor is the balance of energy, here in local form,

(1)

wherein denotes the rate of the internal energy, denotes the heat flux vector and denotes the external heat supply, while represents the material density, [5] . The balance of energy equation represents the first law of thermodynamics. By assuming the entropy flux to be proportional to the heat flux, as proposed by Coleman and Noll [5] , the general entropy inequality results in the so called Clausius-Duhem inequality

(2)

In the above inequality, is the rate of the entropy and is the absolute temperature. The Clausius- Duhem inequality represents the second law of thermodynamics. A Legendre transform of the internal energy results in the free energy density so that and subsequently. Following Liu [7] and references cited therein, the only remaining part of inequality (2) is the thermal dissi- pation, also referred to as Fourier’s inequality which is fulfilled if the heat flux points from a hot spot to a cold spot. To maintain generality, the specific heat is introduced temperature-dependent, i.e.. The energy balance (1) is now restated in the strong form of the temperature field equation valid for a general rigid heat conductor,

(3)

This equation is to be solved within a general thermodynamical initial boundary value problem. For this pur- pose, the boundary of the body is decomposed into three disjoint parts, i.e., with, and. On Dirichlet boundary conditions are prescribed for the temperature, whereas Neumann and Robin boundary conditions are prescribed for the heat flux on respectively:

(4)

Fourier’s law of heat conduction is applied throughout this work. For a review on non-Fourier heat conduction the reader is referred to other recent works, e.g., the works of Atefi and Talaee [8] [9] .

3. Finite Element Discretisation

The different representations of the energy balance equation derived above are given in strong form. To calculate a solution for the desired field of the temperature by means of the finite element method in the context of inhomogeneous initial boundary value problems, the temperature-based balance of energy has to be reformulated in weak form. Therefore, we transfer Equation (3) to a residual form,

(5)

where is a scalar-valued test function, see [10] for detailed background information. This relation also holds under integration over the volume of a body,

(6)

The divergence of the heat flux in Equation (6) can be reformulated by the application of Gauß’s theorem and integration by parts,

(7)

Here, denotes the outward unit vector on the boundary of the body. By interpreting the test function as the virtual temperature, which is 0 at Dirichlet boundaries, Equation (7) represents the virtual temperature problem and can be written in terms of dynamic, volume, internal and surface terms

(8)

which, in turn, are defined by

(9)

Equation (8) is the weak form of the initial boundary value problem of a rigid heat conductor which now has to be discretised in time and space, following the procedure outlined, e.g., in [11] or Kuhl et al. [12] . For time discretisation, a differential-quotient-based Backward Euler integration scheme is applied [13] ,

(10)

where symbolises an arbitrary quantity of interest and denotes the time increment. Subscript denotes a quantity at the actual time step, whereas subscript denotes a quantity at the previous time step. Application of Equation (10) to the temperature results in

(11)

In addition to the Backward Euler method applied in this work, recent works which address thermal problems use the Crank-Nicolson method [14] or energy-momentum consistent schemes in the context of thermo- elastodynamics [15] . At time, Equation (11) leads to the relation

(12)

for the unknown temperature. In view of the discretisation in space, it has to be taken into account that Robin boundary conditions require additional effort, since they represent one type of temperature dependent loads. For a review on the implementation and algorithmic treatment of deformation dependent loads, the reader is referred to the textbook by Bonet and Wood [16] and references cited therein. For that purpose, on the one hand the body is approximated by a finite number of volume elements and, on the other hand, the Robin boundary of the body is approximated by surface elements, i.e.

(13)

Following the spirit of the isoparametric concept, the geometry of the body, in terms of position vectors, as well as the temperature of the body and the virtual temperature are interpolated element-wise by shape functions respectively and discrete node point positions of respectively element nodes, i.e.

(14)

Hence, the gradients, and are approximated as

(15)

The insertion of the approximations given by Equations (14) and (15) into Equation (12) leads to the fully discretised balance of energy in terms of the unknown temperature

(16)

where the discrete inertia, volume, internal and surface fluxes are expressed by

(17)

In the above equations, represents the assembly operator over all volume element contributions and all surface element contributions at respectively element nodes to the global node points. Note, that and are assumed to be temperature- independent throughout this work. With the definitions in Equation (17), Equation (16) represents the discretised nonlinear temperature field equation. The solution is performed by an incremental Newton-Raphson scheme, see [17] for detailed background information in the context of the finite element method. For that purpose, the Jacobian of the residual with respect to the temperature has to be determined. Therefore, the iterative scheme can be expressed as

(18)

wherein is approximated by the linear term of a Taylor series expansion,

(19)

The derivative is assembled to the global tangent operator

(20)

The first term represents the dynamic contribution and therefore characterises the time-dependency of the problem, whereas the second term corresponds to the internal heat flux where the constitutive relation for the heat flux still has to be defined. The third term reflects the surface element contribution by Robin boundary conditions.

4. Implementation

Two aspects are vital for the computer aided planning of thermal spray processes: the prediction of the coating thickness distribution on the surface of a workpiece and the prediction of the temperatures reached in the workpiece during the process. For the simulation of the coating thickness, an efficient GPU-accelerated, C++ based simulation was developed in [6] . In the present work, this framework is modified so as to compute a time dependent temperature distribution on the surface of the discretised workpiece mesh which serves as an input for the finite element framework developed in the previous section. In order to model the energy input of the unloaded flame of the spray gun into the workpiece, the surface heat flux is assumed to be applied by convection and radiation. An unloaded flame means that the problem at hand is restricted to heat transfer only without accounting for mass transport phenomena. One simple approach to mathematically cover the convective heat contribution is the introduction of a so-called film condition

(21)

wherein denotes the convective heat transfer coefficient which is assumed to be a constant, and where denotes the environmental temperature which either represents the HVOF spray gun or the room temperature. It should be noted that represents a rather complex fluid mechanical process, not captured by the finite element approach as this work proceeds. The heat radiation is modelled by

(22)

where is the emissivity and represents the Stefan-Boltzmann constant, cf. Comini et al. [18] . Both, and, are Robin boundary conditions in manner of Equation (4), i.e.. For implementation, the heat capacity is specified by a cubic polynomial

(23)

To guarantee thermal stability as discussed above, it has to be ensured that holds for arbitrary (positive) values of the absolute temperature. Furthermore, the heat conduction is modelled by Fourier’s law

(24)

where has to be positive semi-definite for any, see e.g. [7] , and is thus represented by

(25)

where characterises the second order identity tensor and where is the heat conduction coefficient. The assumption of being proportional to leads to an isotropic heat conduction model. The nonlinear heat conduction model specified via Equations (21)-(25) is implemented in C++ as part of the already mentioned simulation framework. The remaining part of this section describes the computation of the film condition parameters for all boundary triangles of the tetrahedron volume mesh representing the workpiece.

Apart from the triangulated workpiece mesh, the simulation approach takes a robot guided spray gun movement path as an input, which is represented as a sequence of 3-tuples for discrete robot path positions. Therein is the spray gun position in global coordinates and is a quaternion describing the orientation of the spray gun at time. For detailed background information on quaternions the reader is referred to Argyris [19] , Betsch et al. [20] , Altmann [21] and the references cited therein.

During the simulation a virtual spray gun is moved along the given path in discrete time steps, where positions are interpolated linearly and where quaternion slerp interpolation is used for the orientations [22] . For every gun position the following computations, where steps 1-3 are entirely GPU-based and implemented in OpenGL and the OpenGL Shading Language (GLSL), are performed:

1) Compute the subset of nodes of the workpiece mesh that receive a heat flux from the current spray gun position.

This is defined to be the subset of nodes, the coordinates of which are inside the spray plume, approximated by a (spray) cone which expands from the gun center position, see Figure 2. Only those nodes visible from the spray gun center are considered; nodes on the back side of the workpiece or nodes shadowed by other parts of the workpiece do not receive a high temperature load.

This visibility computation for arbitrary meshes with potentially hundreds of thousands of nodes and faces can be efficiently implemented by exploiting the computational power of modern graphics cards, tailored by design to perform visibility computations on large triangle meshes in real-time. For the rendering of virtual scenes consisting of triangle meshes, the graphics card projects the three dimensional nodes of the mesh in a perspective correct manner onto a two-dimensional image plane. In the first visibility computation step, projected nodes outside the rectangular screen area are discarded, limiting the remaining nodes to the ones inside a pyramidal viewing frustum as depicted in Figure 2. In an additional step, the graphics hardware uses the depth-buffer algorithm [23] to remove any nodes that are shadowed by the geometry for the particular viewpoint considered.

Figure 2. Visibility computation to determine the nodes of the workpiece to be coated and to receive a heat flux from the gun (white nodes). Nodes outside the spray cone or nodes that are shadowed from the gun by parts of the workpiece are not coated.

For the current spray gun position, first, the opening angle of the viewing frustum of the virtual camera is adjusted to match the desired opening angle of the spray cone. Subsequently, the workpiece mesh is rendered and the desired nodes that receive a heat flux from the gun are determined. The case that the viewing frustum is rectangular, but the spray cone should be circular is corrected in the next step.

2) For every node in the subset determine the thermal load in order to set up Robin type boundary conditions.

The nodes are projected onto a plane orthogonal to the central axis of the spray cone at the stand off distance from the gun center. The thermal load for every node is then computed based on its distance from the cone center axis. Due to the lack of precise measurement data, the load function is modelled as a rotationally symmetric Gaussian with standard deviation and amplitude which is offset by in order to cover a larger circular area of magnitude:

(26)

with. The threshold is used to limit the function to a circular radius thereby transforming the rectangular viewing frustum into the desired circular spray plume approximation. For the simulations presented in this article an HVOF spraying process is assumed and the values, , , and are used. The shape of the load function is shown in Figure 3.

3) The remaining nodes are set to receive an ambient thermal load in order to set up Robin type boundary conditions.

4) The thermodynamic model is evaluated.

Figure 3. Shape of the radially symmetric load function used in the benchmark simulations; Equation (26).

For each Newton iteration, the C++ implementation of the thermodynamic model assembles the global iterative residual vector and the global tangent matrix according to Equations (19) and (20), and then employs a conjugate gradient iterative solver with Jacobi (diagonal) preconditioning (pCG) to solve the system.

5) The simulation time index is increased by.

This scheme is iterated until the end of the robot guided path is reached. A key aspect for the efficient computation of the boundary conditions is the similarity of the geometric relationships of the spray plume expanding from the gun nozzle towards the workpiece on the one hand, and the geometry-projection process towards a virtual camera center in the visibility computation on the other. Due to this similarity, the first three steps can easily be performed by exploiting GPU acceleration techniques and are implemented in OpenGL and the OpenGL Shading Language (GLSL) within our framework.

5. Examples

The implemented framework is now applied to carry out simulations of a realistic deep drawing tool in order to demonstrate the capabilities of the novel software. The starting temperature of the material is considered to

coincide with the room temperature whereas the temperature of the spray gun is set to. The steel material is St 35.8 (1.0305), the physical properties of which have been investigated by Richter [24] [25] up to a temperature of 873 K. The parameters for Equation (23), the polynomial of the heat capacity, are taken from [24] [25] and summarised in Table 1. Note, that the polynomial is parameterised in ˚C not in K. For the sake of completeness, the polynomial is depicted with a K-scale in Figure 4(a). It can be seen that the parameter set is not valid for temperatures beyond the interval. Since is guaranteed within the temperature range of interest, , the implemented framework remains thermally stable even when the physical valid temperature range is excelled. This turns out to be helpful for testing the algorithmic framework.

Furthermore, the parameters for the heat conduction coefficient, given by Equation (25), are also summarised in Table 1. Because the parameters are the result of fitting a polynomial of Richter, [25] , both functions are depicted in Figure 4(b). The solid line represents the fitted function and, in contrast, the dashed line represents the polynomial taken from Richter, [25] . Obviously, the fit is very good within the physical valid temperature range. Beyond this range, the fitted curve asymptotically tends to a positive limit while the Richter polynomial has a root at. As an advantage, the fitted function guarantees for any θ, thereby avoiding numerical problems while testing the algorithmic framework. To complete the set of parameters, values for the convective heat transfer coefficients and as well as for the emissivity ε and the Stefan-Boltzmann constant σ are summarised in Table 2. We assume the surface surrounded by air to be exposed to natural convection, as the surface loaded by the HVOF spray gun is assumed to be charged with forced convection.

The integrated coating and thermodynamic simulation is employed to predict the workpiece temperature for the manufacturing of thermal sprayed sheet metal forming tools. In this regard, the heat transport into a deep drawing tool of a moving and unloaded HVOF thermal spray gun is simulated using a complex robot tool path at three different movement speeds. The fastest simulation run completes the tool path within 40 s, the medium speed simulation in 80 s and the slowest run within 160 s. The tool path of the spray gun movement is depicted

Table 1. Parameters for introduced in Equation (23) and introduced in Equation (25). The values are valid within; see [25] for further details.

Table 2. Values for the convective heat transfer coefficients and, the emissivity and the Stefan-Boltzmann constant.

(a) (b)

Figure 4. (a) Richter polynomial of the heat conduction; (b) Fitted curve and the Richter polynomial of the heat conduction coefficient; see [25] for further details.

in Figure 5. The colors of the robot path represent the temperature load resulting from the respective gun position. For the robot path, the blue regions on the right and left parts have not been processed yet, since the gun started at the lower end of the path proceeding upwards along the circular opening of the workpiece onto the plain surface on top. For the path already processed, the coloring represents the surface load temperatures induced by the gun, where grey areas represent zero thermal load. In other words, these grey parts represent positions of the spray gun where the spray cone does not intersect with the workpiece. The surface colors represent the surface temperature at different time steps for the medium gun speed simulation run.

The simulation keeps track of the volume-averaged temperature and the peak workpiece temperature both with respect to time to identify critical sections of the path, which may cause an overheating of the workpiece. Figure 6 shows the simulated temperatures averaged over the entire workpiece volume and the peak temperature occurring in the workpiece―both plotted versus time―for the three gun speed simulation runs. The spraying process can be divided into four sections―compare Figures 5(b)-(d): First, the spray gun aims at the front radius working its way up to the plain surface, which is then coated in circular arc movements from left to right and back. During this process, peak temperatures are reached when the gun approaches the border of the workpiece and turns around for the next arc in opposite direction. The fluctuations of the peak temperature curves are due to the gun leaving and entering the workpiece at the borders. The third and fourth distinguishable sections are the right and left border areas of the workpiece including the screw holes for mounting the part. At the beginning of these sections, the peak and average temperatures increase considerably, because the large screw holes provide a significantly larger surface area that is heated by the spray gun. Afterwards the temperatures drop again when the gun gradually leaves the workpiece area. It can be seen that, for the considered robot tool path, the slowest simulated gun movement speed leads to the highest mean temperature of the workpiece as well as to the highest surface peak temperatures. The heat output to the environment through the surface is lower than the heat input by the spray gun.

(a) (b)(c) (d)293 K 443 K

Figure 5. Simulated deep drawing workpiece and the employed robot movement path. The surface colors represent the surface temperatures and the path colors the temperature load after (a) 6 s; (b) 20 s; (c) 30 s; and (d) 40 s of the medium gun speed simulation run. Blue path segments have not been processed yet.

Figure 6. Simulated development of peak and mean workpiece temperature of the workpiece shown in Figure 5. The temperatures are depicted for the three gun speed simulation runs.

6. Summary

This work presents a novel integrated simulation approach for the computation of the workpiece temperature development during thermal spraying which is suited for large and arbitrarily complex shaped components. As its core, a thermodynamically consistent nonlinear heat transfer model based on the fundamental equations of continuum thermodynamics is employed. The equations obtained are transformed into their weak forms that represent the underlying equations of the presented non-linear finite-element-framework, implemented in C++ and the Eigen library. The material parameters are adjusted to experimental data by fitting suitable functions, and the parameters of the Robin boundary conditions are chosen in order to model realistic heat transfer. The presented simulations of a complex shaped deep drawing tool demonstrate a possible field of application for the novel software tool.

Future work will include the consideration of the mass flow of the HVOF gun in order to model the coating process of the substrate with hot particles. Furthermore, the developed framework shall be used to extend the optimisation software described in [26] in order to be able to optimise temperature quantities. Another aspect of future research is the extension towards a thermo-mechanically coupled framework, [27] , which will enable the analysis of residual stresses and damage due to thermal spraying as discussed in, e.g., [28] [29] .

Acknowledgements

This research was funded by the German Research Foundation (DFG) as part of the collaborative research center 708 “3D-Surface Engineering of Tools for Sheet Metal Forming - Manufacturing, Modeling, Machining” within the projects B1 and B6.

References

- Trompeter, M., Franzen, V., Witulski, J. and Tekkaya, A.E. (2009) Thermisch beschichtete Werkzeuge für die Blechumformung. Der Schnitt- & Stanzwerkzeugbau, 5, 6-12.
- Sahraoui, T., Fenineche, N.-E., Montavon, G. and Coddet, C. (2003) Structure and Wear Behaviour of HVOF Sprayed Cr3C2-NiCr and WC-Co Coatings. Materials & Design, 24, 309-313. http://dx.doi.org/10.1016/S0261-3069(03)00059-1
- Thorpe, M. and H.J., R. (1992) A Pragmatic Analysis and Comparison of HVOF Processes. Journal of Thermal Spray Technology, 1, 161-170. http://dx.doi.org/10.1007/BF02659017
- Fauchais, P., Vardelle, A. and Dussoubs, B. (2001) Quo Vadis Thermal Spraying? Journal of Thermal Spray Technology, 10, 44-66. http://dx.doi.org/10.1361/105996301770349510
- Coleman, B.D. and Noll, W. (1963) The Thermodynamics of Elastic Materials with Heat Conduction and Viscosity. Archive for Rational Mechanics and Analysis, 13, 167-178. http://dx.doi.org/10.1361/105996301770349510
- Wiederkehr, T., Müller, H., Krebs, B. and Abdulgader, M. (2009) A Deposition Model for Wire Arc Spraying and Its Computationally Efficient Simulation. Proceedings of the International Thermal Spray Conference (ITSC), Las Vegas, 492-498.
- Liu, I.-S. (2002) Continuum Mechanics. Springer, Berlin.
- Atefi, G. and Talaee, M.R. (2011) Non-Fourier Temperature Field in a Solid Homogeneous Finite Hollow Cylinder. Archive of Applied Mechanics, 81, 569-583. http://dx.doi.org/10.1007/s00419-010-0436-5
- Talaee, M.R. and Atefi, G. (2011) Non-Fourier Heat Conduction in a Finite Hollow Cylinder with Periodic Surface Heat Flux. Archive of Applied Mechanics, 81, 1793-1806. http://dx.doi.org/10.1007/s00419-011-0518-z
- Lewis, R.W., Nithiarasu, P. and Seetharamu, K.N. (2004) Fundamentals of the Finite Element Method for Heat and Fluid Flow. Wiley, Hoboken and N.J. http://dx.doi.org/10.1002/0470014164
- Bergheau, J.M. and Fortunier, R. (2008) Finite Element Simulation of Heat Transfer. Wiley-ISTE, London. http://dx.doi.org/10.1002/9780470611418
- Kuhl, E., Denzer, R., Barth, F.J. and Steinmann, P. (2004) Application of the Material Force Method to Thermo- Hyperelasticity. Computer Methods in Applied Mechanics and Engineering, 193, 3303-3325. http://dx.doi.org/10.1016/j.cma.2003.09.021
- Vujicic, M.R. (2006) Finite Element Solution of Transient Heat Conduction Using Iterative Solvers. Engineering Computations, 23, 408-431. http://dx.doi.org/10.1108/02644400610661172
- Palani, G. and Kim, K.Y. (2010) Numerical Study on a Vertical Plate with Variable Viscosity and Thermal Conductivity. Archive of Applied Mechanics, 80, 711-725. http://dx.doi.org/10.1007/s00419-009-0336-8
- Gross, M. and Betsch, P. (2011) Galerkin-Based Energy-Momentum Consistent Time-Stepping Algorithms for Classical Nonlinear Thermo-Elastodynamics. Mathematics and Computers in Simulation, 82, 718-770. http://dx.doi.org/10.1016/j.matcom.2011.10.009
- Bonet, J. and Wood, R. (2008) Nonlinear Continuum Mechanics for Finite Element Analysis. 2nd Edition, Cambridge University Press, Cambridge.
- De Borst, R., Crisfield, M., Remmers, R. and Verhoosel, C. (2012) Nonlinear Finite Element Analysis of Solids and Structures (Wiley Series in Computational Mechanics). Wiley, Hoboken. http://dx.doi.org/10.1002/9781118375938
- Comini, G., Del Guidice, S., Lewis, R.W. and Zienkiewicz, O.C. (1974) Finite Element Solution of Non-Linear Heat Conduction Problems with Special Reference to Phase Change. International Journal for Numerical Methods in Engineering, 8, 613-624. http://dx.doi.org/10.1002/nme.1620080314
- Argyris, J. (1982) An Excursion into Large Rotations. Computer Methods in Applied Mechanics and Engineering, 32, 85-155. http://dx.doi.org/10.1016/0045-7825(82)90069-X
- Betsch, P., Menzel, A. and Stein, E. (1998) On the Parametrization of Finite Rotations in Computational Mechanics: A Classification of Concepts with Application to Smooth Shells. Computer Methods in Applied Mechanics and Engineering, 155, 273-305. http://dx.doi.org/10.1016/S0045-7825(97)00158-8
- Altmann, S.L. (2005) Rotations, Quaternions, and Double Groups. Dover Books on Mathematics. Dover Publications, Mineola.
- Shoemake, K. (1985) Animating Rotation with Quaternion Curves. SIGGRAPH Computer Graphics, 19, 245-254. http://dx.doi.org/10.1145/325165.325242
- Catmull, E.E. (1974) A Subdivision Algorithm for Computer Display of Curved Surfaces. Ph.D. Dissertation, The University of Utah, Salt Lake City.
- Richter, F. (1973) Die wichtigsten physikalischen Eigenschaften von 52 Eisenwerkstoffen: Mitteilung aus dem Forschungsinstitut der Mannesmann AG. Verl. Stahleisen, Düsseldorf.
- Richter, F. (1983) Physikalische Eigenschaften von Stählen und ihre Temperaturabhängigkeit: Polynome und graphische Darstellungen. Stahleisen, Düsseldorf.
- Kout, A. and Müller, H. (2009) Parameter Optimization for Spray Coating. Advances in Engineering Software, 40, 1078-1086. http://dx.doi.org/10.1016/j.advengsoft.2009.03.001
- Baek, S. and Srinivasa, A.R. (2003) Thermomechanical Constraints and Constitutive Formulations in Thermoelasticity. Mathematical Problems in Engineering, 2003, 153-171. http://dx.doi.org/10.1155/S1024123X03212011
- Boehmer, J., Funk, G., Jordan, M. and Fett, F. (1998) Strategies for Coupled Analysis of Thermal Strain History during Continuous Solidification Processes. Advances in Engineering Software, 29, 679-697. http://dx.doi.org/10.1016/S0965-9978(98)00033-7
- Fagerstrom, M. and Larsson, R. (2008) A Thermo-Mechanical Cohesive Zone Formulation for Ductile Fracture. Journal of the Mechanics and Physics of Solids, 56, 3037-3058. http://dx.doi.org/10.1016/j.jmps.2008.06.002

NOTES

^{*}Corresponding author.