A numerical simulation method of laser-induced phase explosion based on finite element method

Abstract

The present invention discloses a numerical simulation method for laser-induced phase explosion based on the finite element method, which relates to the technical field of high-power laser and material interaction. Laser parameters are preset to establish a numerical simulation model for nanosecond laser ablation of metals. Laser energy is used as a heat source term to generate a two-dimensional numerical grid for the calculation domain of metal ablation. Phase change heat transfer theory is used to consider the metal's phase change, phase interface movement, and particle evaporation backscattering effects to simulate the physical process and obtain the saturated vapor pressure, hot liquid pressure, and material surface temperature. Based on these saturated vapor pressure, hot liquid pressure, and material surface temperature, the numerical simulation model is discretized and solved using the finite element method. Particle tracking is used to simulate the gas-liquid phase splashing process at the moment of explosion. The present invention can accurately predict the dynamic evolution characteristics of phase explosions during laser ablation, providing important theoretical support and technical assurance for fields such as laser processing, laser propulsion, and laser micro-nano processing.

Description

Technical Field

[0001] The present invention relates to the technical field of interaction between high-power laser and materials, and more particularly to a numerical simulation method of laser-induced phase explosion based on finite element method. Background Art

[0002] With the rapid development of laser technology, laser-induced phase explosions (LIPDs) have found widespread application in materials processing, laser propulsion, pulsed laser deposition, and other fields. However, because LIDs involve the coupling of multiple physical fields (such as heat conduction, phase change, vaporization, and particle injection), existing research methods still have limitations in accurately predicting the evolution of the explosion. Traditional experimental techniques are limited by their temporal resolution and spatial capture capabilities, while single theoretical models struggle to fully describe the complex behavior of LIDs.

[0003] Therefore, developing a high-precision numerical simulation method to study the dynamic evolution characteristics of laser-induced phase explosion is an urgent problem that needs to be solved by technicians in this field. Summary of the Invention

[0004] In view of this, the present invention provides a numerical simulation method for laser-induced phase explosion based on the finite element method. By constructing a numerical simulation framework for multi-physical field coupling, the present invention can accurately predict the dynamic evolution characteristics of phase explosion during laser ablation, especially considering the interaction between laser-induced phase explosion and metal ablation and gasification effects, providing important theoretical support and technical guarantees for laser processing, laser propulsion, laser micro-nano processing and other fields.

[0005] In order to achieve the above object, the present invention adopts the following technical solutions:

[0006] A numerical simulation method of laser-induced phase explosion based on finite element method, comprising:

[0007] Preset laser parameters and establish a numerical simulation model for nanosecond laser ablation of metals;

[0008] The laser energy obtained by the numerical simulation model is used as a heat source term to mesh the calculation domain of metal ablation to generate a two-dimensional numerical grid;

[0009] In the two-dimensional numerical grid, the phase change heat transfer theory is used to consider the phase change of the metal, the phase interface movement and the particle evaporation backscattering effect to simulate the physical process and obtain the saturated vapor pressure, hot liquid pressure and material surface temperature;

[0010] According to the saturated vapor pressure, hot liquid pressure and material surface temperature, the numerical simulation model is discretized and solved in combination with the finite element method to obtain the spatial and temporal evolution characteristics of the phase explosion particles, the phase change depth and the particle injection characteristics; finally, the spatial distribution characteristics and temporal evolution laws of the temperature and pressure in the metal liquid phase region under laser induction are obtained, and the gas-liquid phase splashing process at the moment of the explosion is simulated through particle tracking.

[0011] Preferably, the laser parameters include: laser peak power, laser pulse width, laser beam radius, and target material state parameters including metal density, specific heat capacity, and thermal conductivity.

[0012] Preferably, the physical process in the physical process simulation is divided into: heat conduction temperature rise stage, solid-liquid phase change stage, gas-liquid surface evaporation stage and liquid phase superheated boiling nucleation stage.

[0013] Preferably, the heat flow and heat flux density are calculated during the heat conduction temperature rise stage, including:

[0014] Heat flow:

[0015] Q1=AhΔT;

[0016] Where Q1 represents the heat flow, h represents the surface heat transfer coefficient, A represents the heat transfer area, and ΔT represents the surface temperature difference;

[0017] Thermal radiation:

[0018]

[0019] Where T is the surface temperature of the material, T0 is the ambient temperature, q r is the heat flux density, χ≈5.67×10 -8 W / (m 2 ·K 4 ) is a constant and ε is the surface emissivity.

[0020] Preferably, the heat capacity C of the gas-liquid phase change stage is calculated during the solid-liquid phase change stage. p ,include:

[0021] C p =C p0 +L m D m +L v D v ;

[0022] Among them, C p0 is a constant heat capacity, L m represents the latent heat of fusion, The overall latent heat of fusion, T m represents the melting point of the target material, T is the surface temperature of the material, ΔT mIndicates the melting transition zone, L v represents the latent heat of vaporization, represents the latent heat of evaporation coefficient, T b Indicates the melting point of the target, ΔT b Indicates the evaporation transition zone.

[0023] Preferably, the enthalpy variable ΔH is calculated during the evaporation phase of the gas-liquid surface v , saturated vapor pressure P sat and steam recoil pressure P r , specifically including:

[0024] Enthalpy change ΔH v , calculated using the Clausius-Clapeyron equation:

[0025]

