Laser-induced phase explosion numerical simulation method based on finite element method

Abstract

The invention discloses a laser induced phase explosion numerical simulation method based on a finite element method, and relates to the technical field of interaction of high-power laser and materials, laser parameters are preset, and a numerical simulation model of nanosecond laser ablated metal is established; laser energy is used as a heat source item, and a two-dimensional numerical grid is generated for a metal ablation computational domain; a phase change heat transfer theory is adopted, phase change of metal, phase interface movement and a particle evaporation backscattering effect are considered, physical process simulation is carried out, and saturated vapor pressure, hot liquid pressure and material surface temperature are obtained; according to the saturated vapor pressure, the hot liquid pressure and the material surface temperature, discretizing and solving the numerical simulation model in combination with a finite element method, and simulating the gas-liquid phase splashing process at the moment of explosion in a particle tracking mode. According to the method, the phase explosion dynamics evolution characteristics in the laser ablation process can be accurately predicted, and important theoretical support and technical guarantee are provided for the fields of laser processing, laser propulsion, laser micro-nano processing and the like.

Description

Technical Field

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

[0002] With the rapid development of laser technology, laser-induced phase explosion technology has been widely applied in fields such as material processing, laser propulsion, and pulsed laser deposition. However, due to the fact that laser-induced phase explosion involves multi-physical field coupling (such as heat conduction, phase change, gasification, particle ejection, etc.), there are still certain limitations in the existing research methods for accurately predicting the explosion evolution characteristics. Traditional experimental techniques are limited by time resolution and spatial capture ability, while a single theoretical model is difficult to comprehensively describe complex phase explosion behaviors.

[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 those skilled in the art. 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 the laser ablation process, especially considering the interaction between laser-induced phase explosion and metal ablation and gasification effects, providing important theoretical support and technical guarantee for fields such as laser processing, laser propulsion, and laser micro-nano processing.

[0005] To achieve the above object, the present invention adopts the following technical solutions:

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

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

[0008] Take the laser energy obtained from the numerical simulation model as a heat source term, perform mesh division on the computational domain of metal ablation to generate a two-dimensional numerical mesh;

[0009] In the two-dimensional numerical mesh, adopt the phase change heat transfer theory, consider the phase change of the metal, the movement of the phase interface, and the particle evaporation backscattering effect, and perform physical process simulation to obtain the saturated vapor pressure, the thermal liquid pressure, and the material surface temperature;

[0010] Based on the saturated vapor pressure, the hot liquid pressure, and the material surface temperature, the numerical simulation model is discretized and solved by combining the finite element method to obtain the spatial and temporal evolution characteristics of the phase explosion particles, the phase change depth, and the particle ejection characteristics; finally, the spatial distribution characteristics and temporal evolution law of the temperature and pressure in the metal liquid phase region under laser induction are obtained, and the gas-liquid phase splashing process during the explosion instant is simulated by means of particle tracking.

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

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

[0013] Preferably, the heat conduction temperature rise stage calculates the heat flux and heat flux density, including:

[0014] Heat flux:

[0015] Q1 = AhΔT;

[0016] In the formula, Q1 represents the heat flux, 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] In the formula, T is the material surface temperature, 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 solid-liquid phase change stage calculates the heat capacity C p in the gas-liquid phase change stage, including:

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

[0022] Among them, C p0 is the constant heat capacity, L m represents the latent heat of fusion, as a whole represents the latent heat of fusion coefficient, T m represents the melting point of the target material, T is the material surface temperature, ΔT mRepresents the melting transition zone, L v Represents the latent heat of vaporization Represents the latent heat of vaporization coefficient, T b Represents the melting point of the target, ΔT b Represents the evaporation transition zone

[0023] Preferably, the enthalpy change ΔH in the gas-liquid surface evaporation stage v , the saturated vapor pressure P sat and the vapor counterpressure P r , specifically including:

