A multi-physical field coupling plasma recoil pressure prediction method

Abstract

The application relates to the technical field of laser micro-nano manufacturing, and discloses a multi-physical field coupling plasma recoil pressure prediction method, a three-dimensional finite element calculation domain model containing metal nanoparticles and a substrate is constructed; electromagnetic wave time domain solving is carried out; amplified local light field energy density distribution data are extracted as volume heat sources to be substituted into a double-temperature model; material instantaneous gasification is calculated and spatial temperature distribution and atomic vapor parameters are output; a composite ionization rate equation is input; and high-temperature plasma number density distribution is obtained; based on fluid control equation simulation of supersonic expansion process of the plasma, expansion parameters are converted into surface load data corresponding to recoil force, and the surface load data are applied to a solid mechanics model of the substrate surface; the propagation and superposition of internal stress waves are solved; and finally, dynamic pressure prediction values are calculated. The application breaks through the whole-link physical evolution process of "light-heat-flow-solid", and realizes accurate simulation and prediction of super-high dynamic pressure under cross-scale and multi-physical field coupling.

Description

Technical Field

[0001] This invention relates to the field of laser micro / nano manufacturing technology, and in particular to a method for predicting plasma recoil pressure using multi-physics coupling. Background Technology

[0002] With the rapid development of ultrafast laser technology, the interaction between ultrashort pulse lasers (such as femtosecond lasers) and matter induces extreme high-temperature and high-pressure states, which have significant application value in cutting-edge fields such as materials processing, high-energy physics, and inertial confinement fusion. In recent years, by utilizing solid targets with metal nanoparticles attached to their surfaces and exciting local surface plasmon resonance effects, it is possible to break through the conventional optical diffraction limit and achieve extremely high-rate amplification and efficient absorption of laser energy at the nanoscale, thereby generating ultra-high intensity dynamic pressure and shock waves.

[0003] However, the process of generating dynamic pressure in metal nanostructures induced by ultrashort pulse lasers not only involves cross-scale spatial evolution but also includes an extremely complex multi-physics dynamic coupling mechanism. This physical process, in its temporal evolution, sequentially undergoes: localized deposition and amplification of electromagnetic wave energy; non-equilibrium heat transfer between electrons and the lattice inducing instantaneous vaporization of the material; recombination ionization under the influence of a strong laser field, leading to the generation of high-temperature plasma; and the rapid expansion of the plasma into the vacuum, which, through recoil, excites longitudinal compressive stress waves within the target matrix.

[0004] Current theoretical models and simulation techniques are often limited to calculating single or local physical fields. For example, some existing technologies can only focus on calculating the enhanced light field distribution near nanostructures using electromagnetic field models, or can only simulate the temperature rise and phase transition processes of materials based on traditional thermodynamic models. These fragmented simulation methods cannot seamlessly connect and couple the multiple physical mechanisms of "light-heat-fluid-solid". Due to the lack of a complete cross-scale multi-physics coupling prediction model from the evolution of the initial light field to the final mechanical pressure, those skilled in the art find it difficult to accurately calculate the evolution and ejection characteristics of plasma, let alone accurately predict the peak magnitude and spatial distribution of the dynamic pressure formed inside the matrix. This gap between theoretical models and simulation tools severely restricts the mechanism revelation and engineering application design of plasmonic-enhanced high-pressure laser-driven technology. Summary of the Invention

[0005] This invention provides a multi-physics field coupled plasma recoil pressure prediction method. By irradiating gold nanoparticles with a femtosecond laser, localized surface plasmon resonance (LSPR) is excited, forming localized light field "hot spots" at the tips and gaps of the nanoparticles, achieving a significant enhancement of the incident light intensity by 10-1000 times. The enhanced light field is efficiently absorbed by the gold nanoparticles, and energy is rapidly deposited within the femtosecond scale, causing the gold nanoparticles to rapidly increase in temperature and instantly vaporize. The vaporized gold atoms undergo further strong field ionization, forming a high-temperature, high-density plasma. The plasma expands supersonically towards the vacuum, generating a strong recoil force according to the principle of momentum conservation, thereby exciting stress waves internally and forming a huge dynamic pressure on the order of GPa. This method overcomes the defects of high energy reflection and transmission loss and low energy conversion efficiency when directly irradiating metal targets with femtosecond lasers, while avoiding the side effects of superheated electrons and X-ray preheating of the target introduced by simply increasing the laser power density. It achieves an order-of-magnitude increase in recoil pressure under the same laser input power.

[0006] This invention provides a method for predicting plasma recoil pressure using multi-physics coupling, comprising:

[0007] S1. Establish a three-dimensional finite element computational domain model that includes the substrate and the target geometric structure located on the surface of the substrate;

[0008] S2. Obtain the initial electromagnetic wave parameters of the pulsed laser, and perform electromagnetic wave time-domain solution on the three-dimensional finite element calculation domain model based on the initial electromagnetic wave parameters to extract the local light field energy density distribution data of the target geometric surface.