[0026] Where, T b Indicates the boiling temperature at normal pressure, T C is the critical temperature of the target, ΔH v0 Typical temperature T b The phase change enthalpy under , T is the surface temperature of the material;

[0027] Saturated vapor pressure P sat :

[0028]

[0029] Where ΔV is the specific volume, ΔH v Substitute the Clausius-Clapeyron equation and obtain the saturated vapor pressure by iterative integration;

[0030] Steam recoil pressure P r :

[0031]

[0032] Where Ts is the interface temperature, β R is the backscatter coefficient.

[0033] Preferably, the boiling nucleation stage calculates the pressure difference ΔP inside and outside the phase interface and the critical nucleation free energy W c and the nucleation rate per unit volume J, including:

[0034]

[0035] Among them, r is the radius of the gas core, σ is the surface tension coefficient, P v represents the internal critical nucleation pressure, P l represents the external pressure, r crepresents the critical nucleation radius, m is the mass of a single molecule, H v is the molar enthalpy, T is the surface temperature of the material, k B is the Boltzmann constant, and R is the universal gas constant.

[0036] Preferably, the discretization and solution of the numerical simulation model based on the saturated vapor pressure, hot liquid pressure and material surface temperature in combination with the finite element method specifically includes:

[0037]

[0038] in, is the gradient operator, u is the velocity vector, temperature gradient, K is the Darcy drag coefficient, I is the unit matrix, p is the pressure, is the divergence of velocity, is the velocity gradient, u x 、u y with u z are the components of the velocity vector u in the x, y, and z directions, respectively, and C p is the heat capacity, ρ is the target density, k is the target thermal conductivity; Q is the sum of other heat sources, which is composed of the heat taken away by the evaporation of gas-liquid phase change, the energy loss caused by convection heat transfer and thermal radiation; α is the dynamic absorption rate of the target, t p is the laser pulse width, t0 is the initial injection delay, E is the laser energy, r0 is the laser irradiation radius, x is the vertical direction of the laser beam, t is the time, g is the acceleration of gravity, and μ is the dynamic viscosity.

[0039] It can be seen from the above technical solution that, compared with the prior art, the present invention discloses a numerical simulation method of laser-induced phase explosion based on the finite element method, which can provide high-precision numerical simulation of laser-induced phase explosion, significantly improving the accuracy of phase explosion characteristic prediction, especially in the precise description of complex physical phenomena under the action of high-power lasers. BRIEF DESCRIPTION OF THE DRAWINGS

[0040] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the following briefly introduces the drawings required for use in the embodiments or the description of the prior art. Obviously, the drawings described below are merely embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on the provided drawings without paying any creative work.

[0041] Figure 1 A diagram of the method steps provided by the present invention;

[0042] Figure 2 A graph showing the relationship between the homogeneous nucleation rate and temperature provided by the present invention;

[0043] Figure 3 A grid division diagram provided by the present invention;

[0044] Figure 4 The evolution result diagram provided by the present invention;

[0045] Figure 5 The temperature rise curve provided by the present invention;

[0046] Figure 6 The temperature rise rate diagram provided by the present invention;

[0047] Figure 7 A spatiotemporal distribution diagram of particle explosion products provided by the present invention;

[0048] Figure 8 The particle distribution result diagram of the phase explosion space provided by the present invention;

[0049] Figure 9 A physical process simulation diagram provided by the present invention;

[0050] Figure 10 A simulation diagram of the gas-liquid phase splashing process provided by the present invention;

[0051] Figure 11 This is a schematic diagram of the PVT relationship of the van der Waals equation of state provided by the present invention. DETAILED DESCRIPTION

[0052] The following will clearly and completely describe the technical solutions in the embodiments of the present invention in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts are within the scope of protection of the present invention.

[0053] The embodiment of the present invention discloses a numerical simulation method of laser-induced phase explosion based on finite element method, such as Figure 1 Shown, including:

[0054] S1: Preset laser parameters and establish a numerical simulation model for nanosecond laser ablation of metals;

[0055] S2: Using the laser energy obtained from the numerical simulation model as a heat source term, the calculation domain of metal ablation is meshed to generate a two-dimensional numerical grid;

[0056] S3: In a two-dimensional numerical grid, phase change heat transfer theory is used to simulate the physical process, considering the phase change of the metal, phase interface movement, and particle evaporation backscattering effect, to obtain the saturated vapor pressure, hot liquid pressure, and material surface temperature;

[0057] S4: Based on the saturated vapor pressure, hot liquid pressure and material surface temperature, the finite element method is used to discretize and solve the numerical simulation model and the phase transition depth, and finally the spatial distribution characteristics and time evolution law of the temperature and pressure in the metal liquid phase region under laser induction are obtained.

[0058] In a specific embodiment, the numerical simulation model of nanosecond laser ablation of metals takes into account the spatial distribution of the molten pool temperature and pressure in the liquid phase overheating area.

[0059] In a specific embodiment, S4 specifically includes: obtaining the distribution of superheated liquid area, temperature, and pressure in the numerical simulation model of nanosecond laser ablation of metal, and using the inherited decoupling method to use the distribution of superheated liquid area, temperature, and pressure as the initial conditions of this stage. The discrete phase model is used to calculate its nucleation and growth process to obtain the distribution of the discrete phase mass fraction with temperature, pressure and material physical parameters in the current system, and the interface of its transmission effect is given. Further, the distribution of the discrete phase mass fraction with temperature, pressure and material physical parameters is brought into the fluid transport model, and the discrete phase interface is introduced to describe the interaction process between the continuous phase and the discrete phase. Finally, the effect result is used by particle tracking to simulate the gas-liquid phase splashing process at the moment of the explosion, such as Figure 10 shown.