[0024] The enthalpy change ΔH v is calculated by the Clausius-Clapeyron equation:

[0025]

[0026] In the formula, T b represents the boiling temperature under atmospheric pressure, T C is the critical temperature of the target, ΔH v0 is the phase change enthalpy at the typical temperature T b , and T is the surface temperature of the material;

[0027] The saturated vapor pressure P sat :

[0028]

[0029] where ΔV is the specific volume. Substitute ΔH v into the Clausius-Clapeyron equation and obtain the saturated vapor pressure through iterative integration;

[0030] The vapor counterpressure P r :

[0031]

[0032] where Ts is the interface temperature, and β R is the backscattering coefficient

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

[0034]

[0035] where r is the radius of the gas nucleus, σ is the surface tension coefficient, P v represents the internal critical nucleation pressure, P l represents the external pressure, and r cdenotes 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.

[0036] Preferably, the discretization and solution of the numerical simulation model according to the saturated vapor pressure, the thermal liquid pressure and the material surface temperature, combined with the finite element method, specifically include:

[0037]

[0038] Among them, is the gradient operator, u is the velocity vector, is the temperature gradient, K is the Darcy resistance coefficient, I is the identity matrix, p is the pressure, is the divergence of the velocity, is the velocity gradient, u x 、u y and u z are the components of the velocity vector u in the x, y, and z directions respectively, C p is the heat capacity, ρ is the density of the target material, k is the thermal conductivity of the target material; Q is the sum of other heat sources, which is composed of the heat carried away by the gas-liquid phase change evaporation, the energy loss caused by convective heat transfer and thermal radiation; α is the dynamic absorption rate of the target material, 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 direction perpendicular to the laser beam, t is the time, g is the acceleration due to gravity, μ is the dynamic viscosity.

[0039] Through the above technical solutions, it can be seen that compared with the prior art, the present invention discloses a numerical simulation method for 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 predicting the characteristics of phase explosion, especially in the accurate 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 technical solutions in the embodiments of the present invention or the prior art, the following will briefly introduce the drawings required for the description of the embodiments or the prior art. Obviously, the drawings in the following description are only the embodiments of the present invention. For those of ordinary skill in the art, other drawings can be obtained according to the provided drawings without creative efforts.

[0041] Figure 1 is the method step diagram provided by the present invention;

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

[0043] Figure 3 The meshing 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 diagram provided by the present invention;

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

[0047] Figure 7 The spatio-temporal distribution diagram of the particle explosion products provided by the present invention;

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

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

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

[0051] Figure 11 The schematic diagram of the P-V-T relationship of the van der Waals equation of state provided by the present invention. Specific embodiments

[0052] Next, the technical solutions in the embodiments of the present invention will be clearly and completely described in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all the embodiments. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present invention without creative efforts shall fall within the protection scope of the present invention.

[0053] The embodiments of the present invention disclose a numerical simulation method for laser-induced phase explosion based on the finite element method, as Figure 1 shown, including:

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

[0055] S2: Take the laser energy obtained from the numerical simulation model as the heat source term, mesh the computational domain of metal ablation, and generate a two-dimensional numerical mesh;

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

[0057] S4: According to the saturated vapor pressure, the pressure of the hot liquid, and the surface temperature of the material, combine with the finite element method to discretize and solve the numerical simulation model, and the phase change depth, and finally obtain the spatial distribution characteristics and time evolution law of the temperature and pressure in the metal liquid phase region under laser induction.

[0058] In a specific embodiment, the numerical simulation model of nanosecond laser ablation of metal considers the spatial distribution of the temperature and pressure in the molten pool in the liquid phase superheat region.