[0009] S3. The local light field energy density distribution data is used as a volume heat source term and solved based on the dual-temperature model to output the spatial temperature distribution data and atomic vapor parameters of the target geometric structure.

[0010] S4. Input the space temperature distribution data and the atomic vapor parameters into the preset ionization rate equation for solution to obtain plasma number density distribution data;

[0011] S5. Input the plasma number density distribution data into the preset fluid control equation for solution to obtain the plasma expansion velocity and mass ejection rate, and convert the plasma expansion velocity and mass ejection rate into surface load data according to the principle of momentum conservation.

[0012] S6. Apply the surface load data to the substrate surface in the three-dimensional finite element calculation domain model, and perform transient solution based on the solid mechanics equations to obtain the predicted dynamic pressure value in the substrate.

[0013] Furthermore, S1 specifically includes:

[0014] S101. Establish a cuboid solid domain as the base, and construct a unit cell model including spherical solids as the target geometry on the surface of the base based on periodic boundary conditions. Set at least two spherical solids in the unit cell model to form gap spaces between adjacent spherical solids. Establish a fluid expansion domain above the target geometry, and wrap an electromagnetic wave absorbing layer around the base, the target geometry, and the fluid expansion domain to reconstruct a multi-layer geometric computational domain.

[0015] S102. Spatial discretization is performed on the multi-layer geometric computational domain. Different sizes of meshes are used to perform partitioned meshing on the target geometric structure and the gap space, the matrix and the fluid expansion domain, and the electromagnetic wave absorption layer. After the partitioning is completed, the mesh distortion rate is checked to optimize the distorted mesh elements, thus completing the establishment of the three-dimensional finite element computational domain model.

[0016] Furthermore, S2 specifically includes:

[0017] S201. Assign material and electromagnetic properties to the three-dimensional finite element computational domain model, and set pulsed laser source and boundary conditions in the electromagnetic wave time domain solution architecture.

[0018] S202. Set the time-domain solver parameters and run the solution. After the solution is completed, extract the local electric field intensity from the tip and gap space of the target geometry, calculate the local electric field enhancement factor, and output the local light field energy density distribution data.

[0019] Furthermore, in step S202, the formula for calculating the local electric field enhancement factor is:

[0020]

[0021] Where F is the local electric field enhancement factor, and its value ranges from 10 to 1000. For the local electric field strength, The input is the incident laser electric field intensity.

[0022] Furthermore, S3 specifically includes:

[0023] S301. Define electronic temperature and lattice temperature as solution variables, set the parameters of the dual-temperature model, and load the local light field energy density distribution data as the volume heat source term into the region where the target geometry is located; wherein, the parameters of the dual-temperature model include electron-phonon coupling coefficient, electronic heat capacity, and lattice heat capacity;

[0024] S302. Set a segmented time step and perform transient solution for the evolution process of the electron temperature and the lattice temperature based on the dual-temperature model.

[0025] S303. Extract the lattice temperature evolution curve of the target geometry, and when it is determined that the lattice temperature rises above the boiling point, extract and output the temperature spatial distribution result and total mass parameter of the vaporization region, and use the temperature spatial distribution result and the total mass parameter as the spatial temperature distribution data and the atomic vapor parameter, respectively.

[0026] Furthermore, in step S302, the dual-temperature model satisfies the following equation during the solution process:

[0027]

[0028]

[0029] in, For the electronic heat capacity, The lattice heat capacity, The electron temperature, The lattice temperature is given, g is the electron-phonon coupling coefficient, S(t) is the bulk heat source term, and t is time.

[0030] Furthermore, S4 specifically includes:

[0031] S401. Define a preset ionization rate equation in the partial differential equation interface and set the composite model parameters of the ionization rate equation; wherein, the ionization rate equation adopts a composite model of multiphoton ionization and tunneling ionization, and the composite model parameters include the multiphoton ionization order and the tunneling ionization field strength threshold.

[0032] S402. Input the spatial temperature distribution data and the atomic vapor parameters as initial conditions into the preset ionization rate equation to perform transient solution of the strong field ionization process;

[0033] S403. Extract the obtained electron number density and ion number density and perform validity verification. After the verification is passed, output the result containing the electron number density and the ion number density as the plasma number density distribution data.

[0034] Furthermore, S5 specifically includes:

[0035] S501. Set the initial and boundary conditions for simulating the supersonic expansion of plasma in the preset fluid control equation; wherein, the preset fluid control equation adopts the Euler equation for compressible fluid flow, and the initial conditions include the initial plasma temperature and the initial plasma pressure.

[0036] S502. Set time-domain solution parameters to perform transient solution, extract the plasma velocity field and mass flux distribution, and perform surface integration operation on the substrate surface to extract the plasma expansion velocity and mass ejection rate.

[0037] S503. Calculate the total recoil force acting on the substrate surface according to the momentum conservation formula, and convert the total recoil force into the surface load data output.

[0038] Furthermore, S6 specifically includes:

[0039] S601. Input the basic mechanical parameters of the matrix in the solid mechanics interface corresponding to the matrix, apply the surface load data to the surface of the matrix in the form of surface load, and set the load time history and solution boundary conditions.