[0060] In a specific embodiment, the laser parameters include: laser peak power, laser pulse width, laser beam radius, and target material state parameters include metal density, specific heat capacity, and thermal conductivity.

[0061] In a specific embodiment, Figure 9 As shown in FIG, the physical process in the physical process simulation is divided into: heat conduction temperature rise stage, solid-liquid phase change stage, gas-liquid surface evaporation stage and liquid phase superheated boiling nucleation stage.

[0062] In a specific embodiment, calculating the heat flow and heat flux density during the heat conduction temperature rise phase includes:

[0063] Heat flow:

[0064] Q1=AhΔT;

[0065] Where Q1 represents the heat flow, h represents the surface heat transfer coefficient, A represents the heat transfer area, and ΔT represents the surface temperature difference;

[0066] Thermal radiation:

[0067]

[0068] Where T is the surface temperature of the material, T0 is the ambient temperature, q r is the heat flux density, χ≈5.67×10 -8 W / (m 2 ·K 4 ) is a constant and ε is the surface emissivity.

[0069] In a specific embodiment, the heat capacity C of the gas-liquid phase transition stage is calculated from the solid-liquid phase transition stage. p ,include:

[0070] C p =C p0 +L m D m +L v D v ;

[0071] Among them, C p0 is a constant heat capacity, L m represents the latent heat of fusion, The overall latent heat of fusion, T m represents the melting point of the target material, T is the surface temperature of the material, ΔT m Indicates the melting transition zone, L v represents the latent heat of vaporization, represents the latent heat of evaporation coefficient, T b Indicates the melting point of the target, ΔT b Indicates the evaporation transition zone.

[0072] In a specific embodiment, the enthalpy change ΔH is calculated during the evaporation phase of the gas-liquid surface. v , saturated vapor pressure P sat and steam recoil pressure P r , specifically including:

[0073] Enthalpy change ΔH v , calculated using the Clausius-Clapeyron equation:

[0074]

[0075] Where, T b Indicates the boiling temperature at normal pressure, T C is the critical temperature of the target, ΔH v0 Typical temperature T b The phase change enthalpy under , T is the surface temperature of the material;

[0076] Saturated vapor pressure P sat :

[0077]

[0078] Where ΔH v Substitute the Clausius-Clapeyron equation and obtain the saturated vapor pressure by iterative integration, where ΔV is the specific volume.

[0079] Steam recoil pressure P r :

[0080]

[0081] Where Ts is the interface temperature, β R is the backscatter coefficient.

[0082] In a specific embodiment, the pressure difference ΔP inside and outside the phase interface and the critical nucleation free energy W are calculated in the boiling nucleation stage. c and the nucleation rate per unit volume J, including:

[0083]

[0084] Among them, r is the radius of the gas core, σ is the surface tension coefficient, P v represents the internal critical nucleation pressure, P l represents the external pressure, r c represents the critical nucleation radius, m is the mass of a single molecule, H v is the molar enthalpy, T is the surface temperature of the material, k B is the Boltzmann constant, and R is the universal gas constant.

[0085] In a specific embodiment, the discretization and solution of the numerical simulation model based on the saturated vapor pressure, the hot liquid pressure and the material surface temperature in combination with the finite element method specifically includes:

[0086]

[0087] in, is the gradient operator, u is the velocity vector, temperature gradient, K is the Darcy drag coefficient, I is the unit matrix, p is the pressure, is the divergence of velocity, is the velocity gradient, u x 、u y with u z are the components of the velocity vector u in the x, y, and z directions, respectively, and C pis the heat capacity, ρ is the target density, k is the target thermal conductivity; Q is the sum of other heat sources, which is composed of the heat taken away by the evaporation of gas-liquid phase change, the energy loss caused by convection heat transfer and thermal radiation; α is the dynamic absorption rate of the target, t p is the laser pulse width, t0 is the initial injection delay, E is the laser energy, r0 is the laser irradiation radius, x is the vertical direction of the laser beam, t is the time, g is the acceleration of gravity, and μ is the dynamic viscosity.

[0088] In a specific embodiment, during each stage of the physical process simulation:

[0089] Heat conduction temperature rise stage

[0090] When the laser is incident on the target surface, reflection, absorption and transmission will occur. According to the law of energy conservation and Fresnel's formula, the absorption rate of the target to the laser can be expressed as:

[0091]

[0092] Where ε0 represents the dielectric constant of a vacuum, ω represents the angular frequency of the incident light, and c represents the speed of light in a vacuum. ρ represents the material resistivity, and λ represents the wavelength of the incident laser. This is then used as the target's absorptivity for the laser beam and then used in the Gaussian heat source equation.

[0093] After the target material absorbs the laser, processes such as temperature rise, convection and radiation will occur.

[0094] Heat conduction:

[0095] According to Fourier's law:

[0096]

[0097] Where q represents energy density and k is the thermal conductivity of the material.

[0098] convection:

[0099] According to Newton's cooling formula:

[0100] Q1=AhΔT;

[0101] Where Q1 represents the heat flow, h represents the surface heat transfer coefficient, A represents the heat exchange area, and ΔT represents the surface temperature difference.

[0102] Thermal radiation:

[0103] According to the Stefan-Boltzmann law:

[0104]