[0059] In a specific embodiment, S4 specifically includes: obtaining the distribution of the superheated liquid region, temperature, and pressure in the numerical simulation model of nanosecond laser ablation of metal, and using the inheritance decoupling method to use the distribution of the superheated liquid region, temperature, and pressure as the initial value conditions at this stage. Calculate the nucleation and growth process using the discrete phase model to obtain the distribution of the discrete phase mass fraction with respect to temperature, pressure, and material physical property parameters in the current system, and give the interface of its transport effect. Further, substitute the distribution of the discrete phase mass fraction with respect to temperature, pressure, and material physical property parameters into the fluid transport model, introduce the discrete phase interface, and describe the interaction process between the continuous phase and the discrete phase. Finally, use the particle tracking method for the action result to simulate the gas-liquid phase splashing process at the moment of explosion, as Figure 10 shown.

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

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

[0062] In a specific embodiment, the heat conduction temperature rise stage calculates the heat flux and heat flux density, including:

[0063] Heat flux:

[0064] Q1 = AhΔT;

[0065] In the formula, Q1 represents the heat flux, 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] In the formula, T is the surface temperature of the material, T0 is the ambient temperature, and 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 change stage is calculated during the solid-liquid phase change stage p , including:

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

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

[0072] In a specific embodiment, the enthalpy change ΔH v , the saturation vapor pressure P sat and the vapor back pressure P r are calculated during the gas-liquid surface evaporation stage. Specifically, it includes:

[0073] The enthalpy change ΔH v is calculated by the Clausius-Clapeyron equation:

[0074]

[0075] In the formula, T b represents the boiling temperature under normal pressure, T C is the critical temperature of the target material, ΔH v0 is the enthalpy of phase change at the typical temperature T b , and T is the surface temperature of the material;

[0076] The saturation vapor pressure P sat :

[0077]

[0078] Among them, ΔH v is substituted into the Clausius-Clapeyron equation, and the saturated vapor pressure is obtained by iterative integration. ΔV is the specific volume.

[0079] Vapor backpressure P r :

[0080]

[0081] Among them, Ts is the interface temperature, and β R is the backscattering coefficient.

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

[0083]

[0084] Among them, r is the radius of the gas nucleus, σ 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, according to the saturated vapor pressure, the hot liquid pressure and the surface temperature of the material, the discretization and solution of the numerical simulation model are carried out by combining the finite element method, specifically including:

[0086]

[0087] Among them, is the gradient operator, u is the velocity vector, is the temperature gradient, K is the Darcy resistance coefficient, I is the identity matrix, p is the pressure, is the divergence of the velocity, is the velocity gradient, u x 、u y and u z are the components of the velocity vector u in the x, y, and z directions respectively, C pwhere \(C_p\) is the heat capacity, \(\rho\) is the density of the target material, \(k\) is the thermal conductivity of the target material; \(Q\) is the sum of other heat sources, which is composed of the heat carried away by the gas-liquid phase change evaporation, the energy loss caused by convective heat transfer and thermal radiation; \(\alpha\) is the dynamic absorption rate of the target material, \(t\) p is the laser pulse width, \(t_0\) is the initial injection delay, \(E\) is the laser energy, \(r_0\) is the laser irradiation radius, \(x\) is the direction perpendicular to the laser beam, \(t\) is the time, \(g\) is the acceleration due to gravity, and \(\mu\) is the dynamic viscosity.

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

[0089] Thermal conduction temperature rise stage

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

[0091]

[0092] In the formula, \(\varepsilon_0\) represents the vacuum permittivity, \(\omega\) represents the angular frequency of the incident light, \(c\) represents the speed of light in vacuum. \(\rho\) is the resistivity of the material, and \(\lambda\) is the wavelength of the incident laser. The absorption rate of the target material for the laser is substituted into the following Gaussian heat source formula.

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

[0094] Thermal conduction:

[0095] According to Fourier's law:

[0096]

[0097] In the formula, \(q\) represents the energy density, and \(k\) is the thermal conductivity of the material.

[0098] Convection:

[0099] According to Newton's cooling formula:

[0100] \(Q_1 = Ah\Delta T\);