[0040] S602. Set the time-domain solver parameters and perform transient solution of the propagation and superposition process of longitudinal compressive stress wave inside the matrix based on the solid mechanics equations.

[0041] S603. Extract the compressive stress distribution inside the matrix, calculate the dynamic pressure peak value by combining it with the effective area of ​​laser irradiation, and output the data containing the dynamic pressure peak value and its superimposed distribution result as the predicted value of dynamic pressure in the matrix.

[0042] Furthermore, in S603, the formula for calculating the peak dynamic pressure is:

[0043]

[0044] Wherein, P is the peak dynamic pressure, F is the corresponding total recoil force, A is the effective area of ​​laser irradiation, and v is the plasma expansion velocity. This refers to the mass spray rate.

[0045] The beneficial effects of this invention are as follows:

[0046] 1. Significant improvement in optical field enhancement and energy conversion efficiency. When traditional femtosecond lasers directly irradiate targets, the reflectivity of the metal surface to 800nm ​​wavelength lasers is as high as 70%-90%, while there is a transmission loss of 5%-15%. The effective energy actually deposited on the target surface and used for plasma generation is less than 15%. This invention precisely matches the femtosecond laser photon energy with the LSPR peak of gold nanoparticles, forming local optical field "hot spots" at the tips of the gold nanoparticles and in the interparticle gaps (spacing ≤10nm). The local electric field enhancement factor can reach 10-1000 times. The absorption efficiency of the gold nanoparticles to the incident laser is increased to over 85%, and the energy is deposited within a femtosecond scale of 30-100fs with no significant thermal diffusion loss. The overall laser energy utilization rate is significantly improved compared to traditional femtosecond laser direct loading technology.

[0047] 2. The plasma recoil pressure achieves an order-of-magnitude leap, breaking through the GPa level of extreme high pressure threshold. Currently, when nanosecond lasers directly irradiate pure solid targets, the plasma recoil pressure is only in the range of kilopascals to megapascals (e.g., the maximum recoil pressure measured in the current technology of nanosecond laser ablation of bronze diamond grinding wheels is approximately 1203 Pa); the highest recoil pressure achieved by traditional femtosecond laser direct loading is only 0.1-1 GPa, which cannot meet the needs of extreme high-pressure research. This invention, under the same laser input power, can excite electrons with an enhanced local optical field that reaches a temperature of 10... 4 -10 6 K, ion number density 10 20 -10 22 cm -3 High-temperature, high-density plasma, with plasma ejection velocity into the vacuum increased to 10 6 -10 7 m / s. Based on the principle of conservation of momentum, the recoil pressure on the target surface is greater than 1 GPa, representing an order-of-magnitude improvement over traditional femtosecond laser direct loading technology.

[0048] 3. Effectively avoids the harmful side effects of high-power lasers, significantly improving the accuracy and repeatability of high-pressure loading. Traditional techniques require increasing the laser power density to 10 to increase the recoil pressure. 14 W / cm 2 The above process generates a large amount of high-energy superthermal electrons and hard X-rays, leading to premature preheating of the target material, severe distortion of the stress wave waveform, pressure measurement error exceeding 30%, and poor experimental repeatability. This invention addresses this issue with a laser power density of only 10... 11 -10 13 Under the condition of W / cm², GPa-level high pressure loading can be achieved, and the preheating effect of the target material can be ignored; the waveform fidelity of the stress wave generated by loading is ≥95%, which greatly improves the reliability of extreme high pressure experimental data. Attached Figure Description

[0049] Figure 1 This is a schematic diagram of the generation of localized plasmon-assisted femtosecond laser plasma recoil pressure in this invention.

[0050] Figure 2 This is a diagram illustrating the mechanism of femtosecond laser energy deposition and plasma formation of gold nanoparticles in this invention.

[0051] Figure 3 This is a schematic diagram of the three-dimensional mesh division of the nanoparticle computational domain in this invention.

[0052] Figure 4 The image shows the simulation results of LSPR light field enhancement of gold nanoparticles on a gold thin film substrate in this invention.

[0053] Figure 5 This is a graph showing the variation of plasma recoil pressure of Au and Ag nanoparticles with particle spacing in this invention.

[0054] The realization of the objective, functional features and advantages of the present invention will be further explained in conjunction with the embodiments and with reference to the accompanying drawings. Detailed Implementation

[0055] It should be understood that the specific embodiments described herein are merely illustrative of the invention and are not intended to limit the invention.

[0056] like Figure 1 As shown, this invention provides a method for predicting plasma recoil pressure using multi-physics coupling; simultaneously, as... Figure 2 As shown, this method encompasses the physical mechanism of femtosecond laser energy deposition of gold nanoparticles and the formation of high-temperature, high-density plasma. Specific steps include:

[0057] S1. Establish a three-dimensional finite element computational domain model that includes the matrix and the target geometry located on the surface of the matrix.