[0105] Where T is the surface temperature of the material, T0 is the ambient temperature, q r is the heat flux density, σ≈5.67×10-8 W / (m 2 ·K 4 ) is a constant and ε is the surface emissivity.

[0106] Solid-liquid phase transition stage

[0107] After the target absorbs the laser energy, the temperature rises to the solid-liquid phase transition point. According to the phase change heat transfer theory, the heat transfer equation will be modified as follows:

[0108]

[0109] In the formula, according to the laser ablation area, the reverse direction of the laser beam is determined as the z-axis, and the horizontal direction is determined as the r-axis.

[0110] f s is the solid phase ratio, L is the latent heat of phase change. s =1, the target is a solid pure phase. Then:

[0111]

[0112] Substituting the phase change heat transfer equation into the equation, we can get:

[0113]

[0114] Among them, due to the existence of the phase change transition zone, the latent heat of phase change needs to be dealt with, so the specific heat capacity needs to be described using the dynamic specific heat capacity method:

[0115] C p =C p0 +L m D m ;

[0116] Where C p0 is a constant heat capacity, L m represents the latent heat of fusion, The overall latent heat of melting coefficient T m Indicates the melting point of the target, ΔT m Indicates the melting transition zone.

[0117] The heat capacity for the gas-liquid phase change stage will be modified to:

[0118] C p =C p0 +L m D m +L v D v ;

[0119] Where, L v represents the latent heat of vaporization, represents the latent heat of evaporation coefficient, T b Indicates the melting point of the target, ΔTb Indicates the evaporation transition zone.

[0120] The above relationship combined with the formula can be used to obtain the governing equations of the solid-liquid interface, where the mass conservation and momentum conservation formulas are the fluid transport model:

[0121] Conservation of mass:

[0122]

[0123] Conservation of Momentum:

[0124]

[0125] Conservation of Energy:

[0126]

[0127] Where T is temperature, is the velocity vector, ρ is the density, C p is the specific heat capacity, P is the pressure, is the unit matrix, k is the thermal conductivity, and μ is the dynamic viscosity. is the Darcy friction force:

[0128]

[0129] Where K is the Carman-Kozeny permeability coefficient, C is the Carman-Kozeny constant, and b is a minimum constant to prevent the denominator from being equal to 0. ε is the porosity, i.e., the liquid volume fraction during the solid-liquid phase transition:

[0130]

[0131] For the gas-liquid interface, consider the Marangoni effect:

[0132]

[0133] γ is the surface tension coefficient, is the normal vector to the gas-liquid interface.

[0134] Gas-liquid surface evaporation stage

[0135] 1) Saturated vapor pressure prediction

[0136] The more commonly used form in laser ablation is the Clausius-Clapeyron equation:

[0137]

[0138] The present invention focuses on simulating the laser-induced phase explosion process. The conventional simplified saturated vapor pressure formula will deviate significantly under high temperature conditions. Therefore, the Watson equation is used to describe the enthalpy variable ΔH in the phase change process. v

[0139]

[0140] Where, T b Indicates the boiling temperature at normal pressure, T C is the critical temperature of the target, ΔH v0 Typical temperature T b The phase change enthalpy under standard atmospheric pressure is usually taken as the value corresponding to the boiling point. v Substitute the Clausius-Clapeyron equation and obtain the saturated vapor pressure through iterative integration.

[0141] 2) Study on the dynamic process of Knudsen layer

[0142] Since the physical quantities in the Knudsen layer have large jump changes, it is difficult to describe their spatial continuity. Therefore, the distribution function is used to represent the process:

[0143]

[0144] Where f( + ),f (-) Represent the forward scattering and backscattering distribution functions respectively, Vx, Vy, Vz are the components of the particle velocity in each direction, n s is the saturated vapor number density at the interface, Ts is the interface temperature, P sat (T s ) is the interface saturated vapor pressure, u is the backscattered particle velocity, k B is the Boltzmann constant, m is the mass of the metal atom, ρ v , Tv are steam pressure and temperature respectively, β R is the backscattering coefficient. Under the condition of phase equilibrium, that is, when the saturated vapor pressure is equal to the external pressure, the number of particles emitted from the gas-liquid interface is the same as the number of particles returned, and the net flux is zero. As the evaporation intensity increases, the difference in interface imbalance distribution gradually increases, so β ​​is introduced. R is the backscattering coefficient, considering the backscattering in Knudsen during strong evaporation.

[0145] According to Samokhin's research:

[0146]

[0147] Where γ is the gas specific heat ratio, U is the steam flow velocity, Mach is the Mach number.

[0148] The pressure distribution of metal vapor on the gas-liquid interface. This pressure is mainly generated by two parts: 1) the reaction of particle evaporation movement; 2) the effect of backscattered particles on the interface. Introducing the backscattering coefficient can be obtained:

[0149]

[0150] Where P r is the steam recoil pressure.

[0151] Obtain the enthalpy change ΔH during the phase change of the relevant parameters v ; Saturated vapor pressure Psat; Steam recoil pressure P r .

[0152] Liquid superheating stage

[0153] The superheated state of a liquid refers to the phenomenon that the liquid is higher than the boiling point corresponding to the current ambient pressure but does not boil. Based on the ideal gas state equation:

[0154] PV = nRT;

[0155] When considering the actual volume of gas molecules, the gas molar volume V is introduced m and the volume b of a single molar molecule, we get:

[0156] P(V m -b)=RT;