[0101] In the formula, \(Q_1\) represents the heat flow rate, \(h\) represents the surface heat transfer coefficient, \(A\) is the heat transfer area, and \(\Delta T\) represents the surface temperature difference.

[0102] Thermal radiation:

[0103] According to the Stefan-Boltzmann law:

[0104]

[0105] In the formula, \(T\) is the surface temperature of the material, \(T_0\) is the ambient temperature, \(q\) r is the heat flux density, \(\sigma\approx5.67\times10-8 W / (m 2 ·K 4 ) is a constant, and ε is the surface emissivity.

[0106] Solid-liquid phase change stage

[0107] After the target material absorbs the laser energy, its temperature rises to the solid-liquid phase change point. According to the theory of phase change heat transfer, 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 the r-axis.

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

[0111]

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

[0113]

[0114] Among them, due to the existence of the phase change transition zone, the problem of latent heat of phase change needs to be processed. Therefore, the specific heat capacity needs to be described by the dynamic specific heat capacity method:

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

[0116] In the formula, C p0 is the constant specific heat capacity, L m represents the latent heat of fusion, as a whole represents the latent heat of fusion coefficient T m represents the melting point of the target material, and ΔT m represents the melting transition zone.

[0117] For the gas-liquid phase change stage, the specific heat capacity will be modified as:

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

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

[0120] Using the above relationships and combining with formulas, the control equation of the solid-liquid interface can be obtained, where the mass conservation and momentum conservation formulas are the fluid transport model:

[0121] Mass conservation:

[0122]

[0123] Momentum conservation:

[0124]

[0125] Energy conservation:

[0126]

[0127] In the formula, T is the 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] In the formula, K is the Carman-Kozeny permeability coefficient, C is the Carman-Kozeny constant, and b is a very small constant to control the denominator not to be zero. ε is the porosity, that is, the liquid volume fraction in the solid-liquid phase change process:

[0130]

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

[0132]

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

[0134] Gas-liquid surface evaporation stage

[0135] 1) Saturated vapor pressure prediction

[0136] A relatively common form in laser ablation is the Clausius-Clapeyron equation:

[0137]

[0138] This invention focuses on simulating the process of laser-induced phase explosion. The conventional reduced formula for saturated vapor pressure will deviate significantly at high temperatures. Therefore, the Watson equation is used to describe the enthalpy change ΔH during the phase change process. v

[0139]

[0140] In the formula, T b represents the boiling temperature under atmospheric pressure, T C is the critical temperature of the target material, and ΔH v0 is the enthalpy of phase change at the typical temperature T b , which is usually taken as the value corresponding to the boiling point under standard atmospheric pressure. Substituting ΔH v into the Clausius-Clapeyron equation, the saturated vapor pressure is obtained through iterative integration.

[0141] 2) Research on the kinetic process of the Knudsen layer

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

[0143]

[0144] In the formula, f( + ) and f (-) represent the forward scattering and backward scattering distribution functions respectively. Vx, Vy, and 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 velocity of the backward scattered particles, k B is the Boltzmann constant, m is the mass of the metal atom, ρ v , and Tv are the vapor pressure and temperature respectively. β R is the backward scattering coefficient. In the case 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 returning, and the net flux is zero. As the evaporation intensity increases, the difference in the interface imbalance distribution gradually becomes larger. Therefore, β R is introduced as the backward scattering coefficient to consider the backward scattering situation in the Knudsen layer during the strong evaporation process.

[0145] According to the research results of Samokhin:

[0146]

[0147] In the formula, γ is the specific heat ratio of the gas, U is the vapor flow velocity, and Mach is the Mach number.

[0148] The pressure distribution of the metal vapor on the gas-liquid interface. This pressure is mainly generated by two parts: 1) the reaction force of particle evaporation motion; 2) the action of backscattered particles on the interface. By introducing the backscattering coefficient, we can obtain:

[0149]

[0150] In the formula, P r is the vapor counterpressure.

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

[0152] The liquid superheat stage