[0058] In one embodiment of the present invention, step S1, which constructs a mathematical topological model and a discretized space that conform to the boundary conditions of real physical evolution, specifically includes the following sub-steps:

[0059] S101, Reconstruction of Multi-Layer Geometric Computation Domain

[0060] A cuboid solid domain with a length of 200nm, a width of 200nm, and a thickness of 50nm was established as the underlying substrate model. It was set as a gold thin film substrate material and served as the underlying physical space for receiving subsequent recoil pressure and bearing the propagation of stress waves.

[0061] On the upper surface of the underlying substrate model, a target geometric structure for exciting local surface plasmon resonance is constructed. To accurately reflect the macroscopic distribution of the self-assembled system and reduce computational overhead, a unit cell model of the nanoparticle array is constructed using periodic boundary conditions. Specifically, two spherical entities with an average particle size of 20 nm are configured within the unit cell model as gold nanoparticle models. To accurately capture the local light field enhancement effect, the center-to-center distance between the two spherical entities is constrained to the range of 22 nm to 30 nm, thereby constructing a gap space (i.e., a preset light field hotspot region) with a size of 0 to 10 nm between the particles. Furthermore, by adjusting the number and arrangement characteristics of the spherical entities within the unit cell, the surface coverage parameter is preset.

[0062] Above the target geometry, a fluid expansion domain with a thickness of 100 nm is established. During the electromagnetic wave solution stage, this region is given environmental parameters equivalent to a vacuum (refractive index set to 1), and serves as the initial buffer space for the supersonic expansion of plasma in the subsequent fluid dynamics solution stage.

[0063] An electromagnetic wave absorption layer (such as a perfectly matched layer) with a thickness of 20 nm is wrapped around the outermost edge of the entire geometric computation domain to absorb electromagnetic waves emitted outside the computation domain, thereby eliminating the interference of non-physical reflections at the cut-off boundary on the accuracy of optical field calculation.

[0064] S102. After reconstructing the multi-layer geometric computational domain, the multi-layer geometric computational domain is spatially discretized. Different sized meshes are used to perform partitioned meshing on the target geometric structure and the gap space, the matrix and the fluid expansion domain, and the electromagnetic wave absorption layer, respectively. To achieve the best balance between finite element computational resource consumption and final solution accuracy, the specific strategy for partitioned meshing is as follows:

[0065] ① Target geometry and gap space: For the target geometry (i.e., the surface of gold nanoparticles) and the gap space (i.e., the preset area of ​​light field hotspots) which are regions with drastic changes in local electric field gradient, high-density meshing is performed using free tetrahedral meshes, and the maximum mesh cell size in this region is strictly controlled to not exceed 2nm to ensure the accuracy and fidelity of the numerical calculation of the local electric field enhancement factor.

[0066] ② Substrate and fluid expansion domain: For the substrate (i.e., gold thin film substrate) where the internal field strength or physical gradient distribution changes relatively gently, and the fluid expansion domain, conventional free tetrahedral meshes with a size larger than the aforementioned hot spot region are used for partitioning, and the maximum mesh unit size is controlled to be no more than 10nm.

[0067] ③ Electromagnetic wave absorbing layer: For the outermost electromagnetic wave absorbing layer, a swept mesh partitioning technique is used to divide it into 5 mesh units uniformly and continuously along its thickness direction.

[0068] After all the above-mentioned partitioned meshing is completed, a mesh distortion rate check is performed to optimize distorted mesh cells. Specifically, the average quality of the mesh cells is used as the check index. When a cell with an average quality lower than 0.3 is detected, it is identified as a distorted mesh cell, and automatic or manual removal or mesh reconstruction optimization is performed until the average quality of all mesh cells is better than 0.3. This effectively avoids the impact of numerical divergence and distortion on the final calculation results. The specific three-dimensional meshing of the nanoparticle computational domain is as follows: Figure 3 As shown, the three-dimensional finite element computational domain model is finally established.

[0069] S2. Obtain the initial electromagnetic wave parameters of the pulsed laser, and perform electromagnetic wave time-domain solution on the three-dimensional finite element computational domain model based on the initial electromagnetic wave parameters to extract the local optical field energy density distribution data of the target geometric surface.

[0070] In one embodiment of the present invention, step S2 specifically includes the following sub-steps:

[0071] S201. Assign material and electromagnetic properties to the three-dimensional finite element computational domain model, and set the pulsed laser source and boundary conditions in the electromagnetic wave time-domain solution framework. Specifically, this includes:

[0072] The material properties of the target geometry (i.e., gold nanoparticles) and the substrate (i.e., gold thin film substrate) are set to gold, and their optical parameters are obtained by calling the Johnson-Christy dispersion model to automatically match the complex permittivity at the center wavelength of 800 nm.

[0073] The refractive index of the fluid expansion region (i.e., the region above the nanoparticles) is set to n=1 to create an equivalent vacuum environment.

[0074] A plane wave incident boundary perpendicular to the surface of the substrate and pointing downwards (i.e., along the negative z-axis) is established to introduce a Gaussian femtosecond pulse laser source.

