Prediction of the Weld Pool Stability by Material Flow Behavior of the Perforated Weld Pool

This article presents the application of a computational fluid dynamics (CFD) finite volume method (FVM) model for a thermo-mechanical coupling simulation of the weld pool used in variable polarity plasma arc welding (VPPAW). Based on the mechanism of the additional pressure produced through self-magnetic arc compression and the jet generated from mechanical plasma arc compression, and considering the influence of arc height and keyhole secondary compression on arc pressure, a three-dimensional transient model of variable polarity plasma arc (VPPA) arc pressure was established. The material flow behaviors of the perforated weld pools were studied. The results show that three kinds of flow behavior existed in the perforation weld pools and it is feasible to predict the weld pool stability by the material flow behaviors of the perforated weld pools. The weld pools can exist stably if the material flow in the bottom of the perforated weld pools can form confluences with moderate flow velocities of 0.45 m/s, 0.55 m/s and 0.60 m/s. The weld pools were cut when the material flowed downward and outward with the maximum velocity of 0.70 m/s, 0.80 m/s. When the maximum material flow velocity was 0.40 m/s, the weld pool collapsed downward under the action of larger gravity. The thermo-mechanical coupling model was verified by the comparison of the simulation and experimental results.


Introduction
Aluminum alloy 2219 (AA2219) has the characteristics of high specific strength, good mechanical properties at low and high temperatures, high fracture toughness and excellent stress corrosion resistance. It is one of the ideal materials for manufacturing large-scale space tanks [1]. Variable polarity plasma arc welding (VPPAW) has great application potential in the welding of aluminum alloys with thicknesses of 4-20 mm [2], which can be used to realize the double-sided formation of an aluminum alloy plate through single-sided welding without a forming groove, also known as the "zero defect" welding method [3][4][5]. However, the application of VPPAW in the manufacturing industry is restricted by its narrow process range and intense process sensitivity.
Han et al. [6,7] investigated the process characteristics and weld stability of aluminum alloy in VPPAW. However, most of the welding parameters were determined empirically, requiring numerous experiments for verification. With the development of computer science, numerical simulation can be used to understand the mechanisms and physical interactions of various welding processes [8]. Using stainless steel as the welding material, Wu et al. [9][10][11][12][13][14] developed a thermo-mechanical model to gain a better understanding of energy propagation and calculate the dynamic keyhole evolution and the fluid flow in the weld pool. Wu, Shinichi Tashino et al. [15,16] investigated the weld pool convection and keyhole formation mechanism in the keyhole plasma arc welding (PAW). The experimental methods of the above numerical studies were PAW. For VPPAW, Yang et al. studied the change of keyhole size and weld pool material flow with the change of heat input in the process of vertical welding on the premise of taking the average of the positive and negative polarity current [17,18]. Chen et al. [19][20][21] mainly studied the physical properties of 5A06 and A5052P aluminum alloy during VPPAW.
Therefore, numerical simulation analysis regarding to the material flow behavior of AA2219 weld pool has seldom been carried out. However, AA2219 has a high thermal conductivity, so it dissipates heat quickly in the welding process, and the keyhole of AA2219 is more hardly formed than that of steel or 5-Series aluminum alloy in VPPAW [15]. Therefore, it is important to study the material flow behavior of AA2219 during VPPAW.
VPPAW essentially belongs to perforated welding, and thus, the keyhole evolution is extremely important. Generally, the keyhole is formed in several seconds, which is proposed to be the dominant driving force of arc pressure [15]. The arc pressure model established at present mainly considers the change of current and the arc pressure distribution along the thickness direction of the workpiece [9,10,19,22]. For VPPAW, the plasma gas flow rate, nozzle diameter, and welding speed have important impacts on the welding process. Thus, an arc pressure model affected by multiple welding parameters [23] was established based on the continuous plasma impact theory. However, the additional pressure generated by arc self-magnetic compression was not considered. Guo [24] pointed out that when the additional pressure changes along the arc length direction, the medium in the arc column driven by the axial pressure gradient flowing from high to low pressure impacts the workpiece, which cannot be neglected.
In our study, an improved three-dimensional transient arc pressure model affected by multiple parameters was developed based on the principle of additional pressure produced through self-magnetic arc compression and the jet generated by plasma arc compression. In addition, the model considered the influence of arc height and keyhole secondary compression on the arc pressure, which firstly decreased linearly with the increase in the arc height and increased squarely with the increase in the keyhole depth after the keyhole formation. The influence of the perforated weld pool flow behaviors on the weld pool stability was investigated. The simulation results were verified by experiments. The weld pool stability, predicted by the material flow behavior of the perforated weld pool, is feasible and the application of the computational fluid dynamics (CFD) finite volume method (FVM) model for a material flow behavior simulation for VPPAW is practical in formulating and improving the welding technology during production.