[0153] The superheated state of the liquid means that the liquid is above the boiling point corresponding to the current ambient pressure without boiling. Based on the ideal gas state equation:

[0154] PV = nRT;

[0155] When considering the actual gas molecular volume, introduce the gas molar volume V m and the single molar molecular volume b, and we get:

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

[0157] Further modify the ideal gas state equation by introducing a constant a to represent the interaction between molecules, then:

[0158]

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

[0160]

[0161] In the formula, P V represents the saturated vapor pressure corresponding to the equilibrium position, V g is the pure gas phase volume on the equilibrium line, V lis the pure liquid volume corresponding to the equilibrium point. The Maxwell equal - area method can effectively correct the errors near the critical point, but it is difficult to solve its cubic equation. Therefore, the non - root - finding method obtained by Carl W. David using intermediate variables is adopted here. By substituting the pressure at the liquid end of the isotherm and the liquid - vapor conversion pressure point equivalently, the relationship between T and V g 、V l can be obtained as follows:

[0162]

[0163] By further introducing the intermediate variable e d , finally, we can get:

[0164]

[0165] Where:

[0166]

[0167] Finally, the relationship among Tr, Vr, Pr, which is the simplified form of the van der Waals equation of state about T / Tc, V / Vc, P / Pc, can be obtained. As shown in Figure 11 , 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 A - B - D - C is the 0.9T c isotherm. The liquid - phase and gas - phase saturation lines converge smoothly at the critical point. This phenomenon is consistent with the fact that the enthalpy of phase change becomes zero at the critical point described in the previous text. Therefore, the pure gas volume and the pure liquid volume corresponding to the equilibrium point on the equilibrium line in the liquid - phase superheat stage are not used in the simulation of the present invention.

[0168] To further illustrate the superheat problem, the coordinates of points D' and D are extracted from the curve in the figure. It can be seen that the pressures at the two points are equal, and there are differences in temperature and gas volume. At the same time, since D' is a point on the saturated liquid curve, point D can be regarded as the state reached by the saturated liquid system through isobaric heating. Correspondingly, it can be found that point D has not transitioned to the saturated vapor curve at this time. Therefore, point D shows liquid superheat at this time and is in a "metastable state". After that, point D may spontaneously transition along the isotherm to C and return to the gas - phase stable point; it can also reach the new liquid - phase stable point B by external pressure supply.

[0169] From the calculation and study of the superheat problem, it can be seen that superheat can exist at any position before reaching the critical point, and the reason is the lack of nucleation sites. Usually, when the system reaches the critical point, the gas - liquid phase will undergo spontaneous transition. This conclusion is valid for subsequent thermodynamic parameters at T cThe description of the point provides strong evidence; on the other hand, due to the existence of metastable states, during the laser ablation process, as the laser power increases, the overheating problem will become particularly important. It has a crucial impact especially in aspects such as determining the laser ablation threshold and studying the ablation efficiency; finally, the instability of the liquid under metastable states will also provide leading calculation results for research on homogeneous nucleation, explosive boiling, etc.

[0170] The situation of the hot liquid region, temperature, and pressure distribution is obtained therefrom, and the situation of the superheated liquid region, temperature, and pressure distribution is used as the initial value 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 calculation. The formulas included in this stage are the discrete phase model, and the situation of the hot liquid region, temperature, and pressure distribution is used as the initial value of the discrete phase model. According to the equilibrium state theory, the equilibrium temperature of the system depends on the pressure and volume changes. Described by the Young-Laplace equation:

[0173]

[0174] In the formula, ΔP is the pressure difference inside and outside the phase interface (external pressure minus internal pressure), r is the radius of the gas nucleus, and σ is the surface tension coefficient. This formula relates the surface tension to the pressure difference and is used to describe the work required to be done externally during the formation and growth of the gas phase core. To further study the nucleation problem during the phase explosion process, according to the classical nucleation theory:

[0175]