[0075] The initial electromagnetic wave parameters for obtaining the pulsed laser include: a center wavelength of 800 nm, a default pulse half-width at half-maximum (FWHM) of 50 fs (adjustable within the range of 30-100 fs), a laser system repetition frequency of 1-1000 Hz, and a single-pulse loading method. The amplitude of the incident electric field is calculated using the laser power density, with an initial power density set to 10. 12 W / cm² (and can be 10) 11 ~10 13(Adjusted within the W / cm² range). During the solution process, the center wavelength of the laser is adjusted to precisely match the local surface plasmon resonance (LSPR) peak of the gold nanoparticles, in order to accurately simulate the excitation dynamics of the collective oscillation of free electrons on the surface of the gold nanoparticles.

[0076] Periodic Bloch conditions are set on the four sides of the three-dimensional finite element computational domain model to simulate an infinitely large area nanoparticle array structure; electric wall boundary conditions are set on the bottom surface of the substrate; the outermost electromagnetic wave absorption layer (i.e., the perfectly matched layer) retains the electromagnetic wave absorption properties to absorb the emitted electromagnetic waves and eliminate the interference of boundary reflections on the optical field calculation.

[0077] S202. Set the time-domain solver parameters and run the solution. After the solution is completed, perform data extraction and optical field parameter calculation.

[0078] The total simulation duration for transient solutions was set to 200 fs in the solver, and a fixed time step of 0.5 fs was used to completely capture the near-field optical field evolution process of the femtosecond pulsed laser on the surface of the target geometry.

[0079] After the solution is obtained, the local electric field intensity at the tip of the target geometry and the gap space (spacing ≤ 10 nm) between adjacent particles is extracted. Due to the excitation of local surface plasmon resonance, local optical field "hot spots" are formed at the tip and the gap space, and their local electric field enhancement factor F satisfies the following formula:

[0080]

[0081] In the formula, The local electric field intensity extracted at the hotspot. The input is the incident laser electric field intensity. In this embodiment, the calculated local electric field enhancement factor F reaches 10~1000, thereby amplifying the incident light intensity by 10~1000 times, achieving efficient absorption of the incident laser energy by the free electrons of the gold nanoparticles. The LSPR light field enhancement simulation results of the gold nanoparticles on the gold thin film substrate are as follows: Figure 4 As shown. Simultaneously, the local optical energy density distribution data (i.e., the optical energy density distribution in the nanoparticle region) of the target geometric surface is extracted and output.

[0082] S3. The local light field energy density distribution data is used as a volume heat source term and solved based on the dual-temperature model to output the spatial temperature distribution data and atomic vapor parameters of the target geometric structure.

[0083] In one embodiment of the present invention, step S3 specifically includes the following sub-steps:

[0084] S301. Construct partial differential equations and set dual-temperature model parameters to achieve coupling between light field and energy deposition.

[0085] A partial differential equation solving module (using a coefficient-type partial differential equation interface) is added to the simulation architecture, and the electronic temperature corresponding to the two-temperature model is defined. With lattice temperature Two solution variables. Then, the baseline parameters of the two-temperature model of the target geometry (i.e., gold material) are entered: the electron-phonon coupling coefficient g is set to 2 × 10⁻⁶. 16 W / (m 3 K), electronic heat capacity With lattice heat capacity All were characterized using temperature-dependent empirical formulas, and the initial values ​​of the electron temperature and lattice temperature were set to room temperature of 300K.

[0086] After the parameters are set, the local light field energy density distribution data stored in step S2 is retrieved as the volume heat source term S(t) and loaded into the region where the target geometry is located, thereby realizing the unidirectional coupling of light field energy and material energy deposition process.

[0087] S302. Set the segmented solution time step and perform transient solution calculations based on the dual-temperature model.

[0088] Before performing the calculation, the total simulation time was set to 10 ps, ​​and a segmented time step scheme was adopted for the solution. The time step in the interval from 0 to 100 fs was set to 0.5 fs, and the time step in the interval from 100 fs to 10 ps was set to 10 fs, so as to ensure the calculation accuracy while taking into account the simulation running efficiency.

[0089] During the transient solution calculation, the temperature evolution of electrons and the lattice within the target geometry satisfies the following two-temperature model equation:

[0090]

[0091]

[0092] In the formula, , Here, denoted as electronic heat capacity and lattice heat capacity, respectively; g is the electron-phonon coupling coefficient; and S(t) is the input laser thermal source term. During the femtosecond-scale energy deposition process, the electron temperature rapidly increases within 100 fs, and then transfers energy to the lattice within 1–10 ps through electron-phonon coupling, causing a sharp increase in the lattice temperature of the gold nanoparticles.

[0093] S303. Extract temperature evolution data to determine the vaporization state, and output temperature distribution and atomic vapor parameters.

[0094] After the solution is obtained, the lattice temperature evolution curve of the target geometry is extracted. When it is confirmed that the lattice temperature rises to above the boiling point of gold (2856℃, i.e., 3129K) within the range of 1~10 ps, ​​it is determined that the target geometry has achieved complete and instantaneous vaporization, forming high-density gold atom vapor.