[0157] The ideal gas state equation is further modified by introducing a constant a. represents the interaction between molecules, then:

[0158]

[0159] Since the van der Waals equation has obvious misestimation near the critical temperature, Maxwell's equal area method is used to correct it, and finally its simplified form can be obtained:

[0160]

[0161] Where, P V Indicates the saturated vapor pressure corresponding to the equilibrium position, V g is the volume of pure gas phase on the equilibrium line, V lis the volume of pure liquid phase at the corresponding equilibrium point. Maxwell's equal area method can effectively correct the error near the critical point, but it is difficult to solve its cubic equation. Therefore, the root-free method obtained by Carl W. David using intermediate variable substitution is adopted here. By replacing the end pressure of the isotherm liquid with the liquid-gas conversion pressure point, T and V can be obtained. g 、V l The relationship is:

[0162]

[0163] Further introduce the intermediate variable e d , we can finally get:

[0164]

[0165] in:

[0166]

[0167] Finally, we can get the relationship between T / Tc, V / Vc, P / Pc, which is the simplified form of the van der Waals equation of state Tr, Vr, and Pr, as shown in the following example: Figure 11 As shown, it can be seen that the left end of the critical point is the saturated liquid equilibrium curve, the right end of the critical point is the saturated vapor equilibrium curve, and ABDC is 0.9T c Isotherms. The liquid and gas saturation lines converge at the critical point, resulting in a "smooth" transition. This phenomenon is consistent with the fact that the phase change enthalpy becomes zero at the critical point as described above. Therefore, the pure gas phase volume on the equilibrium line during the liquid superheating stage and the pure liquid phase volume at the corresponding equilibrium point are not used in the simulation of the present invention.

[0168] To further illustrate the overheating issue, extracting the coordinates of points D' and D from the curve in the figure reveals that while the pressure at these two points is equal, there are differences in temperature and gas volume. Furthermore, since D' is a point on the saturated liquid curve, point D can be considered the state reached by a saturated liquid system through an isobaric temperature increase. Corresponding to the aforementioned issue, it can be seen that point D has not yet transitioned to the saturated vapor curve, indicating that the liquid is now superheated and in a "metastable state." After this point, point D may spontaneously transition along the isotherm to point C, returning to the gas phase steady state. Alternatively, external pressure may be applied to reach a new liquid phase steady state, point B.

[0169] The computational study of the overheating problem shows that the overheating condition can exist at any position before reaching the critical point. The reason is the lack of nucleation sites. Usually when the system reaches the critical point, the gas-liquid phase will undergo a spontaneous transition. This conclusion has a great impact on the subsequent thermodynamic parameters in T cOn the other hand, due to the existence of metastable states, overheating becomes more important as laser power increases during laser ablation. This has a crucial impact on determining the laser ablation threshold and studying ablation efficiency. Finally, the instability of liquids in metastable states will provide leading calculation results for studies such as homogeneous nucleation and explosive boiling.

[0170] The hot liquid area, temperature, and pressure distribution are obtained, which are used as the initial conditions for the boiling nucleation stage. The discrete phase model is used to calculate the superheated phase nucleation and growth process.

[0171] Boiling nucleation stage

[0172] For the boiling nucleation stage, the discrete phase model is used for calculations. The formulas included in this stage are the discrete phase model, and the hot liquid area, temperature, and pressure distribution are used as the initial values ​​of the discrete phase model. According to equilibrium theory, the equilibrium temperature of the system depends on the pressure and volume changes. It is described by the Young-Laplace equation:

[0173]

[0174] Where ΔP is the pressure difference between the inside and outside of the phase interface (external pressure minus internal pressure), r is the radius of the gas core, and σ is the surface tension coefficient. This equation relates surface tension to the pressure difference to describe the external work required during the formation and growth of the gas core. To further study the nucleation problem during the phase explosion process, according to classical nucleation theory, we know that:

[0175]

[0176] Where r c represents the critical nucleation radius, W c represents the critical nucleation free energy. Considering the nucleation process of superheated liquid under non-equilibrium state, the Young-Laplace relation is introduced to obtain:

[0177]

[0178] Where, P v represents the internal critical nucleation pressure, P l represents the external pressure. For the internal critical nucleation pressure, under non-equilibrium conditions, considering that the saturated vapor pressure is equal to the liquid saturation pressure, we can obtain:

[0179]

[0180] Among them, P s is the saturated vapor pressure in the current equilibrium state, P lis the superheated liquid pressure, combining the above formulas we can get:

[0181]

[0182] From this, we can obtain the relationship between the critical nucleation radius and temperature in the superheated liquid. Considering the nucleation rate distribution in the liquid under the superheated state, it can be seen that, considering the laser ablation target material is aluminum, homogeneous nucleation is different from heterogeneous nucleation in that there is no established nucleation site. Therefore, its nucleation distribution can only be associated with relevant physical quantities in a statistical distribution way, which can be expressed as:

[0183]

[0184] The exponential term e is the critical nucleation free energy mentioned above, which can be used as a whole to represent the equilibrium concentration of critical nuclei; N is the number of liquid molecules per unit volume, and B is an exponential factor used to represent the probability distribution of core growth and collapse caused by vaporization and condensation during the nucleation process. In the Doring-Volmer theory, the B factor is rewritten as two terms, so the nucleation rate per unit volume can be expressed as:

[0185]

[0186] Where σ is the surface tension coefficient, m is the mass of a single molecule, and H v is the molar enthalpy,