[0176] In the formula, r c represents the critical nucleation radius, and W c represents the critical nucleation free energy. Considering the nucleation process of superheated liquid under non-equilibrium states and introducing the Young-Laplace relationship, we can obtain:

[0177]

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

[0179]

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

[0181]

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

[0183]

[0184] The exponential term e in the formula 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 cWhen the homogeneous nucleation rate gradually stabilizes, the reason is that the surface tension coefficient drops rapidly after this moment, the probability of gas nucleus collapse increases, and the effective nucleation probability decreases, ultimately leading to the limited increase and stabilization of the nucleation rate. In the final stage, an "avalanche" drop in the nucleation rate occurs mainly because at the critical temperature, most of the parameters of the target material reach the thermodynamic end point and cannot participate in the operation in the original state in numerical calculations. At the same time, in the actual physical process, near the critical temperature, the liquid phase will spontaneously and smoothly transition to the gas phase, and all the liquid in the molten pool will be completely transformed into gas. This large-scale liquid-phase nucleation will ultimately trigger a more violent explosion and ejection process, and these violent evolution processes also determine the fact that the nucleation rate curve cannot reach the thermodynamic end point.

[0189] Thus, the distribution of the discrete phase mass fraction with respect to temperature, pressure, and material physical property parameters under the current system is obtained, and the result of its transport effect is given as an interface.

[0190] Furthermore, the distribution of the discrete phase mass fraction with respect to temperature, pressure, and material physical property parameters is introduced into the fluid transport model, and a discrete phase interface is introduced to describe the interaction process between the continuous phase and the discrete phase. Finally, the result of the interaction is used to simulate the gas-liquid splashing process at the moment of explosion through particle tracking.

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

[0192] In Specific Example 1, for the establishment of the two-dimensional ablation model, first, according to the task requirements of laser ablation, the key parameters of the laser are set, specifically including:

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

[0194] For high peak power laser ablation, due to the existence of multiple physical field coupling processes and a high degree of time and space non-linearity, in terms of grid meshing, in order to effectively reduce the calculation cost, the grid needs to be reasonably divided, as Figure 3 shown. The selected metal target is aluminum, and the thermophysical parameters of the material are shown in Table 1:

[0195] Table 1 Thermophysical parameters of the material

[0196]

[0197]

[0198] During the pulsed laser ablation calculation process, there are extremely high pressure gradients, temperature gradients, and velocity gradient changes. Therefore, the model has high nonlinear effects both in time and space. To reduce the cost required for numerical calculation while ensuring the calculation accuracy, a two-dimensional model is used to numerically simulate the process.

[0199] In the numerical simulation, the pulsed laser energy is loaded through a volumetric heat source, and its spatio-temporal distribution function is as follows:

[0200]

[0201] In the formula, α is the dynamic absorption rate of the target material, 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 temperatures of the target material boundary and the environment are both 20 °C.

[0202] Solution of the 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, obtained from the analysis of the invention content, ρ is the target material density given by the material library, k is the target material thermal conductivity given by the material library, p is the pressure. When the temperature in the ablation area exceeds the melting temperature, the thermal conductivity k will become:

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

[0207] In the formula, ρ0 is the resistivity of the target material under normal temperature environment, Γ is the metal resistivity temperature coefficient, and L is the Lorentz constant.

[0208] Q is the sum of other heat sources, composed of the heat taken away by the gas-liquid phase change evaporation, the energy loss caused by convective heat transfer and thermal radiation. For the melting phase change of the target material, the Carman-Kozeny equation is used to describe it. For the gaseous change part, it needs to be distinguished according to its phase change intensity. According to theoretical analysis, the gasification phase change process can be divided into two situations: evaporation and boiling.

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

[0210]

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

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

[0213]

[0214] where L v is the latent heat of evaporation. By citing Watson's theory again, the relationship between the phase change enthalpy and temperature can be obtained as follows:

[0215]

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