[0095] After the vaporization is determined to be complete, the spatial temperature distribution of the vaporization region is extracted and output as the spatial temperature distribution data, and the total mass parameter of the vaporization region is output as the atomic vapor parameter.

[0096] S4. Input the spatial temperature distribution data and the atomic vapor parameters into the preset ionization rate equation and solve it to obtain the plasma number density distribution data.

[0097] In one embodiment of the present invention, step S4 specifically includes the following sub-steps:

[0098] S401. Add a solution interface and set the composite model parameters for the ionization rate equation.

[0099] After completing the energy deposition and gasification simulations, a partial differential equation interface was added to the computing system to customize the preset ionization rate equation. Specifically, the ionization rate equation adopts a composite model of multiphoton ionization and tunneling ionization. In the partial differential equation interface, the multiphoton ionization order and tunneling ionization field strength threshold corresponding to the 800nm ​​femtosecond laser are input, and these two are configured as the basic physical parameters for calculating plasma generation.

[0100] S402. Import initial conditions and execute the solution to simulate the strong field ionization and energy conversion process.

[0101] The spatial temperature distribution data and atomic vapor parameters (i.e., the gold atom vapor temperature and density obtained in the previous step) from step S3 are retrieved and used as initial conditions to solve the preset ionization rate equation. During the solution evolution process, the vaporized gold atom vapor undergoes further multiphoton ionization and tunneling ionization (i.e., strong field ionization) under the influence of the residual strong laser field. This strong field ionization process rapidly converts laser light energy into the internal energy of gold nanoparticles and the thermal energy of plasma, thereby forming a high-temperature, high-density plasma.

[0102] S403. Extract ion parameters and perform numerical verification, outputting plasma number density distribution data.

[0103] After solving the ionization rate equation, the electron number density and ion number density of the plasma are extracted. Subsequently, the extracted values ​​are validated to confirm whether the calculated electron and ion number densities fall within the range of 10. 20 ~1022 cm -3 Within the specified range. After successful verification, the extracted results, including the electron number density and ion number density, are output as the plasma number density distribution data, and a storage operation is performed.

[0104] S5. Input the plasma number density distribution data into the preset fluid control equation for solution to obtain the plasma expansion velocity and mass ejection rate, and convert the plasma expansion velocity and mass ejection rate into surface load data according to the principle of momentum conservation.

[0105] In one embodiment of the present invention, step S5 specifically includes the following sub-steps:

[0106] S501, Add the fluid control equation interface and set the initial and boundary conditions for the expansion simulation.

[0107] An Euler equation interface for compressible fluid flow is added to the simulation system as the preset fluid control equation to simulate the supersonic expansion of plasma into vacuum. The plasma number density distribution data output in step S4 is retrieved, and the following initial and boundary conditions are set:

[0108] The initial temperature of the plasma is set to 10. 4 ~10 6 K, while the initial pressure of the plasma is calculated using the ideal gas law. Regarding boundary condition configuration, the upper boundary of the computational domain is set as the vacuum exit boundary, while periodic boundary conditions are applied to the sides of the computational domain.

[0109] S502. Set the solver time step to perform transient solution, extract and verify the expansion velocity and injection rate.

[0110] After setting up, a transient solution was performed, with a total duration of 100 ps and a time step of 1 ps. During the transient solution process, high-temperature, high-density plasma contains high-temperature plasma with a density as high as 10... 12 A huge pressure gradient of Pa / m drives it to expand supersonically toward the vacuum.

[0111] After solving the problem, the plasma velocity field and mass flux distribution are first obtained. Then, in post-processing, the normal plasma ejection velocity and mass ejection rate on the substrate surface (i.e., the upper surface of the substrate) are integrally analyzed using an area integral to extract the plasma expansion velocity and mass ejection rate. After extraction, the values ​​are verified, confirming that the plasma ejection velocity falls within the range of 10–100 km / s.

[0112] S503. Calculate the recoil force based on the principle of momentum conservation and convert it into surface load data.

[0113] After obtaining the plasma expansion velocity and mass ejection rate, the total recoil force acting on the substrate surface (i.e., the substrate surface) is calculated according to the principle of momentum conservation (i.e., the momentum conservation formula). This recoil force F satisfies the following formula:

[0114]

[0115] In the formula, m is the mass of the ejected plasma, v is the expansion velocity of the plasma, and dm / dt is the plasma mass ejection rate. This calculation process converts the thermal energy of the plasma into mechanical energy, and the calculated total recoil force is output and stored as the surface load data.

[0116] S6. Apply the surface load data to the substrate surface in the three-dimensional finite element calculation domain model, and perform transient solution based on the solid mechanics equations to obtain the predicted dynamic pressure value in the substrate.

[0117] In one embodiment of the present invention, step S6 specifically includes the following sub-steps:

[0118] S601. Set the basic parameters of the solid mechanics interface and apply surface load data and boundary conditions.