Experimental Methods
The VPPA welding and data acquisition system was schematically shown in Figure 1, which consisted of a welding power source (VPPA-2), plasma arc welding torch (PAW-300), workpiece, high-speed camera (Red Lake Y4) and an industrial computer.
The workpiece employed in this study was an AA2219 plate with a chemical composition as shown in Table 1. The dimensions of the workpiece was 120 mm × 60 mm × 6 mm (length × width × thickness). The shielding gas used was pure argon with a volume fraction of 99.99% throughout the experiments. The high-speed camera technology was used to capture a real-time evolution image of the back keyhole. Materials 2020, 13, x doi: FOR PEER REVIEW 3 of 19 Figure 1. Schematic of variable polarity plasma arc (VPPA) welding and data acquisition system.

Model Description
The model was developed in an Ansys-Fluent environment based on FVM, which allowed the simulation of the VPPA welding process to be conducted.

Basic Assumptions
(1) Molten metal is considered as a laminar incompressible Newtonian fluid affected by the plasma arc pressure, surface tension, electromagnetic force and gravity [9].
(2) A Boussinesq approximation was used to treat the buoyancy term [27]. The welding parameters are summarized in Table 2 [6,21]. To reduce the burning loss of a tungsten electrode, a large I DCEP and a short t DCEP are suggested to satisfy high energy while reducing the tungsten electrode loss and cleaning the oxidation film [26].

Model Description
The model was developed in an Ansys-Fluent environment based on FVM, which allowed the simulation of the VPPA welding process to be conducted.

Basic Assumptions
(1) Molten metal is considered as a laminar incompressible Newtonian fluid affected by the plasma arc pressure, surface tension, electromagnetic force and gravity [9].
(2) A Boussinesq approximation was used to treat the buoyancy term [27].
(3) The thermal conductivity, specific heat, surface tension and viscosity are related to temperature, whereas the other thermo-physical parameters are constant [18].

Governing Equations
The continuous equation, momentum conservation equation, and energy conservation equation need to be solved to obtain the perforation process and the weld pool flow behavior for VPPAW.

Continuous Equation
The continuous equation can be expressed through Equation (1): where ρ is the constant fluid density, t is the time, → v is the velocity vector.

Momentum Conservation Equation
The momentum conservation equation in the x, y, z directions can be written as follows: where u, v, and w are the velocities in x, y, and z directions, respectively; ρ is the mixture density; and µ is the viscosity. The momentum source terms S x , S y , and S z in Equations (2)-(4) can be expressed as follows: where A x , A y , and A z are the momentum loss in the solid-liquid two-phase region caused by solidification in the x, y and z direction, respectively; F x , F y , and F z are the electromagnetic force in the x, y and z direction, respectively; and P a (x,y) is the arc pressure, which is described individually in Section 3.4. The momentum loss in the solid-liquid two-phase region is described as follows: In the solidification process of the weld pool, the amount of flowing liquid metal gradually decreases to zero, which leads to the momentum loss of the flowing metal. An enthalpy-porosity technique is used to deal with the momentum loss of a solid-liquid zone [28] (mushy zone). The enthalpy-porosity technique is essentially used to calculate the momentum loss in each element based on temperature. To be specific, the liquefaction volume fraction of the unit is calculated based on the temperature of the unit. The momentum loss from solidification is calculated through Equation (8): Here, A x , A y , and A z obtained through the orthogonal decomposition of Equation (8) on the x, y, and z axes are as follows: where A mush is the mushy zone constant, which is used to adjust the damping coefficient; → v is the velocity vector of the fluid, and β is the volume fraction of the liquid phase and can be expressed through the following equation.
According to Maxwell's electromagnetic theory, the electromagnetic forces [29] are expressed in a Cartesian coordinate system as follows: where r = (x − vt) 2 + y 2 and I is the welding current, which is described in Equation (17).