[0187] So far, we can get the relationship between the homogeneous nucleation rate and temperature, such as Figure 2 shown.

[0188] Below about 0.9T c At this temperature, the nucleation rate in the liquid is almost zero, and the liquid will still exist in a metastable state and continue to overheat; until the liquid temperature reaches 0.9T c Homogeneous nucleation occurs instantaneously near 0.9T, and the nucleation rate increases rapidly with the increase of temperature after this moment. This phenomenon is mainly due to the fact that the molar enthalpy and the pressure ratio inside and outside the interface increase with temperature at 0.9T. c Strong changes occur near the metastable system, and the metastable system is disturbed, causing an instantaneous increase in the nucleation rate; when the system temperature reaches about 0.95T c, the homogeneous nucleation rate gradually tends to be stable. The reason is that the surface tension coefficient drops rapidly after this moment, the probability of gas nucleus collapse increases, and the probability of effective nucleation decreases, which eventually leads to the nucleation rate being restricted and stabilizing. The nucleation rate "avalanche" decreases in the final stage, mainly because under the critical temperature, most parameters of the target material have reached the thermodynamic end point, and numerical calculations can no longer be performed in the original state. At the same time, in the actual physical process, the liquid phase near the critical temperature will spontaneously and smoothly transition to the gas phase, and all the liquid in the molten pool will completely transform into gas. This large-scale liquid phase nucleation will eventually trigger more violent explosions and ejection processes. These violent evolution processes also determine the fact that the nucleation rate curve cannot reach the thermodynamic end point.

[0189] In this way, the distribution of the mass fraction of the discrete phase with temperature, pressure and material properties in the current system is obtained, and the transmission effect results are given in the interface.

[0190] Furthermore, the distribution of the discrete phase mass fraction as a function of temperature, pressure, and material properties was incorporated into the fluid transport model, and a discrete phase interface was introduced to describe the interaction between the continuous and discrete phases. Finally, the interaction results were used through particle tracking to simulate the gas-liquid splashing process at the moment of the explosion.

[0191] In this stage, the relevant parameters are obtained: the pressure difference between the inside and outside of the phase interface ΔP; the critical nucleation free energy W c ; Nucleation rate per unit volume J.

[0192] In the first embodiment, the two-dimensional ablation model is established. First, according to the task requirements of laser ablation, the key parameters of the laser are set, including:

[0193] Laser energy: E=1,1.5mJ; laser pulse width t p =15ns; laser beam radius r0 = 150um; laser wavelength 1064nm;

[0194] For high peak power laser ablation, due to the existence of multiple physical field coupling processes, the degree of temporal and spatial nonlinearity is high. Therefore, in order to effectively reduce the computational cost, it is necessary to reasonably divide the grid, such as Figure 3 The selected metal target material is aluminum, and the thermophysical parameters of the material are shown in Table 1:

[0195] Table 1 Thermophysical parameters of materials

[0196]

[0197]

[0198] The pulsed laser ablation process involves extremely high pressure, temperature, and velocity gradients, resulting in highly nonlinear effects in both time and space. To reduce the cost of numerical calculations while ensuring accuracy, a two-dimensional model was used to simulate the process.

[0199] In the numerical simulation, the pulsed laser energy is loaded by a body heat source, and its spatiotemporal distribution function is as follows:

[0200]

[0201] Where α is the target's dynamic absorptivity, tp is the laser pulse width, t0 is the initial injection delay, E is the laser energy, and r0 is the laser irradiation radius. The initial temperature of the target boundary and the ambient environment is 20°C.

[0202] Solution of Metal Ablation Aluminum Metal Overheating Model

[0203] The governing equations are as follows:

[0204]

[0205] Among them C p is the dynamic specific heat capacity, which is obtained by analyzing the content of the invention, ρ is the target density given by the material library, k is the target thermal conductivity given by the material library, and p is the pressure. When the temperature of the ablation area exceeds the melting temperature, the thermal conductivity k will become:

[0206] k=ρ0L(1+ΓT)T

[0207] Where ρ0 is the resistivity of the target at room temperature, Γ is the temperature coefficient of metal resistivity, and L is the Lorentz constant.

[0208] Q is the sum of all other heat sources, including the heat removed by evaporation during the gas-liquid phase transition, convective heat transfer, and energy losses due to thermal radiation. The Carman-Kozeny equation is used to describe the target material's melting phase transition. The gaseous phase transition needs to be differentiated based on its phase transition intensity. Theoretical analysis shows that the gasification phase transition process can be divided into two types: evaporation and boiling.

[0209] For normal evaporation, the Hertz-Knudsen equation is used to describe the mass transfer.

[0210]

[0211] β R is the backscattering coefficient, P sat is the saturated vapor pressure of the system, and M is the atomic mass of the metal.

[0212] The energy loss caused by the gas-liquid phase change is expressed as:

[0213]

[0214] Where L v is the latent heat of vaporization. Referring to Watson's theory again, we can get the relationship between phase change enthalpy and temperature:

[0215]

[0216] ΔH in the above formula v0 is the phase change enthalpy required for the gas-liquid transition at atmospheric boiling temperature. This result describes the change in phase change enthalpy at different temperatures. It can be seen that as the temperature approaches the critical temperature, the phase change enthalpy gradually approaches zero. The corresponding enthalpy variation at each phase change temperature is obtained, and thus the latent heat of vaporization is obtained.