[0119] After obtaining the recoil force parameters (i.e., the surface load data) from step S5, a solid mechanics interface is added to the simulation system, and its effective range is limited to the substrate (i.e., the gold substrate region). The basic mechanical parameters of the substrate are input, specifically including: density 19300 kg / m³. 3 Elastic modulus 79 GPa, Poisson's ratio 0.42.

[0120] Regarding the setting of loads and boundary conditions, the recoil force calculated in the previous step is used as the surface load data and applied to the substrate surface (i.e., the upper surface of the substrate) in the three-dimensional finite element calculation domain model in the form of a surface load. The time history of this load is set to correspond completely with the time curve of plasma jet. At the same time, a low-reflection boundary is set on the bottom surface of the substrate to eliminate the interference of stress wave bottom surface reflection. Periodic boundary conditions are retained on the side of the calculation domain.

[0121] S602. Set the time-domain solver parameters and solve the stress wave propagation process based on the solid mechanics equations.

[0122] After configuring the boundary and load parameters, a transient solution based on solid mechanics equations is executed. The total duration of the transient solution is set to 500 ps, ​​and the time step is 5 ps. During the solution process, the propagation and superposition of longitudinal compressive stress waves inside the substrate (gold substrate) are simulated and solved.

[0123] S603. Extract the compressive stress distribution and combine it with the effective area of ​​laser irradiation to calculate and output the predicted dynamic pressure value.

[0124] After the transient solution is completed, the compressive stress distribution inside the matrix is ​​extracted through post-processing. Subsequently, the dynamic pressure peak is calculated based on the effective area A of the laser irradiation. The dynamic pressure P generated inside the target material satisfies the following formula:

[0125]

[0126] In the formula, F is the recoil force, and A is the effective area of ​​the laser irradiation. The load calculation term is derived from the conservation of momentum.

[0127] During the calculation process, the pressure distribution after the superposition of stress waves from multiple adjacent hot spots was simultaneously statistically analyzed to verify and confirm that the peak value reached the GPa level. Finally, the dynamic pressure peak value and distribution results calculated by combining the compressive stress distribution and the effective area were output as the predicted dynamic pressure value in the matrix, completing the full-process simulation from light field enhancement to ultra-high pressure generation. Furthermore, as... Figure 5 The curves showing the variation of plasma recoil pressure of Au and Ag nanoparticles with particle spacing further verify the regulatory effect of different particle spacing sizes on the final high pressure intensity.

[0128] It should be noted that those skilled in the art should understand that the above-mentioned process of solving multiphysics coupling and partial differential equations can be implemented not only by writing computer program code (e.g., calling finite element algorithm libraries based on languages ​​such as C++ and Python), but also by directly using the corresponding physics interfaces in existing commercial multiphysics numerical simulation platforms (e.g., but not limited to COMSOL Multiphysics, ANSYS, and other finite element analysis software). This invention does not impose any unique limitation on the specific underlying implementation tools.

[0129] It should be noted that, in this document, the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, apparatus, article, or method that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such process, apparatus, article, or method. Unless otherwise specified, an element defined by the phrase "comprising one..." does not exclude the presence of other identical elements in the process, apparatus, article, or method that includes that element.

[0130] The above description is merely a preferred embodiment of the present invention and does not limit the patent scope of the present invention. Any equivalent structural or procedural transformations made based on the content of the present invention's specification and drawings, or direct or indirect applications in other related technical fields, are similarly included within the patent protection scope of the present invention.

Claims

1. A method for predicting plasma recoil pressure using multiphysics coupling, characterized in that, include: S1. Establish a three-dimensional finite element computational domain model that includes the substrate and the target geometric structure located on the surface of the substrate; S2. Obtain the initial electromagnetic wave parameters of the pulsed laser, and perform electromagnetic wave time-domain solution on the three-dimensional finite element calculation domain model based on the initial electromagnetic wave parameters to extract the local light field energy density distribution data of the target geometric surface. S3. The local light field energy density distribution data is used as a volume heat source term and solved based on the dual-temperature model to output the spatial temperature distribution data and atomic vapor parameters of the target geometric structure. S4. Input the space temperature distribution data and the atomic vapor parameters into the preset ionization rate equation for solution to obtain plasma number density distribution data; S5. Input the plasma number density distribution data into the preset fluid control equation for solution to obtain the plasma expansion velocity and mass ejection rate, and convert the plasma expansion velocity and mass ejection rate into surface load data according to the principle of momentum conservation. S6. Apply the surface load data to the substrate surface in the three-dimensional finite element calculation domain model, and perform transient solution based on the solid mechanics equations to obtain the predicted dynamic pressure value in the substrate.

2. The method for predicting plasma recoil pressure using multiphysics coupling according to claim 1, characterized in that, S1 specifically includes: S101. Establish a cuboid solid domain as the base, and construct a unit cell model including spherical solids as the target geometry on the surface of the base based on periodic boundary conditions. Set at least two spherical solids in the unit cell model to form gap spaces between adjacent spherical solids. Establish a fluid expansion domain above the target geometry, and wrap an electromagnetic wave absorbing layer around the base, the target geometry, and the fluid expansion domain to reconstruct a multi-layer geometric computational domain. S102. Spatial discretization is performed on the multi-layer geometric computational domain. Different sizes of meshes are used to perform partitioned meshing on the target geometric structure and the gap space, the matrix and the fluid expansion domain, and the electromagnetic wave absorption layer. After the partitioning is completed, the mesh distortion rate is checked to optimize the distorted mesh elements, thus completing the establishment of the three-dimensional finite element computational domain model.