Energy Conservation Equation
According to the energy conservation code, the increased rate of energy in the microelement is equal to the net heat flux in the microelement, added to the mass force and the work conducted on the surface of the microelement, which can be expressed as follows [30]: where S H = h + ∆H; h = h re f + T T re f C P dT; ∆H = β L m . S H is the energy change caused by the liquid phase transition; h re f is the reference enthalpy; T re f is the reference temperature; C P is the specific heat; β is the volume fraction of the liquid phase, which can be expressed through the Equation (12), and L m is the latent heat of fusion. S v is the heat source equation and can be expressed as Equations (18) and (19).

Energy Input for VPPAW Process
The welding current in VPPAW, which changes in stages, has a significant impact on energy input into the weld, and thus the current expression should first be determined. As shown in Figure 2 below, Materials 2020, 13, 303 6 of 20 the current change stages include the preheating stage, current rising stage, variable polarity stage, and current descending stage [31]. The heat transfer model was constructed and consists of a double ellipsoid heat source on the top and a conical heat source below [32].
The heat flux distribution of the double ellipsoid heat source at the (x, y, z) points is expressed as follows: where m1 is the distribution coefficient of the double ellipsoid heat source, and are the distribution coefficients of the front and rear heat inputs of the double ellipsoid heat source, and + = 1.
The heat flux distribution of the conical heat source at the (x, y, z) points is expressed as follows [33]: where ( ) = − ( − ) ; and are the radius of the conical heat source for the top and bottom workpiece surfaces; and are the z component coordinates for the top and bottom workpiece surfaces, respectively; m2 is the distribution coefficient of the conical heat source, and m1 + m2 = 1.

Arc Pressure for VPPAW Process
The arc pressure is an important factor in the keyhole evolution. Based on the mechanism of additional pressure produced through arc self-magnetic compression and a jet generated by plasma arc mechanical compression, a three-dimensional transient arc pressure model was developed. In addition, considering that the arc pressure first decreases with the increase in the arc height, and then increases with the additional constraints of the keyhole after the keyhole formation, the arc pressure model is further optimized. The plasma arc pressure can be expressed as follows [34]: where P is the total arc pressure, is the static arc pressure, and is the dynamic arc pressure.

Static Arc Pressure
The arc column expanding from the nozzle to the workpiece is assumed to have a conical shape, and the current distribution along the solid angle of the cone is uniform. A schematic diagram of the electromagnetic static arc pressure distribution is shown in Figure 3. The current variation in Figure 2 is represented through Equation (17): where tan α = (I DCEN − I 0 )/(t S − t 0 ); tan β = I DCEN − (t N − t M ); I 0 is the preheating current, t 0 is the preheating time; t S is the current rising time from the beginning of the welding process to I DCEN ; t M is the time from the beginning of the welding process to the start of the keyhole closure phase; t N is the total welding time; and T C is the positive and negative polarity current cycle time during the variable polarity stage.
The heat transfer model S v was constructed and consists of a double ellipsoid heat source on the top and a conical heat source below [32].
The heat flux distribution of the double ellipsoid heat source at the (x, y, z) points is expressed as follows: where m 1 is the distribution coefficient of the double ellipsoid heat source, f 1 and f 2 are the distribution coefficients of the front and rear heat inputs of the double ellipsoid heat source, and The heat flux distribution of the conical heat source at the (x, y, z) points is expressed as follows [33]: where ; r e and r i are the radius of the conical heat source for the top and bottom workpiece surfaces; z e and z i are the z component coordinates for the top and bottom workpiece surfaces, respectively; m 2 is the distribution coefficient of the conical heat source, and m 1 + m 2 = 1.

Arc Pressure for VPPAW Process
The arc pressure is an important factor in the keyhole evolution. Based on the mechanism of additional pressure produced through arc self-magnetic compression and a jet generated by plasma arc mechanical compression, a three-dimensional transient arc pressure model was developed. In addition, considering that the arc pressure first decreases with the increase in the arc height, and then increases with the additional constraints of the keyhole after the keyhole formation, the arc pressure model is further optimized. The plasma arc pressure can be expressed as follows [34]: where P is the total arc pressure, P 0 is the static arc pressure, and P v is the dynamic arc pressure.

