The numerical modeling of the physical process of manufacturing parts using additive technologies is complex and needs to consider a variety of thermomechanical behavior. This is connected with the extensive use of the finite element computer simulation by means of specialized software packages that implement mathematical models of the processes. The algorithm of calculation of nonstationary temperature fields and stress-strain state of the structure during the process of 3D deposition of wire materials developed and implemented in ANSYS is considered in the paper. The verification of the developed numerical algorithm for solving three-dimensional problem of the production of metal products using arc 3D deposition of wire materials with the results of the experiment is carried out. The data obtained from calculations on the developed numerical model are in good agreement with the experiment.
- wire deposition
- thermomechanical behavior
- additive technologies
Additive technologies are a breakthrough solution of this century. At the same time, when we speak about additive technologies, we generally mean the manufacturing of products of small sizes and irregular shapes. Large-sized products are still characterized by the use of traditional casting and forging shops. If mass production, such as the motor vehicle industry, is generally satisfied with traditional solutions, then small-scale manufacturing, such as the aircraft industry, requires a more advanced approach. A manufacturing cycle of relatively simple and large-sized products, such as longerons, frame elements, etc., may in certain cases take 2 years at a cost of 1.2 million dollars. The situation at hand is one of those reasons which result in unreasonably long terms and high costs of creation and introduction of new products to the market.
The use of powder additive technologies sometimes gives rise to problems associated with low productivity of existing methods, high costs of equipment being used, limited types of materials, which is caused by the fact that powder systems melted by a powerful thermal source [1, 2, 3] are traditionally used as initial materials for additive formation of products. The formation of products from many aluminum alloys and active metal alloys, such as titanic and magnesium alloys, leads to increased porosity of materials of resulted products with considerably decreased mechanical properties [4, 5, 6, 7]. The productivity of formation of components from powder materials in traditional additive technologies is extremely low, which practically excludes any prospects of the use of these technologies for the purpose of manufacturing large-sized products.
Hybrid manufacturing technologies combine the best characteristics of additive formation of workpieces and those of subsequent mechanical removal of materials in the course of creation of metal products [8, 9]. This process can be implemented on one platform with the hybrid layer-by-layer application of wire materials and the processing by CNC machines and is an optimum solution for the manufacturing of large-sized components of molds of low and average complexity.
One of the first companies that is engaged in promoting wire material-related technologies is Sciaky (USA) [1, 2] that specializes in the development of electron-beam welding technologies and equipment. Additive manufacturing machines made by Sciaky produce components with the use of the layer-by-layer build-up welding method for materials in melts created by an electronic beam. This technology is known as EBDM (Electron Beam Direct Manufacturing). High performance levels (3–9 kg/h) demonstrated by the EBDM technology allow us to prepare components whose sizes are expressed in meters, which is impossible or extremely expensive when using any other additive technologies. This component formation principle is responsible for a low-quality surface of a synthesized component. However, the EBDM technology combined with traditional machining technologies allows us to obtain results with acceptable costs. Model materials used in this technology are additions (metal bars or wires), which is also an advantage as there are many available materials of this kind: nickel alloys, stainless and instrument steels, Co-Cr alloys and many others whose prices are considerably lower than those of these materials in their powdered condition . Today, the company does not make standard machines and is generally involved in producing basic Sciaky’s DM models with the sizes of a formation zone of 5700 × 1200 × 1200 mm, and all modifications are created according to customer requirements. Machines allow to consistently form up to 10 various components in automatic mode and during one vacuumization cycle of a working chamber. The price of one machine is more than $2.0 million.
The use of arc and plasma sources for the purpose of melting metal wire materials in the course of implementation of hybrid additive manufacturing technologies has been actively developed in recent years. In 2016, Norsk Titanium, a Norwegian start-up, attracted additional investments of 25 million dollars in order to certify materials resulted from the plasma layer-by-layer deposition with titanium wire and used to produce components for Boeing and Airbus aircrafts. It is also necessary to mention a company called WAAM (at Cranfield University) that is engaged in developing large-sized product formation technologies with the use of plasma technologies or consumable electrode surface welding processes with impulsive wire feeds and cold metal transfers (CMT) developed by Fronius. In the middle of 2016, Europe witnessed a 3-year LASSIM project with a budget of about 5 million euros that united 16 companies. The purpose of this project is to create a stand in order to implement several processes within a single space in the course of hybrid manufacturing of large-sized workpieces: additive manufacturing; multi-coordinate machining; layer-by-layer work hardening, measurement; non-destructive testing.
When modeling heat and mass transfer processes, it is necessary to consider that additive technologies are technologically closest to multi-layer deposition that also makes an initial material (filler wire) interact with a heat source, gradually builds up layers and superposes thermal cycles as and when new layers are added and transitional changes occur in the geometry. Temperature fields in a product to be processed are in most cases difficult to determine with the use of experimental methods. The most frequently used temperature measurement method consists in placing thermocouples directly next to the area of influence by a thermal source [10, 11, 12, 13]. Thermocouples can measure temperatures only in places where they are installed, and it is difficult to have a general picture of temperature distribution even when several thermocouples are used. Infrared thermography [14, 15, 16, 17] can measure only surface temperatures and cannot ensure the distribution of transient temperature processes in volume. The volume distribution of temperatures in a workpiece can be determined by mathematical modeling. However, the modeling of three-dimensional temperature fields is a rather demanding and complex problem due to accompanying difficult physical processes, their velocities, and many previously made calculations are connected with simplifications [18, 19, 20, 21].
In order to model heat and mass transfer processes, we usually use the following equation system: mass, momentum and energy conservation equations [22, 23, 24, 25]. Solutions of these equations allow us to obtain temperature fields in all projections of a sample to be formed, melt flow rates, cooling rates and crystallization parameters that specify structures and properties of components. At the same time, information on computational temperature fields obtained by using adequate models allows us to predict microstructures and properties of workpieces resulted from different additive manufacturing methods.
The process of manufacturing of components with the use of the additive manufacturing method is followed by complex thermo-mechanical phenomena resulting in the formation of technological residual stresses and possible contraction of components [26, 27, 28, 29, 30, 31, 32, 33]. The appearance of internal stresses in an object to be produced is connected with the essential spatiotemporal heterogeneous distribution of temperature and conversion fields.
A standard approach to the numerical solution of temperature deformations and residual stresses in additive manufacturing is to systematically and consistently analyze thermal conductivity in a transient mode and elastoplastic deformations . Further, a transient heat conduction problem is firstly solved in a numerical manner; then, a temperature field is imported to a mechanical model as “thermal loads” in order to calculate stresses and deformations.
In case of additive technologies with a rather small number of deposition passes, it is acceptable to thoroughly model each pass in the course of production of components [35, 36]. With this modeling method, the supply of heat brought by a thermal source is generally used as a volume thermal flow whose center moves along a deposition trajectory, thereby representing a moving source of heat. However, the additive formation of products usually has a large number of layers, which makes it unreasonable to model each separate pass in the course of creation of components. In order to make calculations more effective, we use a principle where successive melting steps and even layers are grouped together for subsequent simultaneous activation [37, 38]. This method provides that a stationary thermal flow is assigned to a discrete area for a period of time specified by users. It is obvious that the way of building-up of materials and of application of thermal loads in an individual manner is more correct, but requires higher computational efforts.
Therefore, the possibility of modeling of technological processes related to the layer-by-layer synthesis of products by multi-layer deposition of wire materials is a considerable reserve for the purposes of optimizing technological modes of production of components, developing control programs, minimizing defects and increasing manufacturing quality. At the same time, a lot of works are aimed at modeling selective laser melting or laser gas powder deposition welding processes. The modeling of arc methods of deposition welding of wire materials has specific features connected with a large volume of materials to be deposited and, as a result, with great possible deformations of products, as well as with specific aspects of the description of a thermal source .
General variable parameters for mathematical models of heat and mass transfer processes with deposition welding with wire materials are as follows:
distribution of density of an energy flow of a heat source;
initial temperature of a sample;
distribution of an additional volume source (in case of additional induction heating);
dependence of thermophysical characteristics of materials on temperature;
characteristics of phase transitions;
velocity of a heat source;
orientation and feed rate of wire and its diameter.
As a part of the mathematical models described below, the following assumptions have been accepted:
In the process in question, the strength of a plasma arc, the feed rates of a support material Vн and wire Vп are constant or time intervals of their change considerably exceed a typical time of relaxation of temperature, velocity and concentration fields.
The ambient temperature is constant.
An energy source is characterized by Gaussian distribution of a thermal flow onto surfaces of a support material and a bed to be deposited.
Support and wire materials have the same chemical composition.
Molten metal is considered to be an incompressible Newtonian liquid whose physical parameters (density, viscosity, thermal conductivity, etc.) do not depend on temperature.
When describing thermal effects in the course of melting and hardening of materials, effective thermal capacity within a quasi-equilibrium model is used.
In order to describe the influence of a two-phase zone on the motion of a melt, Darcy term, representing the damping force when fluid passes through a porous media dendrite structures, is introduced to the motion equation .
An impact on the motion of a melt of electric and magnetic fields generated by a plasma flow is not considered.
Filler wires are fed into the zone of influence of a plasma arc, thereby causing a mass inflow. A plasma arc causes melting and evaporation of a filler material and a base material. The surface of a melt can be somehow deformed under the influence of the arc pressure, of falling droplets (when melting a wire) or of a local increase in metal vapors pressure.
Technological process parameters affect the nature of a metal transfer to a melt. We can distinguish three typical modes:
Continuous metal transfer (Figure 1a). This mode is carried out at a low strength of a thermal source and is rather low-temperature. This mode slightly changes sizes and a form of cross section of a weld bed, as well as a structure and physical and mechanical properties of a metal to be deposited. This regime is optimal for the additive process. But their change may be caused only by disturbing factors (instable strength and position of a thermal source, deviation of a wire feed rate or displacement of a product, as well as influence of reheating zones). These factors are not considered at the current stage of works. A further decrease in energy to be delivered results in the fact that wires are melted only partially. This gives rise to formation defects and, therefore, this mode is not satisfactory.
A coarse-droplet transfer (Figure 1b) takes place when strength is increased above some critical value. The form of a weld bed and its cross dimensions range in length, metal splashes are present, and a crystallization process cannot be considered as a stationary one. This transfer mode may be acceptable for the consumable electrode welding technology. However, its use in most cases results in decreased quality when additive shape-forming processes are implemented.
A spray transfer (Figure 1c) is carried out in a mode with a higher energy of a thermal source. Metals being constantly fed experience a thermal influence sufficient for strong boiling. This mode is characterized by high spraying and uneven surfaces of weld beds to be formed.
2. Mathematical statements of problems (models) of a heat and mass transfer in the course of additive formation of products by melting of wire materials with a plasma arc
The geometry of a computational domain of the model in question is presented in a three-dimensional arrangement (Figure 2) and in the X-Z section (Figure 3). A surface source of a mass to be deposited onto a solid support material is fed to an interface at the velocity of Vп with some distribution in the X, Y plane (in case of a droplet transfer approximation is possible by Gaussian distribution). A thermal flow from an energy source is determined by Gaussian normal distribution. A minimum temperature in a source is supposed to be higher than the solidus temperature (TS). A mass source moves along a solid support material at the velocity of Vн (deposition welding velocity). There is the air above a solid support material.
This problem is supposed to be tackled by the shock-capturing method where motion and temperature distribution equations are solved within the whole domain presented in Figure 2. In order to describe the motion of a material interface, the level set method is used. The position of an interface and values of material parameters of media (density, viscosity, thermal conductivity, thermal capacity) are determined according to a value of a special remote function ϕ, to which a separate equation is assigned. The level set method helps to determine an interface position between a metal (a molten metal) and the air. The position of a boundary of a phase transition between a solid metal and a melt is established according to a position of an isotherm corresponding to a melting temperature.
Courses in metal and in air are described by free thermal convection equations for incompressible media in Boussinesq approximation :
where ρ—density, —melt flow velocity vector, μ—dynamic viscosity, P—pressure, —total volume force being as follows:
where g—gravitational acceleration, β—thermal-expansion coefficient, T—absolute temperature, Tref—initial temperature taken as the solidus temperature (TS).
The first equation member (2) represents thermo-gravitational convection, the second one does energy dissipation of in a two-phase zone (in the zone of a phase transition from a melt to a solid metal), according to the Kozeny-Carman equation where B—a small computational constant used to avoid division by zero, C—a constant reflecting the morphology of a two-phase zone (in these studies it is possible to use values around 104 … 106), fL—a function determining the position of a boundary between a liquid phase and a solid one in a metal (fusion zone boundary) :
where TL and TS—liquidus temperature and solidus temperature.
The third equation member (2) is a volume approximation of forces acting at the boundary of a phase interface.
The thermal energy distribution in the computational domain is described by means of the differential energy transfer equation:
where T—absolute temperature, —thermal diffusivity coefficient, u—melt flow rate, c—thermal capacity, ρ—density.
The approach in question provides that the location of a welding source of heat and a source of mass at the “metal-gas environment” interface are arbitrary functions of time. Besides, it is necessary to consider that a wire material is fed to the “metal-gas environment” interface in an already melted state, which also requires some part of energy . In an approximation that the temperature of a melted wire material to be fed is equal to the temperature of a melt at the interface, thermal capacity Q to be brought in the Eq. (4) is expressed as follows:
where Q0—maximum thermal power in a source, r—radius of a thermal spot from a source, —efficiency of a heat source, qh—thermal flow coefficient determined by a convective current, σсб—Stefan-Boltzmann constant, εч—metal blackness degree, ρL and ρS—densities of a liquid phase and a solid one of a metal, CpL and CpS—specific thermal capacities at a constant pressure of a liquid phase and a solid one of a metal, Lm—specific metal melting heat, NS(r)—function determining the distribution of a material to be delivered to a deposition zone, to be determined by deposition conditions (see Figure 3), δ(ϕ) and nz—delta function and normal projection to the boundary of the metal-air interface in the direction of the Z axis setting the distribution of heat along the boundary of the metal-air interface, will be determined below as a part of the description of the level set method.
The latent melting and crystallization heat was taken into consideration by introducing the effective thermal capacity:
where C0—thermal capacity depending on temperature, Hf—latent melting heat, Tmelt—melting temperature that is taken as an average one within a range from the solidus temperature to the liquidus one.
When in transition across the boundary of the interface of a liquid phase and a solid one in a metal, viscosity, first of all, suddenly changes from some final value in a melt to an actually infinite value in a solid metal . The value of viscosity in (2) was determined by means of function (4) as :
In calculations, will be final, but it will be such as , which will ensure the “freezing” of a liquid in a solid phase. Sudden changes in density, thermal diffusivity and thermal metal capacity when in transition across the boundary of a phase transition may be described in a similar way.
Surface forces are approximated as volume ones so that:
volume force is concentrated in a narrow transitional layer where density and viscosity gradients markedly differ from zero;
direction coincides with the direction of remote function gradient ϕ;
value is proportional to the gradient of a remote function and the curvature of a surface;
in a limiting case, when a transitional layer turns into a jump in typical parameters of a liquid, the consideration of this volume force leads to an ordinary dynamic condition on the surface of an interface.
Based on the mentioned assumptions, an expression for the volume force is written as:
where —marker function, —interface surface curvature,
—Dirac delta function, —melt surface tension coefficient.
The first addend in (8) describes capillary forces acting normally to the surface of the metal-air interface, the second addend describes thermo-capillary Marangoni forces acting tangentially to the surface where —temperature gradient component tangential to a curved surface. The third addend specifies the pressure force of plasma arc :
where Ia—arc current, kI—electrodynamic constant.
In the level set method, the position of an interface boundary is determined by a zero value of a level set or remote function ϕ. As for the remote function, we solve a transfer equation that will be as follows:
Density, viscosity or concentration values are restored according to the remote function. For example, in case of density we have:
where —media density, —Heaviside function set by the following expression:
For numerical implementation, expression (13) is written as
The thickness of a transitional layer between the phases is equal to . It is usually supposed that where is a distance between nodes of a spatial grid or a grid cell size around a boundary, α is an integer.
Far from a fusion zone, at all boundaries of the computational domain, except for an upper boundary, a media velocity is supposed to be equal to zero. At the upper boundary, pressure is supposed to be set. For the remote function, at all boundaries the normal derivative is supposed to be equal to zero. Thermal conditions at the boundaries of the computational domain are determined by external technological conditions, but for model calculations all boundaries may be considered to be thermally insulated.
3. Description of a thermal source when using a plasma arc at various polarity
When modeling heat and mass transfer processes with the use of a plasma arc as a thermal source, there is a factor concerning the description of a thermal source, especially when using a plasma arc with current reverse polarity. A heat transfer to a product, under the influence of a plasma direct arc, is carried out by two mechanisms: convection from a plasma spray (a plasma flow) and heat emission in discharge (electrode) spots. Work  establishes that at an identical current and under all other conditions being equal a heat input to a product, during the operation of a plasmatron with current reverse polarity, is 1.3… 1.6 times higher than with current direct polarity, which is explained by a higher arc voltage. Unlike a constricted arc of direct polarity, a constricted arc of reverse polarity is characterized by a thermal power distributed more uniformly along the surface of a product (Figure 4).
It should be noted that at plasma processing with a current of direct polarity values dpf and das are proportional, and it is impossible to control their size separately. The diameter of a plasma flow dpf is actually determined by the diameter of a nozzle dpf ≈ dn . The distribution of density of a full thermal flow at direct polarity deposition welding is described by Gaussian distribution.
With current reverse polarity, a plasma arc belongs to a type of arcs with non-stationary cathode spots traveling along the surface of a cathode. A travel width depends on the design of a plasmatron and the material of a product. One of distinctive features of non-stationary spots is their short-term existence and high current density (j∼105–106 A/m2), and local specific thermal flows in the area of short-term influence at a spot reach values (q∼106–107 W/cm2) (Figure 5). The time-averaged thermal influence of cathode spots can be approximated evenly by a thermal source distributed along the area restricted to dcs. The distribution of density of a thermal flow delivered by a plasma flow at reverse polarity deposition welding is described by Gaussian distribution. The diameter of a plasma flow dpf is supposed to be the diameter of a nozzle dcs ≈ dn . For the deposition of metal during the layer-wise formation, the following equipment which has been developed at Perm National Research Polytechnic University was used: the plasma welding unit BPS-350; universal plasmatron PM1-15.
A value dcs under the influence of an arc with current reverse polarity, can be actively controlled.
The size of a cathode region depends on parameters of a technological mode and requires determination. The travel of cathode spots leads to a known phenomenon of cathode cleaning under the influence of an electric arc of reverse polarity  on the surface of metals. It is generally believed  that the destruction and removal of an oxide film from a product surface in the area of influence of a reverse polarity arc result from attacking the surface of a metal by positive ions.
When a plasmatron operates with a reverse polarity current on a surface subjected to cathode cleaning, there is a visible mark  that allows us to visually estimate the area of traveling of cathode spots (Figure 6).
A phenomenon of cathode cleaning can be used to create a statistical dependence of the influence of technological parameters of reverse polarity plasma arcs on the diameter of a cathode region. At the first stage, it is possible to use an empirical model .
where Ia—arc current, Qpg—protective gas consumption (l/min), Qp—plasma-forming gas consumption (l/min), V—plasmatron displacement speed, m/h, H—plasmatron nozzle height above a product (mm).
Therefore, a plasma source at reverse polarity deposition welding is represented by a combination of a source of heat delivered by a plasma flow with Gaussian distribution, of diameter dpf, and an evenly distributed source of diameter dcs. An energy ratio between these sources is determined experimentally according to the method suggested in the work .
In order to analyze the applicability of selected approaches, the preliminary numerical implementation in the two-dimensional statement of a problem of melting of wire materials with concentrated and arc energy sources was carried out. A Comsol 4.4 application software package (Heat Transfer, Level Set and Laminar Flow modules) was used. The geometry of the computational domain is presented in Figure 7.
The sizes of the computational domain are 50 mm (length) and 10 mm (height). A process simulation area was covered with a two-dimensional grid included in the computational domain. The grid had an even pitch. The size of a cell was 0.3 mm.
A 304-L stainless steel was accepted as a model material. The stainless steel workpiece had a composition of 18.2% Cr, 8.16% Ni, 1.71% Mn, 0.02% C, 0.082% N, 0.47% Mo, 0.44% Si, 0.14% Co, 0.35% Cu, 0.0004% S, 0.03% P, and balance Fe. Many thermophysical characteristics of materials are functions of temperature [51, 52, 53, 54]. At this stage, these nonlinearities were not taken into consideration. Table 1 specifies accepted values that were used in calculations. Parameters of the deposition welding mode presented in Table 2. As boundary conditions for a model example all boundaries, except for a lower one, were taken as thermally insulated. At the lower boundary, a constant temperature of 273 K was maintained.
|Specific heat capacity||C||[J·kg−1·K−1]||500|
|Specific melting heat||Hf||kJ/kg||84|
|Source strength, W||Source displacement velocity, mm/s||Wire feed rate, mm/s|
Figure 8 shows the distribution of temperature in a sample to be deposited 2 s later after the beginning of a process. A melt area is marked by a line.
The deposition welding results obtained upon completion of a 4-s process are given in Figure 9. This model appropriately helps to determine distributions of temperatures, melt flow rates, pressures, components of heat flow densities, forms and sizes of weld beds to be deposited.
4. Development of the models intended to form fields of residual stresses and deformations in products resulted in the course of additive formation by melting wire materials with a plasma arc
The process of manufacturing of components with the use of the additive manufacturing method is followed by complex thermo-mechanical phenomena resulting in the formation of technological residual stresses and possible contraction of components [28, 30, 31, 32, 55, 56, 57, 58]. The appearance of internal stresses in an object to be produced is connected with the essential spatiotemporal heterogeneous distribution of temperature and conversion fields. The appearance of residual stresses is caused due to the fact that inelastic deformations  are not consistent, first of all, temperature shrinkage deformations at cooling, structural shrinkage due to the course of phase transformations (melt crystallization) that is notable for a deformation history of various material points because of heterogeneous temperatures, temperature gradients and temperature velocities.
In most cases, such measures as preliminary natural validations of technologies, selection of certain locations of components, selection of technological modes, are taken in order to eliminate and minimize any geometrical defects arising in the course of manufacturing of modern components, which is a serious obstacle to organizing comprehensive digital manufacturing processes.
An approach associated with creating and using mathematical and computer models of thermomechanical behaviors of components in the course of additive manufacturing allows us not to carry out natural experiments at a stage of design of technological process parameters and structural parameters and to predict qualitative and quantitative characteristics of stress conditions and contractions of future components, as well as study the effectiveness of possible measures aimed at decreasing residual stresses and deformations (heat treatment).
Calculation of residual stresses and warpage for the wire fusion process remains the most difficult aspect in numerical simulation. Often, the addition of material is modeled by adding and (or) activating new elements to previously placed ones. The growth of new elements adds additional rigidity to the existing structure and requires the gradual addition of more and more new elements to the solution area. There are three most commonly used methods of modeling material deposition—the so-called birth element (1), a quiet element (2) and hybrid activation (3) [59, 60]. In the method of the birth element, the elements for the material that has not yet been created are deactivated (and thus not included in the solution area), and then gradually regenerated and included in the solution area. In the method of sleeping elements, all elements are present in the calculation model from the very beginning and have artificial properties with very low rigidity. As the details grow, the properties of these elements gradually switch to real physical properties. Finally, the hybrid activation method combines the methods of emerging and sleeping elements, where only the current deposition layer is activated and set to a sleeping state, and all subsequent layers are deactivated .
For products with a relatively small number of surfacing passages, detailed modeling of each passage in the construction of a part is permissible [35, 36]. With this method of modeling, the heat input from the beam energy is usually applied as a volumetric heat flux whose center moves along the trajectory of surfacing, thus representing a moving heat source. Such a method with a moving heat source was used to model both Directed Energy Deposition (DED)  and Powder Bed Fusion (PBF)  additive production processes. Nevertheless, usually the product has a large number of layers, which makes it impractical to simulate each individual passage when creating a part. To ensure greater efficiency of calculations, the principle is used, in which the successive melting steps and even the layers are grouped together for subsequent simultaneous activation [37, 38].
We know studies aimed at modeling additive manufacturing processes, optimizing thermal cycles, estimating the influence of process parameters on changes in forms of finished products [59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80].
A large volume of research has been published for PBF from various materials, ranging from stainless steel [19, 63, 64, 65], carbon steels [78, 79, 80], nickel alloys [69, 70, 71, 72, 73, 74, 75] and titanium alloy Ti-6Al-4V [34, 69, 70, 71, 72, 75, 76]. By modeling residual stresses and deformations for additive technologies using wire, much less work has been published than for powder systems. The normalized maximum residual stress is presented in [19, 64] as a function of the power, scanning speed, and preheating of the substrate in the construction of a thin-walled structure. In  3D sequential temperature and elastic calculations were performed using the COMSOL and MATLAB packages. In , a 3D sequential temperature and elastic-plastic analysis was performed in the SYSWELD package. In works [66, 67] the ABAQUS is used for sequential temperature and elastic-plastic 3D analysis, deformation due to a phase transition is taken into account. In [62, 73], the COMET program and element activation technology are used, a thermo-elasto-viscoplastic material model is used. Additive technology Wire Arc Additive Manufacturing (WAAM) is suitable for manufacturing parts that require a large amount of surfacing [69, 72, 77]. In [69, 70, 71, 72] for WAAM technology in ABAQUS software, 3D analysis is performed, assuming the elastic-plastic behavior of the material.
To improve the efficiency of the calculations, Li et al.  proposed a method that displays the local residual stress field calculated at the mesoscale level for the rapid prediction of warping of a part. Another approach to efficient modeling warping is to use the inherent-strain method, developed by Yuan and Ueda to calculate the deformation in welding large-sized parts . The method directly applies a known proper (initial) deformation to calculate buckling and does not require a numerical solution of non-stationary temperature and elastoplastic problems in a step-by-step formulation. Finally, Mukherjee et al. [71, 75] constructed an analytical expression for the special deformation parameter for estimating the maximum residual strains as a function of the linear heat release, substrate stiffness, maximum temperature, the thermal expansion coefficient of the fusion alloy and the Fourier number, which expresses the ratio between the speed changes in the thermal conditions in the environment and the speed of the adjustment of the temperature field inside the system under consideration.
The numerical modeling of residual stresses and thermoshrinkable deformations at additive formation of products with the use of wire materials melted by a plasma arc is considered below.
In view of small deformations and negligibly small dissipative heat emission, it is possible to divide a boundary problem of non-stationary thermal conductivity and a boundary problem of thermomechanics with regard to a stressed-deformed state that are unrelated in such statement. Death and birth technologies (Elements Birth and Death in ANSYS) with regard to material parts (originally absent in a model and then added in the course of deposition) can be used in order to solve them. At the same time, an area occupied by an already finished product is considered as a computational one. The continuous building-up of a material is carried out discretely, at each substage of calculation corresponding to the “birth” of a next subarea from “dead” elements, a boundary problem of heat conductivity and thermomechanics is solved, and a result of solution of a previous substage serves as initial conditions for a subsequent one.
At substage , the statement of a boundary problem of non-stationary heat conductivity according to temperatures fields in area with boundary Sk, with accepted hypotheses taken into consideration, includes :
Thermal conductivity equation:
where , , —thermal capacity, thermal conductivity and density of an unevenly alloyed material, —specific strength of an external heat source.
where the first addend of a right part describes a convective heat transfer, and the second one—radiation (the Stefan-Boltzmann law); the third one—radiation from a plasmatron nozzle, —blackness coefficient, —Stefan-Boltzmann constant, —source radiation capacity, —distance between a surface point and a source, and —angle that is formed by a surface normal directed to a source, —heat transfer coefficient, —ambient temperature, —external single normal to the boundary of a body being cooled.
where —initial temperature distribution for substage , —temperature determined at the end of a previous one.
These relations consider that an area of study remains invariable throughout a substage. Here and are used to mark zones occupied by “live” and “dead” elements respectively. At the same time, thermophysical properties of the material in the area of “dead” elements are subject to degradation:
, , .
where —stress tensor;
Cauchy geometrical relations:
where —displacement vector, —total deformation tensor.
Displacement boundary conditions:
Stress boundary conditions:
where —parts of a boundary border with set displacements and loads respectively.
Thermomechanical properties of the material in the area of “dead” elements exclude physical nonlinearity and are perfectly elastic with degraded values:
where —the fourth-rank tensor of elastic constants of the material.
The general equation system of the deformable solid body mechanics boundary value problem includes constitutive relations as well. To describe the viscoelastoplastic behavior of the alloy under study, we used the Anand  model included to the list of ANSYS physical models. The model has the following form:
where —the fourth-rank tensor of elastic constants of the material; —the tensor of elastic deformations, —the tensor of total deformations; —the tensor of temperature deformations; —the tensor of viscous deformations; —the coefficient of temperature expansion of the material; —the onset temperature of temperature deformation; —the second-rank identity tensor. The viscous deformation rate is calculated according to the formula:
where —the stress deviator; —the average stress; —the equivalent stress (the following symbol “:” means a scalar product, the symbol “··”, set in earlier, means a tensorial product); —the equivalent viscous deformation rate connected to the other parameters by the following formula:
The model also includes the evolutionary equation:
The equation involves a possibility of work hardening and work softening. The designations are as follows:
; ; —the stress intensity; —the resistance to deformation (scalar function of hardening); —the saturation value of the hardening function; —the absolute gas constant; —the activation energy; —the absolute temperature; , , , , , , —the empirical coefficients.
Taking into account behavior characteristics of the elements activated according to the technology used at ANSYS, the formulas (25) are transposed to the following form:
where —the total deformation calculated up to the end of -th sub-step (the current status at the moment of activation of the element is its natural state).
As the setup of the mechanical problem (Eqs. (21)–(28)) imply, the ANSYS system of hypotheses does not include separation of the stress tensor into any components. The input of separate factors to the stress field generation is covered by the corresponding deformation types (see the formula (25)). At that, shrinkage due to phase transformation is covered by temperature deformation with the help of corresponding adjustment of the thermal coefficient of linear expansion within the solidus-liquidus interval.
The algorithm of calculations of temperature fields in numerical simulation of the plasma deposition process in the finite element ANSYS package stipulates performance of the following computational procedures:
To create the finite element model including volumes to be occupied by the product divided to separate horizontal layers as well as the platform—the basis with suitable thermal and physical properties.
To “kill” (the EKILL command) a part of the elements that are not involved in the real process of building up before its beginning.
In the cycle of layer building up beginning from the bottom layer:
To define conditions of convective heat exchange at the upper boundary of the layer according to the formula (19).
In the cycle of sub-steps of passing the deposition areas of the successive computational layer:
To remove a part of the layer at its bottom boundary under the area of previously defined conditions of convective heat exchange where available.
To activate (the EALIVE command) all elements of the -th area.
To heat a part of the area’s elements distributed around the volume using a thermal energy source (see the formula (18)) for a time of plasma arc impact. .
To remove the heat source and wait for the holding time .
To wait until the system cools down completely (partially). Solution of a non-stationary problem with known boundary conditions for a thermo-mechanical model, all of whose elements are active.
The impact time of the plasma arc , taking into account possible overlapping of plasmatron motion paths, is calculated according to the formula
where —the size of the spot, —the movement speed of the plasmatron, —the distance between adjacent tracks.
The holding time for the -th building up area is defined by the following expression:
where —the size of the area.
The specific capacity of heat emission on the surface of the heated spot of the deposited material is defined according to the following formula:
where —the capacity of the plasma arc, —the heating fraction of the capacity.
The specific volumetric capacity of the heat source used in the Eq. (18) has the following form:
where —the volume of the deposited area, —the thickness of the layer. Thickness of the deposited layer is defined on the basis of experimentally observed total height of the deposited part of the product and actual number of layers.
The discrete model of the computational domain is similar to the model of the heat conductivity problem above. Before calculation, the Solid279 heat element is replaced with Solid286, which uses displacements as the degrees of freedom, and thermomechanical properties are added. The algorithm is similar to that of Item 4.1.:
To “kill” (the EKILL command) a part of the elements that are not involved in the real process of building up before its beginning. To define boundary conditions in displacement (e.g. fixation along the cutting plane, vertical fixation in entrapment zones etc.).
In the cycle of layer building up beginning from the bottom layer:
In the cycle of sub-steps of passing the deposition areas of the successive computational layer:
To activate (the EALIVE command) all elements of the -th area.
To implement temperatures, previously calculated for that moment of time, to the units of the model for the impact time of the ray .
To implement temperatures, previously calculated for the end moment of the holding time, and hold for seconds.
To read the temperature field for the period of partially cooldown of the system, to calculate SSS.
To implement the ambient temperature, to release the entrapment, to calculate SSS.
After that we performed preliminary verification of the mathematical model and numerical algorithm for solution of the three-dimensional problem of metal products manufacturing using plasma deposition of wire materials. For the reference we took the results of the experiment for wire deposition on a metal base using arc surfacing in shielding gases with a tungsten (nonconsumable) electrode as described in the Article . Figure 10a shows the plate in question with the dimensions of 275 х 100 х 12 mm. We build up a metal rib with the length of 165 mm and the cross section showed in Figure 10b on the upper surface of the object.
The arc deposition process is performed in 2 steps. The first step involves preliminary heating of the base plate that is performed by torch back and forth movement with the speed of and the impacted area that is 1.5 times larger than the base area of the deposited rib. Power consumption of the torch is , diameter of the spot is . Subsequent deposition is performed in 10 layers along a continuous path (the process continues backward).
Figure 11 shows the finite element model of the problem. 20-unit isoparametric Solid 279 elements were used for solution.
The Inconel 718 nickel alloy was used as the material for calculation .
Thermal and physical properties of the material were taken from the Article . In addition to the above, according to the data, represented in the Article, the density and thermal conductivity were deemed as constants. Temperature dependence of the enthalpy is presented in Figure 12. Heat capacity of the material, included in the thermal conductivity equation, was calculated according to the formula . It should be noted that variations at the arc in the range of 1500–1600 К in Figure 12 arise from heat emission during solidification of the alloy. Heating fraction of the plasma arc capacity was set to 0.2 for calculation.
Identification of the chosen Anand model (25)–(28) for the material under study was performed according to the data of the stretch experiment with the defined speed at different temperatures represented in the Article .
The required constants of the model were defined in several steps by the downhill simplex method in the Matlab package. At that, relative disparity of computational and experimental values of voltage in Diagram (Figure 13) was minimized.
As the result, the following values of material constants were obtained by the authors of this publication:
= 2.0447e + 05 Pa, = 2.0864e + 03, = 0.3335, = 2.5008e + 04, = 0.4363, = 4.0352e + 06, = 3.2148e + 07 Pa, = 0.2273, = 0.4461.
Figure 13 shows that relative error of the calculation does not exceed 10%. The material was considered as isotropic. The Young’s modulus of elasticity and the Poisson ratio were equal to 153 GPa and 0.32 respectively and not dependent on temperature . The average value of CLTE at the temperature range up to the solidus temperature is 6.5 × 10−6 К−1. The lineal shrinkage on the interval between solidus and liquidus temperatures is 0.3% and included in the temperature dependence . Hardening and softening effects, accompanying the building up process, are included into the Anand model.
The main results of the calculations were given in . Verification was performed using the experimental results obtained in the paper [86, 87]. Figure 14 shows the layout of temperature sensors in the investigational studies mentioned above.
In this case sensors No. 2, 3, 5, 6 were placed on the upper surface of the base plate, and No. 1 and 4 were placed on its bottom surface. Initial direction of the torch movement is from right to left. Figures 15–17 show graphs of temperature evolution at the given points.
The figures show that the qualitatively retrieved data comply with the experiment. The worst quantitative match is observed at points on the down plane of the base close to the heating area (points 1 and 4, Figure 15). In addition, the curves have oscillations and the temperature is below the ambient temperature at some points. All in all, the experimental graphs are below the computational graphs. The oscillations can be explained by proximity of the heating area and, consequently, heavy temperature gradients substantially altering the result at a fairly coarse mesh. Higher heating level in the experiment may be also caused by failure to take account of the torch heat radiation expanded far beyond the computational size of the spot.
Figure 19 shows the results of solution of the deformable body mechanics problem and residual displacement at central cross sections of the plate. The figure shows that the residual contraction of the structure can be estimated with acceptable accuracy. It was predicted that the plate will be bowl-shaped with the height of the ribs about 1 mm. The minimum computational accuracy is at the cross section where the relative disparity of the height runs up to 20%.
Figure 20 represents the distribution pattern of residual characteristics in the structure.
Figure 20a shows that the longitudinal bend dominates. The transverse bend is almost 10 times lower. The maximum intensity of residual stresses is observed at the contact area of the deposited material and the base.
The stresses are calculated at every step of solution of the quasipermanent boundary deformable solid body problem (Eqs. (21)–(28)). Figure 20b represents distribution of residual stresses, i.e. the stresses observed in the structure at the end of the production process and complete cooldown to the ambient temperature.
The model of thermomechanical behavior of the sample, created by the additive fabrication method, is designed for evaluation of strain-stress state of the sample and its contraction as well as for definition of parameters of the production process ensuring low contraction and residual stress. The model can be combined with other mathematical and numerical models including the temperature and conversion model specifying distribution of the variable temperature field inside the sample.
The model may be included into a computer-aided expert system to forecast the results of additive fabrication of large-sized components.
We defined the mathematical problem (model) of heat and mass transfer in the process of additive fabrication of products by fusion of wire material using plasma (electric) arc and concentrated energy sources with asymmetrical wire feed. The model describes transient nonequilibrium adjoint heat and mass transfer processes in molten metal with free surface, includes differential equation of viscous fluid movement (Navier-Stokes equations) with convective terms in laminar current, takes into account the Marangoni effect on the surface of melt, and allows to calculate volume distribution of temperatures, melt flow speed, pressure, shape and dimensions of the molten pool, shape of the molten metal free surface, shape and dimensions of the deposited bead. We also defined potential uses of the models.
We suggested approximation of the wire feed process by mass inflow on the metal-gas edge. To describe movement of the molten-solid metal interface, one can use any of the two following methods: the Level-Set method or the Arbitrary Lagrangian-Eulerian method (ALE).
We described the heat source when using a plasma arc of various electrical polarity. In the process of deposition with reverse polarity, the plasma source is represented by a combination of the source of the heat, transferred by the plasma flow, with Gaussian distribution and a uniformly distributed source in the area of cathode spots impact. The energy relation between the sources is defined by experiment. We preselected the approximating formulas which describe influence of technological parameters of plasma arcs on distribution of heat flow density in the source.
We suggested the mathematical 3D-model of the additive product fabrication process by deposition of wire material.
We developed and implemented in the form of an APDL program in the ANSYS package the algorithm for estimation of unsteady temperature fields and strain-stress state of the structure in the process of its generation by arc 3D-deposition of wire materials.
The experimental information of other authors was used to identify the thermal-physical and thermomechanical properties of the material—the Inconel 718 nickel alloy.
Verification of the model showed required accuracy of the results.
The represented example of application of the model and its numerical implementation in the ANSYS package shows that the chosen approach is reasonable and allows to obtain acceptable numerical results properly describing real effects and processes in the objects created by the method of additive volume generation. More accurate verification procedure for the suggested mathematical model showing its accuracy and limits of applicability on the basis of experimental studies of residual stresses and thermoshrinkable deformations of model products manufacturing using the method of additive part fabrication by fusion of wire material with an electric arc and concentrated energy sources will be represented in future studies.
The work was supported by the Ministry of Education and Science of the Russian Federation (RFMEFI58317X0062) and MOST (No. 2017YFE0100100) under the BRICS project.