3. The method for predicting plasma recoil pressure using multiphysics coupling according to claim 1, characterized in that, S2 specifically includes: S201. Assign material and electromagnetic properties to the three-dimensional finite element computational domain model, and set pulsed laser source and boundary conditions in the electromagnetic wave time domain solution architecture. S202. Set the time-domain solver parameters and run the solution. After the solution is completed, extract the local electric field intensity from the tip and gap space of the target geometry, calculate the local electric field enhancement factor, and output the local light field energy density distribution data.

4. The method for predicting plasma recoil pressure using multiphysics coupling according to claim 3, characterized in that, In step S202, the formula for calculating the local electric field enhancement factor is as follows: Where F is the local electric field enhancement factor, and its value ranges from 10 to 1000. For the local electric field strength, The input is the incident laser electric field intensity.

5. The method for predicting plasma recoil pressure using multiphysics coupling according to claim 1, characterized in that, S3 specifically includes: S301. Define electronic temperature and lattice temperature as solution variables, set the parameters of the dual-temperature model, and load the local light field energy density distribution data as the volume heat source term into the region where the target geometry is located; wherein, the parameters of the dual-temperature model include electron-phonon coupling coefficient, electronic heat capacity, and lattice heat capacity; S302. Set a segmented time step and perform transient solution for the evolution process of the electron temperature and the lattice temperature based on the dual-temperature model. S303. Extract the lattice temperature evolution curve of the target geometry, and when it is determined that the lattice temperature rises above the boiling point, extract and output the temperature spatial distribution result and total mass parameter of the vaporization region, and use the temperature spatial distribution result and the total mass parameter as the spatial temperature distribution data and the atomic vapor parameter, respectively.

6. The method for predicting plasma recoil pressure using multiphysics coupling according to claim 5, characterized in that, In step S302, the dual-temperature model satisfies the following equations during the solution process: in, For the electronic heat capacity, The lattice heat capacity, The electron temperature, The lattice temperature is given, g is the electron-phonon coupling coefficient, S(t) is the bulk heat source term, and t is time.

7. The method for predicting plasma recoil pressure using multiphysics coupling according to claim 1, characterized in that, S4 specifically includes: S401. Define a preset ionization rate equation in the partial differential equation interface and set the composite model parameters of the ionization rate equation; wherein, the ionization rate equation adopts a composite model of multiphoton ionization and tunneling ionization, and the composite model parameters include the multiphoton ionization order and the tunneling ionization field strength threshold. S402. Input the spatial temperature distribution data and the atomic vapor parameters as initial conditions into the preset ionization rate equation to perform transient solution of the strong field ionization process; S403. Extract the obtained electron number density and ion number density and perform validity verification. After the verification is passed, output the result containing the electron number density and the ion number density as the plasma number density distribution data.

8. The method for predicting plasma recoil pressure using multiphysics coupling according to claim 1, characterized in that, S5 specifically includes: S501. Set the initial and boundary conditions for simulating the supersonic expansion of plasma in the preset fluid control equation; wherein, the preset fluid control equation adopts the Euler equation for compressible fluid flow, and the initial conditions include the initial plasma temperature and the initial plasma pressure. S502. Set time-domain solution parameters to perform transient solution, extract the plasma velocity field and mass flux distribution, and perform surface integration operation on the substrate surface to extract the plasma expansion velocity and mass ejection rate. S503. Calculate the total recoil force acting on the substrate surface according to the momentum conservation formula, and convert the total recoil force into the surface load data output.

9. The method for predicting plasma recoil pressure using multiphysics coupling according to claim 1, characterized in that, S6 specifically includes: S601. Input the basic mechanical parameters of the matrix in the solid mechanics interface corresponding to the matrix, apply the surface load data to the surface of the matrix in the form of surface load, and set the load time history and solution boundary conditions. S602. Set the time-domain solver parameters and perform transient solution of the propagation and superposition process of longitudinal compressive stress wave inside the matrix based on the solid mechanics equations. S603. Extract the compressive stress distribution inside the matrix, calculate the dynamic pressure peak value by combining it with the effective area of ​​laser irradiation, and output the data containing the dynamic pressure peak value and its superimposed distribution result as the predicted value of dynamic pressure in the matrix.

10. The method for predicting plasma recoil pressure using multiphysics coupling according to claim 9, characterized in that, In S603, the formula for calculating the peak dynamic pressure is: Wherein, P is the peak dynamic pressure, F is the corresponding total recoil force, A is the effective area of ​​laser irradiation, and v is the plasma expansion velocity. This refers to the mass spray rate.

Images