Static Arc Pressure
The arc column expanding from the nozzle to the workpiece is assumed to have a conical shape, and the current distribution along the solid angle of the cone is uniform. A schematic diagram of the electromagnetic static arc pressure distribution is shown in Figure 3.
where I is the welding current calculated through Equation (17), is half of the conical angle, is the angle between the z-axes and the axis of observation point A, and is the distance (cm) between point A and the vertex of the conical angle.
The static arc pressure at the center of the plasma arc is the largest when is a constant. Namely, the maximum static arc pressure ′′ can be obtained if = 0 and is the distance from the torch to the workpiece. K' is an adjusting coefficient which is introduced to consider the influence of the change in plasma gas flow rate on . Assuming that the workpiece under the arc is completely in a melting state, the molten pool will be concave [24] under the static arc pressure.

Dynamic Arc Pressure
The distribution of energy and force at the interface between the keyhole and the weld pool is fatal for the keyhole formation and maintenance [9,20]. During the VPPAW process, the effects of the plasma jet promote the keyhole penetration. The density of the plasma arc is approximately the same. According to Bernoulli's theorem, the dynamic arc pressure is obtained as follows: where V is the velocity of the plasma jet, is the plasma density, and is the pressure difference The static arc pressure at any position in the arc column (ι, ψ) can be expressed through the following Equation (21) [35] where I is the welding current calculated through Equation (17), θ is half of the conical angle, ψ is the angle between the z-axes and the axis of observation point A, and l is the distance (cm) between point A and the vertex of the conical angle.
The static arc pressure at the center of the plasma arc is the largest when l is a constant. Namely, the maximum static arc pressure P 0 can be obtained if ψ = 0 and l is the distance from the torch to the workpiece. K' is an adjusting coefficient which is introduced to consider the influence of the change in plasma gas flow rate on θ. Assuming that the workpiece under the arc is completely in a melting state, the molten pool will be concave [24] under the static arc pressure.

Dynamic Arc Pressure
The distribution of energy and force at the interface between the keyhole and the weld pool is fatal for the keyhole formation and maintenance [9,20]. During the VPPAW process, the effects of the plasma jet promote the keyhole penetration. The density of the plasma arc is approximately the same. According to Bernoulli's theorem, the dynamic arc pressure P v is obtained as follows: where V is the velocity of the plasma jet, ρ is the plasma density, and P v is the pressure difference between Sections A and B. Figure 4 shows the internal structure of the nozzle. We assume that the plasma in the nozzle satisfies the law of mass conservation of a fluid flow, where α is a adjusting coefficient, which was treated as a constant. Introducing into Equation (24), the dynamic arc pressure at the nozzle outlet can be obtained as follows: Because the velocities in the arc axes decrease gradually along the downstream direction, owing to the additional pressure, the gradient decreases from the nozzle to the weld pool surface. Here, can be considered as the maximum velocity in the arc column. We assume that the plasma density is a constant. The velocity of the nozzle outlet is approximately equal to that of the workpiece surface, and thus ′ can be regarded as the largest dynamic pressure on this surface.
To summarize, in VPPAW, the total arc pressure P can be obtained as follows: the arc pressure inside the workpiece conforms to a Gaussian distribution [23]: where d is the nozzle diameter; r is the radius of an arbitrary cross section; r is the radius of the plasma arc force impacting on the workpiece. Significantly, in the above arc pressure model, the current varies according to the different loading stages.