[0217] In the case of strong evaporation, there will be a huge impact on the air. Under the impact of high-speed flow, the surrounding gas is compressed, resulting in obvious interface stratification. Therefore, for the problem of strong evaporation under laser action, it is necessary to introduce a gas phase transport model. The model control equation can be expressed as:

[0218]

[0219] The vapor diffusion transport process on the target surface is described by Fick's law:

[0220]

[0221] Where D F is the reduced diffusion factor in Fick's law. The density during the transport process is described by the ideal gas state equation combined with the mass concentration obtained from Fick's law, which can be obtained as follows:

[0222]

[0223] M=ωM v +(1-ω)M air

[0224] In the formula, M represents the average molar mass of the gas system, and the subscripts v and air represent the metal vapor phase and the air phase, respectively. Considering the possibility of small spatial structures during the ablation process, the Knudsen diffusion correction is applied to the diffusion factor in Fick's law, taking into account the collision of gas molecules with the wall surface, to obtain:

[0225]

[0226] Where D K is the diffusion coefficient of the Knudsen diffusion model, and its value is related to the mean free path of the gas:

[0227]

[0228] Where λ path is the mean free path length of the gas molecules, and M is the molar mass of the gas. Calculations using the diffusion equation reveal the actual driving force of the transport process at the gas-liquid interface and the interaction between gas molecules and the surrounding gas after entering the environment.

[0229] For phase transport and tracking problems during gas-liquid phase change, the level set method can be used to capture the phase interface. The governing equation is as follows:

[0230]

[0231] is a level set variable, with 0 and 1 representing the gas and liquid phases, respectively. To prevent sharp phase changes during the calculation process, a phase interface thickness control parameter ε is introduced, and γ is a reinitialization parameter. By studying the mass transfer problem on both sides of the phase interface during evaporation, according to the conservation of mass and momentum, we obtain:

[0232]

[0233] Where ρ v and ρ l Represent the density of gaseous and liquid metal on both sides of the interface, V v and V l The spatial distribution of the phase interface is described by defining delta functions for the vapor phase and the liquid phase respectively:

[0234]

[0235] The last two terms in the formula represent the spatial distribution of the interface and the normal direction of the interface, respectively.

[0236] Bring various parameters, physical quantities and laser heat source into the ablation model, and obtain the laser ablation metal aluminum overheating model by solving the above control equations. Figure 4 The temperature evolution of the target center when the laser energy is 1mJ is shown in Figure 2. Figure 5 The temperature rise curve shows that under the action of nanosecond pulse laser, the temperature rise rate is extremely high; at the same time, when the temperature reaches the normal pressure boiling point of 2793K, the temperature rise curve does not stop or a phase change platform appears, and finally reaches 4750K, which is lower than the critical temperature of the material of 6700K.

[0237] Explosion particle jet modeling

[0238] There are several difficulties in establishing a model for the superheated nucleation stage: the nucleation process and nucleation sites are difficult to predict; the rate of increase in nucleation rate will bring about extremely strong nonlinear effects; and the initial core growth is difficult to achieve through phase interface tracking.

[0239] For the above problems, the corresponding numerical calculation methods adopt the following solutions:

[0240] The prediction of the nucleation process and nucleation sites is carried out by measuring the relationship between the critical nucleation radius, critical nucleation Gibbs free energy, and related thermophysical parameters associated with the nucleation rate and temperature and pressure. The dynamic process of the boiling nucleation rate as a function of the system state is comprehensively obtained, and the spatial distribution of the target nucleation rate is also obtained.

[0241] The highly nonlinear growth of the nucleation rate is one of the characteristics of homogeneous nucleation. This process is bound to lead to a series of situations such as a sharp change in energy transfer at the target interface, obvious changes or even jumps in interface temperature and pressure gradients, and instability of thermophysical parameters. The highly discontinuous distribution of time and space in this process needs to be optimized during the numerical calculation process. In the model, characteristic variables are used to mark the time distribution, and pressure, temperature, and nucleation rate are used as threshold reference points. The "event" capture interface is enabled to capture its sensitive processes, amplify the highly nonlinear process in time, and reduce the computational pressure. On the one hand, the mesh quality is optimized for the spatial distribution, and on the other hand, high-order partial differential equations are used to perform spatial discretization to reduce the spatial gradient distribution of variables and optimize convergence.

[0242] Under the action of nanosecond laser energy of 1.5mJ, the target material reaches the solid-liquid phase transition plateau at about 10ns, during which the temperature rise rate is relatively slow; within 10-20ns, the target surface temperature rises rapidly, and the temperature rise gradient increases significantly, which is mainly due to the changes in thermal physical parameters such as absorptivity, thermal conductivity and specific heat capacity after the solid-liquid phase transition of the target material. After a short temperature peak appears in the target material at 21ns, it is maintained at around 6350K (about 0.95Tc). According to the classical nucleation theory, after the target surface temperature rises to 0.9Tc, the nucleation rate per unit volume increases in a sudden manner. Driven by a large number of gas nuclei, the target material will expand rapidly in a short time and then explode. After a brief equilibrium of 1-2ns, the target surface temperature drops rapidly to about 3000K under the action of the strong current of the explosion products. Figure 6 shown.