[0217] For strong evaporation conditions, it will have a great impact on the air. Under the impact of the high-speed flow, the surrounding ambient gas is compressed, resulting in an obvious interface stratification phenomenon. Therefore, for the problem of strong evaporation under laser action, a gas-phase transport model needs to be introduced. The governing equations of the model can be expressed as:

[0218]

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

[0220]

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

[0222]

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

[0224] where 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. At the same time, considering that there may be relatively small spatial structures during the ablation morphology process, the collision of gas molecules with the wall surface is considered, and the diffusion factor in Fick's law is corrected by Knudsen diffusion 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] In the formula, λ path is the mean free path length of gas molecules, and M is the molar mass of the current gas. Through the calculation of the diffusion equation, the actual driving situation of the gas-liquid interface transport process can be obtained; on the other hand, the interaction between gas-phase molecules and environmental gas after entering the environment can be obtained.

[0229] For the phase transport and tracking problems in the gas-liquid phase change process, the level set method can be used to capture the phase interface, and the control equations are as follows:

[0230]

[0231] is the level set variable, representing the gas phase and liquid phase with 0 and 1 respectively. To prevent sharp phase mutations in the calculation process, a phase interface thickness control parameter ε is introduced, and γ is the re-initialization parameter. Through the study of the mass transfer problem on both sides of the phase interface during the evaporation process. According to the mass conservation and momentum conservation relationships, we get:

[0232]

[0233] In the formula, ρ v and ρ l represent the gaseous and liquid metal densities on both sides of the phase interface respectively, V v and V l are the defined δ functions of the vapor phase and liquid phase respectively to describe the spatial distribution of the phase interface:

[0234]

[0235] In the formula, the latter two terms represent the interface spatial distribution and the interface normal direction respectively.

[0236] Substitute each parameter, physical quantity, and the laser heat source into the ablation model, and obtain the laser ablation aluminum overheating model by solving the above control equations. Figure 4 is the evolution result of the central temperature of the target irradiated by 1 mJ of laser energy. Through Figure 5 the temperature rise curve, it can be found that under the action of nanosecond pulsed laser, the temperature rise rate is extremely high; at the same time, when the temperature reaches the normal boiling point of 2793 K, the temperature rise curve does not stop or show a phase change plateau, and finally reaches 4750 K, which is less than the critical temperature of the material of 6700 K.

[0237] Modeling of explosive particle ejection

[0238] For the model establishment of the superheat nucleation stage, there are the following difficulties: the nucleation process and nucleation sites are difficult to predict; the rising rate of the nucleation rate will bring a very strong nonlinear effect; the initial core growth is difficult to achieve by the phase interface tracking method.

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

[0240] For the prediction of the nucleation process and nucleation sites, it is measured by the critical nucleation radius, critical nucleation Gibbs free energy associated with the nucleation rate, and the relationships between relevant thermal physical properties and temperature and pressure. The dynamic process of the boiling nucleation rate with respect to the system state is comprehensively obtained, and at the same time, the spatial distribution of the nucleation rate of the target material is obtained;

[0241] The highly non - linear growth of the nucleation rate is one of the characteristics of homogeneous nucleation. This process will inevitably lead to a sharp change in the energy transfer at the target interface, obvious changes or even jumps in the interface temperature and pressure gradients, and unstable thermal physical properties. The highly discontinuous distribution of this process in time and space needs to be optimized during the numerical calculation. In the model, the time distribution is marked with characteristic variables, and the pressure, temperature, and nucleation rate are used as threshold reference points. The "event" capture interface is enabled to capture its sensitive processes, magnifying the highly non - linear process in time to reduce the computational pressure. For the spatial distribution, on the one hand, the grid quality is optimized, and on the other hand, high - order partial differential equations are used for spatial discretization to reduce the spatial gradient distribution of variables and optimize the convergence;