Arc Pressure Model Improvement
The plasma arc pressure was decreased at first due to the increase in the arc length, and then increased due to the additional constraint of the keyhole [20]. For simplicity, in this study, from the top to the one-sixth of the base metal, the arc pressure is assumed to be linearly decreased as Equation (27), and the arc pressure increases quadratically along the direction of the workpiece thickness, as Equation (28) from the keyhole appearance to penetration. Based on Pan's study [21], the keyhole begins to appear about 1 mm below the workpiece. The improved arc pressure model P' can expressed as following: where φ is more than zero and less than one and is introduced to consider the variation of the arc pressure with the change of weld pool depth. is the weld pool depth, ℎ is the real-time keyhole depth and L is the workpiece thickness. It was found that φ equal to 0.4 was suitable throughout the tests. After igniting the plasma arc, the velocity in Section B is calculated through the following equation: where M A and M B are the plasma molecular weight flowing through Sections A and B, respectively; B is the section area at the nozzle outlet; P A , T A , and ρ A , and P B , T B , and ρ B , are the dynamic pressure, plasma temperature, and plasma density in Sections A and B, respectively; Setting M A = M B = 39.95 g/mol and P A = P B = 1.01 × 10 5 N/m 2 . T A = 3000 K, T B = 18,000 K [36], and K is an adjusting coefficient which is introduced to consider the influence of the change in plasma gas flow rate on the ratio of T B /T A . V B can be obtained through Equation (24): where α is a adjusting coefficient, which was treated as a constant. Introducing V B into Equation (24), the dynamic arc pressure at the nozzle outlet can be obtained as follows: Because the velocities in the arc axes decrease gradually along the downstream direction, owing to the additional pressure, the gradient decreases from the nozzle to the weld pool surface. Here, V B can be considered as the maximum velocity in the arc column. We assume that the plasma density ρ Ar is a constant. The velocity of the nozzle outlet is approximately equal to that of the workpiece surface, and thus P v can be regarded as the largest dynamic pressure on this surface. To summarize, in VPPAW, the total arc pressure P can be obtained as follows: the arc pressure inside the workpiece conforms to a Gaussian distribution [23]: where d is the nozzle diameter; r is the radius of an arbitrary cross section; r p is the radius of the plasma arc force impacting on the workpiece. Significantly, in the above arc pressure model, the current varies according to the different loading stages.

Arc Pressure Model Improvement
The plasma arc pressure was decreased at first due to the increase in the arc length, and then increased due to the additional constraint of the keyhole [20]. For simplicity, in this study, from the top to the one-sixth of the base metal, the arc pressure is assumed to be linearly decreased as Equation (27), and the arc pressure increases quadratically along the direction of the workpiece thickness, as Equation (28) from the keyhole appearance to penetration. Based on Pan s study [21], the keyhole begins to appear about 1 mm below the workpiece. The improved arc pressure model P can expressed as following: where ϕ is more than zero and less than one and is introduced to consider the variation of the arc pressure with the change of weld pool depth. Z k is the weld pool depth, h kh is the real-time keyhole depth and L is the workpiece thickness. It was found that ϕ equal to 0.4 was suitable throughout the tests.

Computational Domain
The computational domain for the VPPAW model is shown in Figure 5. The domain is symmetric about the x-z plane. To improve the computational efficiency, the length (x-axes), width (y-axes) and height (z-axes) of the workpiece are 30 mm, 15 mm, and 10 mm, respectively. The domain is divided into three layers, the middle of which is a 6-mm fluid layer. The upper and lower layers are 2-mm air layers. The minimum grid size is 0.3 mm.

Computational Domain
The computational domain for the VPPAW model is shown in Figure 5. The domain is symmetric about the x-z plane. To improve the computational efficiency, the length (x-axes), width (y-axes) and height (z-axes) of the workpiece are 30 mm, 15 mm, and 10 mm, respectively. The domain is divided into three layers, the middle of which is a 6-mm fluid layer. The upper and lower layers are 2-mm air layers. The minimum grid size is 0.3 mm.

Keyhole Tracking
A two-phase (gas and aluminum alloy) flow model was used [37]. In this case, volume fraction function F = 1 corresponds to cells full of metal fluid, while F = 0 corresponds to cells empty of metal fluid. Cells with F lying between 0 and 1 locate in the free surface.

Keyhole Tracking
A two-phase (gas and aluminum alloy) flow model was used [37]. In this case, volume fraction function F = 1 corresponds to cells full of metal fluid, while F = 0 corresponds to cells empty of metal fluid. Cells with F lying between 0 and 1 locate in the free surface. Figure 6 shows a schematic diagram of the boundary conditions, which were set as indicated in Table 3 below.

Keyhole Tracking
A two-phase (gas and aluminum alloy) flow model was used [37]. In this case, volume fraction function F = 1 corresponds to cells full of metal fluid, while F = 0 corresponds to cells empty of metal fluid. Cells with F lying between 0 and 1 locate in the free surface. Figure 6 shows a schematic diagram of the boundary conditions, which were set as indicated in Table 3 below.

Boundary
Momentum Thermal Energy where ∂T/∂ → n is the temperature gradient; ∂T/∂ → y is the variance ratio of temperature along the y axes; and ∂u/∂ → y and ∂w/∂ → z are the velocity variation rates along the y and z axes, respectively.
The surface tension of the gas-liquid interface is adopted using the continuum surface force (CSF) model [38] and calculated according to the following equation [39].
where T is the real-time temperature of the weld pool.