[0243] The inherited solution of the front-end superheat model is adopted, and the distribution of the superheated liquid area, temperature, and pressure is used as the initial conditions of this stage. The discrete phase model is used to calculate its nucleation and growth process, and the distribution of the discrete phase mass fraction with temperature, pressure, and material properties in the current system is obtained, and the interface of its transmission effect is given. Furthermore, the initial value is brought into the fluid transport model to obtain the spatiotemporal distribution of the particle explosion products, such as Figure 7 As shown, (A is the time of 21ns, (B is the time of 22ns, (C is the time of 23ns, and (D is the time of 24ns.

[0244] The numerical simulation results are coupled to particle tracking by using the distribution coupling method, and the phase explosion core generation process is simulated by setting the release conditions to obtain the phase explosion space particle distribution results, such as Figure 8 shown.

[0245] The various embodiments in this specification are described in a progressive manner, with each embodiment focusing on the differences from other embodiments. Reference can be made to the common and similar parts between the various embodiments. For the devices disclosed in the embodiments, since they correspond to the methods disclosed in the embodiments, the description is relatively simple, and the relevant parts can be referred to the method description.

[0246] The above description of the disclosed embodiments is intended to enable one skilled in the art to implement or use the present invention. Various modifications to these embodiments will be readily apparent to one skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of the present invention. Therefore, the present invention is not limited to the embodiments shown herein but is intended to conform to the widest scope consistent with the principles and novel features disclosed herein.

Claims

1. A numerical simulation method of laser-induced phase explosion based on finite element method, characterized in that: include: Preset laser parameters and establish a numerical simulation model for nanosecond laser ablation of metals; The laser energy obtained by the numerical simulation model is used as a heat source term to mesh the calculation domain of metal ablation to generate a two-dimensional numerical grid; In the two-dimensional numerical grid, the phase change heat transfer theory is used to consider the phase change of the metal, the phase interface movement and the particle evaporation backscattering effect to simulate the physical process and obtain the saturated vapor pressure, hot liquid pressure and material surface temperature; According to the saturated vapor pressure, hot liquid pressure and material surface temperature, the numerical simulation model is discretized and solved in combination with the finite element method to obtain the spatial and temporal evolution characteristics of the phase explosion particles, the phase change depth and the particle injection characteristics; finally, the spatial distribution characteristics and temporal evolution laws of the temperature and pressure in the metal liquid phase region under laser induction are obtained, and the gas-liquid phase splashing process at the moment of the explosion is simulated through particle tracking.

2. The numerical simulation method of laser-induced phase explosion based on finite element method according to claim 1, characterized in that: The laser parameters include: laser peak power, laser pulse width, laser beam radius, and target material state parameters including metal density, specific heat capacity, and thermal conductivity.

3. The numerical simulation method of laser-induced phase explosion based on finite element method according to claim 1, characterized in that: The physical process in the physical process simulation is divided into: heat conduction temperature rise stage, solid-liquid phase change stage, gas-liquid surface evaporation stage and liquid phase superheated boiling nucleation stage.

4. The method for numerical simulation of laser-induced phase explosion based on finite element method according to claim 3, characterized in that: The heat conduction temperature rise stage calculates the heat flow and heat flux density, including: Heat flow: Q1=AhΔT; Where Q1 represents the heat flow, h represents the surface heat transfer coefficient, A represents the heat transfer area, and ΔT represents the surface temperature difference; Thermal radiation: Where T is the surface temperature of the material, T0 is the ambient temperature, q r is the heat flux density, χ≈5.67×10 -8 W / (m 2 ·K 4 ) is a constant, and ε is the surface emissivity.

5. The method for numerical simulation of laser-induced phase explosion based on finite element method according to claim 3, characterized in that: The heat capacity C of the gas-liquid phase change stage is calculated in the solid-liquid phase change stage p ,include: C p =C p0 +L m D m +L v D v ; Among them, C p0 is a constant heat capacity, L m represents the latent heat of fusion, The overall latent heat of fusion, T m represents the melting point of the target material, T is the surface temperature of the material, ΔT m Indicates the melting transition zone, L v represents the latent heat of vaporization, represents the latent heat of evaporation coefficient, T b Indicates the boiling temperature at normal pressure, ΔT b Indicates the evaporation transition zone.

6. The method for numerical simulation of laser-induced phase explosion based on finite element method according to claim 3, characterized in that: The enthalpy variable ΔH is calculated during the evaporation stage of the gas-liquid surface v , saturated vapor pressure P sat and steam recoil pressure P r , specifically including: Enthalpy change ΔH v , calculated using the Clausius-Clapeyron equation: Where, T b Indicates the boiling temperature at normal pressure, T C is the critical temperature of the target, ΔH v0 Typical temperature T b The phase change enthalpy under , T is the surface temperature of the material; Saturated vapor pressure P sat : Where ΔV is the specific volume, ΔH v Substitute the Clausius-Clapeyron equation and obtain the saturated vapor pressure by iterative integration; Steam recoil pressure P r : Among them, T s is the interface temperature, β R is the backscatter coefficient.

7. The method for numerical simulation of laser-induced phase explosion based on finite element method according to claim 3, characterized in that: The boiling nucleation stage calculates the pressure difference ΔP inside and outside the phase interface, the critical nucleation free energy W c and the nucleation rate per unit volume J, including: Among them, r is the radius of the gas core, σ is the surface tension coefficient, P v represents the internal critical nucleation pressure, P l represents the external pressure, r c represents the critical nucleation radius, m is the mass of a single molecule, H v is the molar enthalpy, T is the surface temperature of the material, k B is the Boltzmann constant, R is the universal gas constant, N is the number of liquid molecules per unit volume, and B is the exponential factor.

Images