[0242] Under the action of a nanosecond laser energy of 1.5 mJ, the target material reaches the solid - liquid phase transition plateau at about 10 ns. During this period, the temperature rise rate is relatively slow; within the time range of 10 - 20 ns, the temperature on the target surface rises rapidly, and the temperature rise gradient increases significantly. This is mainly due to the changes in thermal physical properties such as the absorption rate, thermal conductivity, and specific heat capacity after the solid - liquid phase transition of the target material. After a short - lived temperature peak appears at 21 ns, the target material is maintained near 6350 K (about 0.95 Tc). According to the classical nucleation theory, after the target surface temperature rises to 0.9 Tc, the nucleation rate per unit volume undergoes a jump - type growth. Driven by a large number of gas nuclei, the target material will expand rapidly in a short time and then undergo a phase explosion. After a short - lived equilibrium at 1 - 2 ns on the target surface, the temperature rapidly drops to about 3000 K under the action of the strong flow of the explosion products, as Figure 6 shown.

[0243] The front - end superheat model inheritance solution is adopted, and the superheated liquid region, temperature, and pressure distribution are used as the initial value conditions for this stage. The discrete phase model is used to calculate the nucleation and growth processes, and the distribution of the discrete phase mass fraction with respect to temperature, pressure, and material physical properties in the current system is obtained, and the results of its transport effect are given at the interface. Further, the initial values are introduced into the fluid transport model to obtain the spatio - temporal distribution of the particle explosion products, as Figure 7 shown, where (A) is at 21 ns, (B) is at 22 ns, (C) is at 23 ns, and (D) is at 24 ns.

[0244] By adopting the distributed coupling method, the numerical simulation results are coupled to particle tracking, and the generation process of the phase explosion core is simulated by setting the release conditions, and the spatial particle distribution results of the phase explosion are obtained, as Figure 8 shown.

[0245] The various embodiments in this specification are described in a progressive manner. The key point of each embodiment is to illustrate the differences from other embodiments. For the same or similar parts among the various embodiments, reference can be made to each other. For the devices disclosed in the embodiments, since they correspond to the methods disclosed in the embodiments, the description is relatively simple, and reference can be made to the description in the method part for relevant parts.

[0246] The above description of the disclosed embodiments enables those skilled in the art to implement or use the present invention. Various modifications to these embodiments will be obvious to those skilled in the art, and the general principles defined herein can be implemented in other embodiments without departing from the spirit or scope of the present invention. Therefore, the present invention will not be limited to the embodiments shown herein, but will be accorded 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 adopted to consider the phase change of the metal, the movement of the phase interface and the backscattering effect of the particle evaporation, and the physical process simulation is performed to obtain the saturated vapor pressure, the hot liquid pressure and the surface temperature of the material; 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 method for numerical simulation 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 include metal density, specific heat capacity, and thermal conductivity.

3. The method for numerical simulation 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; In the formula, Q1 represents heat flow, h represents surface heat transfer coefficient, A represents heat exchange area, and ΔT represents surface temperature difference; Heat radiation: In the formula, T is the material surface temperature, T0 is the ambient temperature, q r is the heat flux, χ≈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 solid-liquid phase change stage calculates the heat capacity C of the gas-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 evaporation, represents the latent heat of evaporation, 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 change ΔH is calculated during the evaporation stage of the gas-liquid surface v , saturated vapor pressure P sat With steam back pressure P r , specifically including: Enthalpy change ΔH v , calculated using the Clausius-Clapeyron equation: Where, T b represents the boiling temperature at normal pressure, T C is the critical temperature of the target, ΔH v0 The 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 and 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, and R is the universal gas constant.

8. The method for numerical simulation of laser-induced phase explosion based on finite element method according to claim 1, characterized in that: The discrete 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: in, is the gradient operator, u is the velocity vector, temperature gradient, K is Darcy's drag coefficient, I is the unit matrix, p is the pressure, is the velocity divergence, 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 thermal conductivity of the target; 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 gravitational acceleration, and μ is the dynamic viscosity.

Images