Properties of AA2219
Some thermo-physical properties of AA2219 were set out in Table 4, as follows: Table 4. Thermo-physical properties of AA2219 [17].

Parameter Symbol Number
Solid temperature T s 816 K Liquid temperature T l 916 K Boiling temperature T b 2740 K Latent heat of fusion L m 391,000 J/(kg·K) For AA2219, the temperature dependence of density is represented by the following linear equation [39]: The specific heat capacity and viscosity of AA2219 are set through the following equations [17]: The thermal conductivity of AA2219 is calculated by the following Equation [38]: In Equations (31)-(33), T is the temperature of the weld pool. In Equation (34), R 0 is the ideal gas constant, N 0 is the Avogadro constant, and M is the molar mass of aluminum.

Results and Discussion
Firstly, the numerical results of case A were analyzed, including the dynamic change of the keyhole size, the arc pressure distribution on the keyhole surface, the temperature field and the flow behavior of the molten metal in the weld pool. Figure 7 shows the arc pressure distribution during the keyhole evolution for case a. In Figure 7a-h, W is half of the keyhole width on the upper workpiece surface, and D is the keyhole depth. The maximum arc pressure appears at the bottom of the keyhole, reaching 3000 Pa. Jiang et al. [40] found that the arc pressure is between 2500 and 4500 Pa, and thus the arc pressure calculated in the present study is reasonable. At 1.72 s, a small pool is formed with a half width of 1.2 mm and a half depth of 0.7 mm. Based on the thermo-mechanical coupling effect, at 2.43 s, the arc pressure penetrates into the weld pool to form a keyhole with a half width of 2.0 mm on the upper surface and a depth of 1.8 mm. Subsequently, the keyhole size increases continuously, and the increasing rate of the keyhole depth is greater than that of the keyhole width. This is because of the secondary compression of the keyhole to the arc pressure. The workpiece is penetrated at 4.55 s. Figure 8 shows the material flow and temperature distribution during the keyhole evolution for case a. The arrows represent the directions of the material flow. At first, the workpiece surface is heated by the plasma and depressed by the arc pressure. As the heat accumulates, at 1.72 s, the molten metal forms eddies, as shown in Figure 8a    The simulation results show that the maximum arc pressure during the DCEN and DCEP phase are 2600 Pa and 3000 Pa as shown in Figure 9a,b respectively.
In order to verify the accuracy of the arc pressure model, an arc pressure measure and data acquisition system were set up schematically, as shown in Figure 10.
The system consisted of a welding power source (VPPA-2), a plasma arc welding torch (PAW-300), a pressure transmitter (YB005-01) developed by General Engineering Research Institute of China Academy of Sciences, an Agilent memory oscilloscope (6422D) and an industrial computer data acquisition card (PCL1800) produced by Advantech. The welding power source was the self-developed VPPA-2 variable polarity plasma inverter with an 80C196KC single chip microcomputer as the control core and the main circuit as the double inverse circuit topological structure. The plasma arc pressure was transmitted to the pressure transmitter through the copper plate with the keyhole and the copper tube. The output voltage signal of the pressure transmitter was collected synchronously by the industrial computer and the memory oscilloscope. The simulation results show that the maximum arc pressure during the DCEN and DCEP phase are 2600 Pa and 3000 Pa as shown in Figure 9a,b respectively. In order to verify the accuracy of the arc pressure model, an arc pressure measure and data acquisition system were set up schematically, as shown in Figure 10.  The simulation results show that the maximum arc pressure during the DCEN and DCEP phase are 2600 Pa and 3000 Pa as shown in Figure 9a,b respectively. In order to verify the accuracy of the arc pressure model, an arc pressure measure and data acquisition system were set up schematically, as shown in Figure 10.   The simulation results show that the maximum arc pressure during the DCEN and DCEP phase are 2600 Pa and 3000 Pa as shown in Figure 9a,b respectively. In order to verify the accuracy of the arc pressure model, an arc pressure measure and data acquisition system were set up schematically, as shown in Figure 10.  The signal converter converts the voltage signals into digital signals. The waveform measured by the oscilloscope corresponds to the arc pressure. As Figures 11 and 12 show, the vertical voltage axis is 2000 Pa/grid, and the horizontal time axis is 30 ms/grid. The waveform shows that the average values of the arc pressure during the DCEN phase and DCEP phase are about 2400 Pa and 2700 Pa, respectively. The measured results are consistent with the simulation results, and the accuracy of the arc pressure model is verified.
The signal converter converts the voltage signals into digital signals. The waveform measured by the oscilloscope corresponds to the arc pressure. As Figures 11 and 12 show, the vertical voltage axis is 2000 Pa/grid, and the horizontal time axis is 30 ms/grid. The waveform shows that the average values of the arc pressure during the DCEN phase and DCEP phase are about 2400 Pa and 2700 Pa, respectively. The measured results are consistent with the simulation results, and the accuracy of the arc pressure model is verified.  The calculated keyhole diameter is 1.86 mm, which is consistent with the measured value of 1.82 mm.
The calculated fusion line (FL) is consistent with the measured fusion line obtained from the experiment, as shown in Figure 13. The signal converter converts the voltage signals into digital signals. The waveform measured by the oscilloscope corresponds to the arc pressure. As Figures 11 and 12 show, the vertical voltage axis is 2000 Pa/grid, and the horizontal time axis is 30 ms/grid. The waveform shows that the average values of the arc pressure during the DCEN phase and DCEP phase are about 2400 Pa and 2700 Pa, respectively. The measured results are consistent with the simulation results, and the accuracy of the arc pressure model is verified.  The calculated keyhole diameter is 1.86 mm, which is consistent with the measured value of 1.82 mm.
The calculated fusion line (FL) is consistent with the measured fusion line obtained from the experiment, as shown in Figure 13. The calculated keyhole diameter is 1.86 mm, which is consistent with the measured value of 1.82 mm.
The calculated fusion line (FL) is consistent with the measured fusion line obtained from the experiment, as shown in Figure 13. For case b-f, simulation tests and welding experiments were successively carried out with the welding parameters in Table 2. The comparison results in Table 5 show that the calculated results were in good agreement with the experimental results and the accurate thermo-mechanical coupling model lays the foundation for the subsequent results analysis.  For case b-f, simulation tests and welding experiments were successively carried out with the welding parameters in Table 2. The comparison results in Table 5 show that the calculated results were in good agreement with the experimental results and the accurate thermo-mechanical coupling model lays the foundation for the subsequent results analysis. For case a-f, the simulation results of the material flow of the perforated weld pool and phase distribution are shown in Figure 14, and there are three kinds of flow behavior. Since the temperature gradients inside the weld pools, as shown in Figure 15, have significant impacts on the material flow behaviors of the perforated weld pools, the two figures will be discussed relevantly below.   For cases a-f, the temperature gradient inside the weld pool was studied to further explore the mechanism of the weld pool material flow behavior to weld pool stability. The temperature value of the z = −5 mm straight line on the xoz plane was extracted. Since the weld pool is symmetrical, the temperature value within 3.5 mm (from the keyhole center to its right side) is shown in Figure 15.   Due to the different horizontal distances from the weld pool boundary to the x-axis for cases a-f, the boundaries of the perforated weld pool were shifted to x = 6 mm in order to compare the slope of the temperature gradients. As shown in Figure 15b, the slope of case d was the largest, the angle between the slope line of cases a-b and the x negative direction horizontal line are 35 • (Tan (35 • ) = 0.70) and 41 • (Tan (41 • ) = 0.87), respectively, and the slope of case f was the smallest.
In Figure 14a-f, under the action of the surface tension and the digging effect of the arc pressure, the upper molten metal of the weld pool flows outwards and then forms clockwise, as shown in Figure 14a 1 ; The molten metal in the middle of the weld pool flows upwards and inwards, as shown in Figure 14a 2 , under the digging action of the arc pressure, the electromagnetic force, and the effect of clockwise eddy flow in Figure 14a 1 ; in particular, the material flow behavior in the lower part of the weld pool is quite different. In case a, case c and case e, the molten metal forms confluences towards the keyhole center, as shown in Figure 14a 3 , at the maximum velocities of 0.45 m/s, 0.55 m/s and 0.6 m/s, respectively, with moderate temperature gradients (the slopes values being greater than 0.7 and less than 0.87), as shown in Figure 15b. In case b and case d, the molten metal of the weld pool bottom flows downwards and outwards, as shown in Figure 14b 1 , at the maximum speed of 0.7 m/s and 0.8 m/s, respectively, under the action of the large surface tension caused by the relatively large temperature gradients (slope value ≥ 0.87) in the weld pool, as shown Figure 15b. In case f, the molten metal flows as shown in Figure 14f 1 2 , which makes the molten metal of the weld pool increase under the joint action of small arc pressure and surface tension. Although the weld pool metal can still form the confluence, as shown in Figure 14 f 3 , the weld pool with a large amount of molten metal collapsed under the action of gravity.
In order to summarize the simulation results, according to the material flow behavior in the lower part of the weld pool, the equilibrium coefficient η is introduced.
where Q is the plasma gas flow rate, its unit being mL/s. When 1/4 ≤ η < 1/3, such as in case a, case c and case e, the material flows to form confluences towards the keyhole center with the moderate velocities (0.45 m/s, 0.55 m/s and 0.6 m/s), and the weld pool can exist stably. When η ≥ 1/3, such as in case b and case d, the material flows downwards and outwards with the maximum velocities of about 0.7 m/s or 0.8 m/s, and the weld pools are cut. When η < 1/4, the maximum downward flow velocity of the perforated weld pool was 0.4 m/s and the weld pool with too large an amount of molten metal, caused by the long-time metal refluxing phenomenon, collapsed downward under the action of a larger gravity. Figure 16a-c are the weld formations obtained from the experiments classified by η. When η ≥ 1/3 or η < 1/4, the welds were cut or collapsed, as shown in Figure 16a,c, respectively, caused by the weld pool instability. When 1/4 ≤ η < 1/3, the weld pools can stably exist and the welds can form well, as shown in Figure 16b.
When 1/4 ≤ < 1/3, such as in case a, case c and case e, the material flows to form confluences towards the keyhole center with the moderate velocities (0.45 m/s, 0.55 m/s and 0.6 m/s), and the weld pool can exist stably. When ≥ 1/3, such as in case b and case d, the material flows downwards and outwards with the maximum velocities of about 0.7 m/s or 0.8 m/s, and the weld pools are cut. When < 1/4, the maximum downward flow velocity of the perforated weld pool was 0.4 m/s and the weld pool with too large an amount of molten metal, caused by the long-time metal refluxing phenomenon, collapsed downward under the action of a larger gravity. Figure 16a-c are the weld formations obtained from the experiments classified by ƞ. When ≥ 1/3 or < 1/4, the welds were cut or collapsed, as shown in Figure 16a,c, respectively, caused by the weld pool instability. When 1/4 ≤ < 1/3, the weld pools can stably exist and the welds can form well, as shown in Figure 16b. The different weld formations further verify that it is feasible to predict the weld pool stability by the material flow behavior of the weld pool. (2) Based on the mechanism of the additional pressure produced through self-magnetic arc compression and the jet generated from mechanical plasma arc compression, and meanwhile considering that the arc pressure first decreases with the increase in arc height, and then increases with the additional constraints of the keyhole after the keyhole formation, a three-dimensional transient model of the VPPA arc pressure which first decreases linearly with the increase in arc height and then increases squarely with the increase in keyhole depth was established. The different weld formations further verify that it is feasible to predict the weld pool stability by the material flow behavior of the weld pool. (2) Based on the mechanism of the additional pressure produced through self-magnetic arc compression and the jet generated from mechanical plasma arc compression, and meanwhile considering that the arc pressure first decreases with the increase in arc height, and then increases with the additional constraints of the keyhole after the keyhole formation, a three-dimensional transient model of the VPPA arc pressure which first decreases linearly with the increase in arc height and then increases squarely with the increase in keyhole depth was established.

Conclusions
(3) The measured arc pressure value, the shape of the fusion line, and the keyhole size on the back of the weld were in good agreement with the simulation results.

Conflicts of Interest:
The authors declare no conflict of interest. Radiation emissivity -L m Latent heat of fusion J kg −1 C p Specific heat J kg −1 K −1 x, y, z

Nomenclature
Coordinate of x, y and z axes m L Workpiece thickness m k Thermal conductivity W m −1 k −1 m 1 , m 2 , f 1 , f 2 Energy distribution coefficient a f , a r b, c, r i , r e Distribution radius of heat source m