An ultrasound enhancement device modeling and optimization method

Abstract

This invention discloses a modeling and optimization method for an ultrasonic enhancement device. It analyzes the device's acoustic field signal, designs a three-dimensional geometric configuration including a concave energy-concentrating structure, a metamaterial microstructure, and an acoustic lens structure, sets corresponding material and geometric parameters for each structure, evaluates the performance of the current configuration through nonlinear acoustic field simulation, and accurately simulates the simulation results of the heterogeneous interface. Then, it inputs the main field and adjoint field dual-channel solutions to obtain the sensitivity information of the design variables. Furthermore, it constructs a deep neural network model, using the geometric parameters of the metamaterial microstructure as input and the equivalent acoustic response as output, to achieve real-time performance prediction and dynamic impedance matching optimization. Combined with objective function evaluation and constraint projection mechanisms, a simulation optimization feedback loop is formed, iteratively outputting the optimal manufacturable configuration. The device designed in this invention has advantages such as being passive, low-cost, and highly adaptable, and is suitable for complex industrial applications such as partial discharge monitoring of power equipment and fault diagnosis in rail transit.

Description

Technical Field

[0001] This invention relates to the fields of computational acoustics and deep learning technology, and in particular to a method for modeling and optimizing an ultrasonic enhancement device. Background Technology

[0002] With the rapid development of Physics-Informed AI (PIAI) technology, deep learning and multiphysics simulation have made significant progress in fields such as engineering design, material inverse design, and intelligent sensing. Particularly in acoustic metamaterial design, graph neural networks, surrogate models, and topology optimization methods have been used to automatically generate microstructural units with specific acoustic responses. However, existing ultrasonic enhancement device design methods still primarily rely on simplified linear acoustic models, ideal interface assumptions, and isolated parameter scanning strategies, lacking the ability to model the complexities of real physical environments. These traditional methods generally suffer from the following key problems:

[0003] First, the physical model is oversimplified. Most simulations ignore nonlinear effects such as harmonic generation, waveform distortion, and thermoviscous absorption in high-intensity sound fields, resulting in simulation results that deviate significantly from measured performance.

[0004] Secondly, the interface behavior modeling is coarse. For the transmission, reflection and phase shift at concave surfaces, curved surfaces or heterogeneous material interfaces, ideal impedance matching or one-dimensional transmission models are often used, which cannot accurately reflect the real scattering characteristics.

[0005] Third, the lack of a multi-scale coupling mechanism leads to the separation of subwavelength modulation of metamaterial microstructures from macroscopic sound field propagation. An end-to-end differentiable mapping from unit geometry to system performance has not been established, making it difficult to achieve global optimization.

[0006] For example, existing designs based on linear wave equations cannot predict focal waveform distortion at high sound pressure levels, leading to distorted sensor signals. Devices using fixed-geometry metamaterial arrays, however, often suffer from sidelobe enhancement, gain reduction, or even structural failure in practical deployments due to a lack of consideration for manufacturing constraints and material nonlinear responses. More critically, current intelligent optimization methods often directly replace physical simulations with deep learning models, lacking mechanisms to correct prediction errors, thus disrupting gradient consistency and rendering optimization results unreliable or unmanufacturable.

[0007] How to construct a closed-loop design framework that integrates high-fidelity nonlinear modeling, precise interface processing, multi-scale collaborative optimization, and manufacturability constraints while ensuring acoustic physical authenticity, so as to achieve efficient, stable, and passive enhancement of weak ultrasonic signals in the air, has become a core challenge in the field of industrial intelligent sensing and nondestructive testing. Summary of the Invention

[0008] To address the aforementioned technical problems in the existing technology, this invention proposes an ultrasonic enhancement device and its optimization method, the specific technical solution of which is as follows:

[0009] A method for modeling and optimizing an ultrasonic enhancement device includes the following steps:

[0010] Step 1: Analyze and process the acoustic field signal of the operating equipment, design the initial geometric configuration of the ultrasonic enhancement device, and generate a parameter data package of the initial geometric configuration. The structure of the initial geometric configuration includes a concave energy-concentrating structure, a metamaterial microstructure, and a lens structure.

[0011] Step 2: Based on the parameter data package, create a three-dimensional acoustic calculation domain and set the corresponding boundary conditions and material properties for each structure;

[0012] Step 3: Perform nonlinear sound field simulation to evaluate the performance of the current geometric configuration, in order to determine whether the current geometric configuration should be optimized;

[0013] Step 4: For different interface types of the current geometric configuration, perform heterogeneous interface hybrid modeling and solve to obtain interface simulation results;

[0014] Step 5: Based on the boundary conditions, material properties, and interface simulation results, perform full-wave simulation on the current geometric configuration, and conduct sensitivity analysis on the full-wave simulation results according to the preset target results; use the parameter inversion method to obtain the true acoustic response of the metamaterial microstructure from the full-wave simulation results;

[0015] Step 6: Introduce a deep neural network surrogate model, take the geometric parameters of the metamaterial microstructure as input, quickly predict the equivalent acoustic response, and use the real acoustic response to compare the error, thereby triggering a local correction mechanism to optimize and update the surrogate model;

[0016] Step 7: Combining the sensitivity analysis results from Step 5 and the predicted equivalent acoustic response results from Step 6, and using an objective function evaluation and constrained projection mechanism, dynamically update the design parameters of the current geometric configuration and output the optimal device configuration that can be used for additive manufacturing.

[0017] Furthermore, in step 1, the parameters of the device's acoustic field signal include the ultrasonic operating frequency, the target focal position, and the medium properties;

[0018] The concave energy-concentrating structure is used to achieve primary energy convergence of ultrasound waves. This structure adopts a parabolic or spherical geometric configuration.

[0019] The acoustic metamaterial microstructure is used for mesoscale phase compensation and impedance matching. The structure adopts a cross-shaped, H-shaped, labyrinth-shaped or spiral topology.

[0020] The acoustic lens structure achieves wavefront shaping and focusing through refractive index distribution. This structure adopts the form of a gradient refractive index lens or a planar superlens.

[0021] Furthermore, in step 1, the manufacturability constraints of the initial geometric configuration of the design are checked, specifically: the minimum feature size is checked; the volume fraction is checked, and the proportion of solid material in the total volume of the metamaterial microstructure is calculated; the topological connectivity is checked, and the physical connections of the metamaterial microstructure are verified by image morphology methods.

[0022] Furthermore, step 2 specifically includes:

[0023] Step 2.1: Extract the geometric parameter information of the concave energy-concentrating structure, metamaterial microstructure and lens structure from the parameter data package, convert them into boundary representations and embed them into a cubic computational domain. Then perform Boolean operations to separate each solid structure from the background medium, forming an air medium region and a solid region, and finally generate a multi-connected computational domain.

[0024] Step 2.2: Call the 3D mesh generator based on Delaunay triangulation to perform unstructured tetrahedral meshing on the multi-connected computational domain; among which, adaptive local refinement is performed in three key regions: first, in the neighborhood of the target focus; second, in the gaps between metamaterial elements; and third, near the solid-air interface.

[0025] Step 2.3: Traverse all meshes and automatically assign corresponding material parameters based on the geometric configuration region to which the centroid of a mesh belongs;

[0026] Step 2.4: Apply three types of physical boundary conditions to the boundary of the computational domain: First, specify a circular region on the back side of the concave energy-concentrating structure as the ultrasonic incident surface and apply plane wave or focusing source excitation; Second, set a perfectly matched layer absorbing boundary on the outer surface of the cubic computational domain to simulate an infinite open space and suppress false reflections; Third, apply symmetric boundary conditions on the symmetric plane if the geometry is symmetric.

[0027] After the boundary settings are completed, perform the final mesh quality verification.

[0028] Furthermore, the multi-connected computational domain includes a concave energy focusing domain, a metamaterial equivalent domain, and a lens-air interface domain. In the metamaterial equivalent domain, a multi-scale homogenization technique is introduced to transform the metamaterial microstructure and equivalent medium method into a continuous medium with equivalent refractive index and impedance.

[0029] Furthermore, step 3 specifically includes:

[0030] Step 3.1: Set the corresponding excitation source according to the input ultrasonic working frequency, and use the direct sparse matrix solution method to describe the propagation behavior of high-intensity ultrasonic waves in the multi-connected computational domain using the Westervelt nonlinear acoustic equation;

[0031] Step 3.2: In the neighborhood of the target focus, firstly, the actual peak sound intensity point is located by searching for local maxima; then, sound intensity profiles are intercepted along the three orthogonal directions (x, y, z), and the full width at half maximum (FWHM) is calculated as the lateral and axial resolution indicators; next, the first local maximum outside the main lobe is identified, and its ratio to the main lobe peak value is calculated to obtain the sidelobe suppression ratio; finally, the volume in the focus region where the sound intensity exceeds 50% of the peak value is counted.

[0032] Step 3.3: Assume that the heat source power is proportional to the sound intensity, calculate the steady-state heat deposition distribution based on the sound field, and estimate the local temperature rise;

[0033] Step 3.4: Based on the focus positioning accuracy, lateral focus resolution, sidelobe suppression ratio obtained in Step 3.2 and the maximum temperature rise obtained in Step 3.3, compare them with the preset thresholds one by one. If all comparison results meet the threshold requirements, the geometric configuration design is successful; otherwise, proceed to the optimization process.

[0034] Furthermore, step 4 specifically includes:

[0035] Step 4.1: Automatically scan all material interface regions. For the curved interface between the concave energy-concentrating structure and the air, use the boundary element method (BEM). For the metamaterial microstructure, use the finite element method (FEM), with the nonlinear weak form expression of the sound field as the governing equation. Solve the spatiotemporal distribution of sound pressure under a multi-scale grid to obtain the main field results including nonlinear propagation effects. For the interface between the metamaterial microstructure and the acoustic lens that is approximately planar or layered, use the transfer matrix method (TMM).

[0036] Step 4.2: For the highly curved surface interface between the concave energy-concentrating structure and the air, the Boundary Element Method (BEM) is used to solve the frequency domain Helmholtz equation to calculate the scattered sound field and extract the sound pressure boundary data.

[0037] Step 4.3: For the interface between the metamaterial microstructure and the acoustic lens that is approximately planar or layered, the Transmission Matrix Method (TMM) is used. Each layer is regarded as a homogeneous equivalent medium. The characteristic matrix chain product is extracted and constructed based on the equivalent acoustic impedance and thickness of each layer. The transmission and reflection coefficients and the cumulative phase delay are calculated. The transmitted wave is subjected to layered transmission, phase compensation and wavefront shaping so that the energy distribution in the focusing area meets the design target.

[0038] Step 4.4: For interfaces with strong curvature surfaces and interfaces with approximate planes or layered stacks, implement unified coupling between the Boundary Element Method (BEM) region and the Transmission Matrix Method (TMM) region in the frequency domain, and establish unified boundary conditions: Force the sound pressure continuity and normal particle velocity matching constraints at the boundary between the regions of the two methods, and embed the above interface conditions in the global linear system by introducing Lagrange multipliers or penalty functions.

[0039] Finally, differential calculations are performed on the unified boundary conditions to extract the energy comparison information between the focal and sidelobe regions, and the error and sound field distribution results are output.

[0040] Furthermore, step 5 specifically includes:

[0041] Step 5.1: Construct the master field solver, input the boundary conditions, material properties and interface simulation results, and output the distribution of sound pressure in the spatiotemporal domain;

[0042] Step 5.2: Simultaneously construct the adjoint field solver, which solves for sensitivity information under the backpropagation condition of the given objective function, and obtains the design variables of the current geometric configuration;

[0043] Specifically, using the adjoint method theory, the gradient of any design variable with respect to the objective function can be expressed as the inner product of the main field and the adjoint field solvers in the perturbation region, without the need for repeated forward simulations. It automatically identifies the geometric or material perturbation operators corresponding to each variable and calculates the sensitivity.

[0044] Step 5.3: Package the gradient vectors of all design variables with the current performance metrics for subsequent optimization feedback of the geometric configuration;

[0045] Step 5.4: For metamaterial microstructures, based on their local incident and transmitted acoustic field data, the parametric inversion method is used to calculate their true acoustic response at the operating frequency.

[0046] Furthermore, step 6 specifically includes:

[0047] Step 6.1: Introduce a deep neural network model. The input is the normalized geometric parameters of the metamaterial microstructure, and the output is the predicted equivalent acoustic response. The deep neural network surrogate model mainly consists of an input layer, a DNN surrogate model body, an output layer, and a monitoring and correction layer. The DNN surrogate model body includes a feature extraction block, an activation layer, a pooling layer, and a fully connected layer connected in sequence.

[0048] Step 6.2: The deep neural network model is pre-trained using 5000 sets of offline-generated full-wave simulation datasets before optimization begins. At the beginning of each optimization iteration, the model is called to quickly predict the macroscopic acoustic response of candidate metamaterial microstructures.

[0049] Step 6.3: Each time the surrogate model compares the real acoustic response with the equivalent acoustic response, if the relative error of any output component exceeds 5%, a local correction mechanism is triggered by monitoring the correction layer: the input sample is added to the training set, fine-tuning is performed only on the weights of the last two layers of the DNN surrogate model, and the DNN is restricted from participating in the adjoint gradient calculation, serving only as an acceleration tool to provide the initial configuration and rapid evaluation; finally, the corrected DNN surrogate model is reconnected.

[0050] Furthermore, step 7 specifically includes:

[0051] Step 7.1: Define the objective function, using the weighted difference between maximizing the energy in the focal region and minimizing the energy in the sidelobe region as the optimization guide, where the focal energy weight coefficient and the sidelobe suppression weight coefficient can be dynamically adjusted;

[0052] Step 7.2: The optimizer integrates the gradient information provided in Step 5 with the equivalent acoustic response predicted by the model in Step 6 to update the design variables.

[0053] Step 7.3: Immediately after each parameter update, perform three hard constraint checks on the new parameters: 1. Minimum feature size; 2. Overall integral; 3. Structural connectivity; If any constraint is violated, pull the parameters back into the feasible region through a projection operation;

[0054] Step 7.4: Continuously monitor the improvement rate of the objective function. In each round, check the relative rate of change of the objective function to determine convergence. Output the final three-dimensional geometric model, material distribution map and process parameter package for use in additive manufacturing equipment driving.

[0055] The advantages and beneficial effects of this invention are as follows:

[0056] 1. Deeply Coupled Physics-Algorithm Closed-Loop Design: This invention achieves deep integration of physical hardware structure and intelligent optimization algorithm by constructing a closed-loop iterative mechanism of "high-fidelity simulation—accompanied sensitivity analysis—manufacturability constraint optimization—parameter feedback update". Compared with the traditional unidirectional design process, this closed loop can dynamically correct model bias, suppress overfitting, and ensure that the final configuration has both high performance and engineering feasibility.

[0057] 2. High-fidelity modeling of multi-scale nonlinear sound fields: This invention innovatively integrates multi-scale homogenization technology with pressure acoustic transient equations containing a nonlinear parameter β, enabling accurate characterization of harmonic generation, waveform distortion, and thermoviscous dissipation effects of high-intensity ultrasound in composite metamaterials. Compared to linear models that ignore nonlinearity, this significantly improves the prediction accuracy of real focusing behavior, making it particularly suitable for strong excitation industrial scenarios ranging from 20 to 500 kHz.

[0058] 3. Hybrid Numerical Solution Strategy for Heterogeneous Interfaces: For complex interfaces between concave energy-concentrating structures, metamaterial arrays, and acoustic lenses, this invention employs an adaptive hybrid solution framework combining the boundary element method and the transfer matrix method, automatically matching the optimal numerical method based on the interface's geometric characteristics. This strategy significantly reduces computational overhead while ensuring high accuracy in curved surface scattering and multilayer transmission, achieving efficient and seamless integration of the entire sound field.

[0059] 4. Synergistic Mechanism for Accelerating and Ensuring Robustness of the Surrogate Model: This invention introduces a deep neural network surrogate model to rapidly predict the equivalent acoustic response of metamaterials, significantly shortening the exploration cycle in the initial optimization phase. Simultaneously, through error monitoring and local fine-tuning mechanisms, the surrogate model maintains physical consistency throughout the entire iteration process. This mechanism improves overall optimization efficiency by more than three times without sacrificing accuracy and effectively avoids convergence failures caused by model inaccuracies.

[0060] 5. Embedded Manufacturability Constraints Throughout the Entire Process: This invention embeds engineering constraints such as minimum feature size, total integral, and structural connectivity in real time during optimization iterations. Projection correction ensures that each generation of design schemes can be directly used for additive manufacturing. This bridges the last mile between "digital design and physical realization," solving the industry pain point of traditional acoustic metamaterial design being "excellent in simulation but unable to be manufactured." Attached Figure Description

[0061] Figure 1 This is a flowchart illustrating the principle of a modeling and optimization method for an ultrasonic enhancement device according to an embodiment of the present invention.

[0062] Figure 2 This is a schematic diagram of multi-scale nonlinear sound field modeling and computational domain partitioning according to an embodiment of the present invention;

[0063] Figure 3 This is a structural diagram of the heterogeneous interface hybrid numerical solution module according to an embodiment of the present invention;

[0064] Figure 4 This is an input-output structure diagram of the deep neural network model in an embodiment of the present invention. Detailed Implementation

[0065] To make the objectives, technical solutions, and technical effects of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and embodiments.

[0066] like Figure 1As shown, this embodiment discloses a modeling and optimization method for ultrasonic enhancement devices based on physical-algorithm coupled feedback and multi-scale nonlinear modeling. It constructs a closed-loop collaborative framework integrating high-fidelity sound field simulation, accompanying sensitivity analysis, neural proxy acceleration, and manufacturability constraint projection. Furthermore, it integrates multi-scale homogeneous modeling, nonlinear acoustic equations, and hybrid numerical solution techniques for heterogeneous interfaces to perform multi-level passive focusing and phase correction on the raw ultrasonic waves radiated by the operating equipment. This solves the problems of physical distortion, rough interface processing, unreliable optimization, and unmanufacturable configurations in existing design methods. Specifically, the method includes the following steps:

[0067] Step 1: Analyze and standardize the multi-source physical input parameters in the scene to generate and initialize the physical hardware structure link that meets the manufacturing constraints. The multi-source physical input parameters include the ultrasonic working frequency, the target focal position, and the medium properties. The physical hardware structure link is composed of a concave energy-concentrating structure, an acoustic metamaterial microstructure, and an acoustic lens structure cascaded in sequence.

[0068] Specifically, the system design first includes an input analysis module and a geometry initialization module. The input analysis module analyzes and standardizes key parameters from multi-source physical data. These key parameters include the ultrasonic operating frequency, the target focal point location, and the fundamental physical properties of the medium, such as sound velocity, density, nonlinear coefficient, and thermoviscosity parameters. The input analysis module performs unit consistency and data integrity checks and establishes a standardized parameter set to provide a unified input for subsequent calculations.

[0069] Subsequently, the geometry initialization module is invoked to generate the initial geometric configurations of the concave energy-concentrating structure, the acoustic metamaterial microstructure, and the acoustic lens structure based on the target focal position and medium properties. The curvature radius of the concave energy-concentrating structure, the period and geometric configuration of the acoustic metamaterial microstructure, and the refractive index distribution of the acoustic lens structure are initialized.

[0070] The concave energy-concentrating structure is used to achieve primary energy focusing of ultrasound waves. This structure adopts a parabolic or spherical geometric configuration, and its radius of curvature R satisfies the relationship between the working wavelength λ and the curvature R. The boundary opening diameter D satisfies After optimization, the primary sound pressure gain can exceed 15 dB, and the focal size is less than half the wavelength λ / 2.

[0071] Acoustic metamaterial microstructures are used for mesoscale phase compensation and impedance matching. These structures employ cross-shaped, H-shaped, labyrinthine, or spiral topologies. In a single embodiment, they are fabricated using homogeneous materials, with the material type selected based on the application scenario: photosensitive resin is used for lightweight applications via stereolithography and 3D printing; aluminum alloy is used for high acoustic impedance matching via micromilling or laser sintering; and piezoelectric composite materials are used for active control applications via multi-material inkjet printing. Different materials are not mixed within the same unit array to ensure both fabrication feasibility and acoustic performance consistency. Unit periodicity. This is to satisfy the subwavelength modulation conditions and suppress higher-order diffraction.

[0072] The acoustic lens structure achieves wavefront reshaping and high-precision focusing through refractive index distribution. This structure employs either a gradient refractive index lens or a planar superlens. The gradient refractive index lens achieves a gradual refractive index change through spatially varying porosity or material filling rate. The planar superlens is composed of a phase-abrupt metasurface, with each metaatomic unit providing a phase coverage of 0 to 2π for achieving far-field or near-field subwavelength focusing. Through near-field evanescent wave coupling and phase gradient modulation, subwavelength focusing is achieved in the near-field region less than λ / 4 from the lens surface, with a focal size reaching 0.42λ, breaking the traditional far-field diffraction limit.

[0073] After the structure is generated, manufacturing constraint verification is performed to ensure that the geometric parameters meet the requirements of minimum feature size, material ratio and topological connectivity.

[0074] More specifically, the execution details of the system's input parsing module and geometry initialization module include the following steps:

[0075] Step 1.1: The system first receives multi-source physical input data from the user interface or external sensors, which constitute the initial boundary conditions for the entire acoustic focusing device design.

[0076] Specifically, the inputs include the ultrasound operating frequency f and the target focal position. and the set of fundamental physical parameters describing the properties of the propagation medium. ,in Represents the speed of sound of a small signal , static density , For dimensionless nonlinear parameters, The thermal viscosity absorption coefficient Because these parameters may originate from different devices or historical databases, their unit systems, precision levels, and numerical ranges can vary significantly. Therefore, the system first enters the parameter parsing phase. In this phase, the input parsing module performs unit consistency checks on each input item, such as automatically converting kHz to Hz. Convert to And for values ​​that clearly exceed the physically reasonable range, such as negative frequencies or sound speeds greater than [a certain value], [the following applies]. Marking is performed. For missing fields, the system calls the built-in working condition knowledge base to interpolate and complete the data based on historical data from similar scenarios; if there is no reliable reference, a preset default value is used, such as air quality. , After cleaning, all parameters are normalized to a set of predefined reference scales to eliminate the interference of dimensional differences on subsequent numerical calculations, ultimately forming a standardized input vector. The mathematical expression of this process is as follows:

[0077] ,

[0078] In this formula, This represents the normalized input parameter vector; f is the user-specified ultrasonic operating frequency; fref is the preset reference frequency, with a value of 100kHz, used for frequency normalization; xf is the position coordinate vector of the target focus in three-dimensional space; Lref is the reference length scale, with a value of 0.5m, used for spatial coordinate normalization; c0 is the small-signal sound velocity of the medium; c ref The reference speed of sound is taken as 343 m / s, corresponding to the speed of sound in standard air; ρ0 is the static density of the medium; ρ ref For reference density, a value of 1.2 kg / m³ is used. 3 , corresponding to standard air density; β is the nonlinear parameter of the medium; β ref For reference, the nonlinear parameter is set to 1.2, representing the nonlinearity level of a typical gas or liquid; δ is the thermoviscous absorption coefficient; δ ref For reference thermoviscosity coefficient, a value of 1.0 × 10⁻⁶ is used. -3 m 2 / S is used to normalize the absorption characteristics; the superscript T indicates the vector transpose operation, ensuring that the output is in column vector form.

[0079] Step 1.2: After obtaining the standardized parameter set, the system activates the geometry initialization module, which automatically generates the initial 3D model of the three-level passive acoustic structure based on physical principles. The core task of this module is to transform abstract acoustic requirements into concrete geometric entities, while also taking into account the functional positioning of each substructure.

[0080] First, to address the primary energy focusing requirements, the system constructs a concave energy-concentrating structure in the form of a parabolic or spherical cap, with its focal point precisely aligned with the user-specified target location X. f The radius of curvature R is dynamically set according to the working wavelength λ, satisfying the relationship This ensures effective convergence in the near-field region.

[0081] Secondly, to compensate for phase distortion and impedance mismatch during mesoscale propagation, the system deploys an acoustic metamaterial microstructure behind the concave surface. This microstructure consists of an array of periodically arranged cross-shaped units, hence it can also be called an acoustic metamaterial array. The geometry of each unit is defined by its characteristic dimension d, thickness h, and period. and rotation angle The initial values ​​are determined according to the subwavelength design criteria.

[0082] Finally, to achieve high-precision focusing, the system deploys a planar acoustic lens at the very front. Its surface is not made of a uniform material, but rather actively shapes the incident wavefront through a spatially varying equivalent refractive index or phase distribution. This phase distribution... The calculations are performed using the ideal spherical wave principle, ensuring that sound waves from all paths are superimposed in phase at the focal point. The initial settings of the above geometric parameters follow these rules:

[0083] In this formula, This represents the radius of curvature of the concave energy-concentrating structure; This is an empirical scaling factor, with a value range of [interval]. , used to adjust the primary focusing intensity; c0 is the small signal velocity of the medium; f is the ultrasonic operating frequency; d represents the period of the acoustic metamaterial unit; d represents the characteristic dimension of the metamaterial unit, such as the width of the cross arm; h represents the thickness of the metamaterial unit. Let be the rotation angle of the metamaterial unit, initially set to 0 radians; Indicates the spatial position of the acoustic lens The required phase delay, Let be the position vector of any point on the lens plane; The target focus location; This serves as a reference point on the lens, typically taken as the center of the lens; This represents the Euclidean norm, which is the straight-line distance between two points. The operating angular wavenumber is used to convert path difference into phase difference.

[0084] Step 1.3: After the initial geometric configuration is generated, immediately assign physical properties to it and perform manufacturability constraint verification. This is a key step to ensure that the design scheme has engineering feasibility.

[0085] First, the environmental medium, typically air, is assigned standardized c0 and ρ0; the solid portion of the metamaterial and lens utilizes a typical photosensitive resin material with a density set to... Speed ​​of sound set to Therefore, the equivalent refractive index can be preliminarily estimated.

[0086] Subsequently, the system initiates three hard constraint checks: The first check examines the minimum feature size, ensuring that all geometric details, such as the metamaterial arm width d or thickness h, are not less than the current additive manufacturing process's resolution limit, typically 50 micrometers; the second check examines the volume fraction, calculating the proportion of solid material in the total volume of the acoustic metamaterial array to prevent printing failure or acoustic performance degradation due to overfilling; the third check examines topological connectivity, verifying through image morphology methods whether each unit maintains a physical connection with adjacent units or the substrate, avoiding unsupported suspended structures. Only when all three checks pass is the configuration considered a valid initial value. The formal expression of the constraint criteria is as follows:

[0087] ,

[0088] ,

[0089] ,

[0090] In this formula, Check1, Check2, and Check3 represent the Boolean judgment results of the three manufacturability constraints; d is the feature size of the metamaterial element; h is the thickness of the metamaterial element; min(d, h) represents the smaller of the two, representing the most critical manufacturing feature size; d min The minimum feature size threshold allowed by the process is set to [value]. ; The volume fraction of the acoustic metamaterial array; This represents the total volume of the solid material in the array; The total envelope volume occupied by the metamaterial array; The upper limit for volume fraction is set at 0.6, which is 60%. This is a topological connectivity judgment function, whose input is the geometric parameters of the metamaterial element. The output is a boolean value; true indicates that the structure is connected and there are no isolated components.

[0091] Step 1.4: The complete design data that has passed constraint verification is encapsulated into a structured initialization data package. This data package serves as the final output of this stage. It not only contains geometric information but also integrates material properties, constraint states, and metadata, ensuring that downstream modules do not need to repeatedly parse the original input. Specifically, the geometry section records the definitions of concave surfaces, metamaterial arrays, and lenses in parametric form; the material section clearly distinguishes the acoustic parameters of the environmental medium and structural materials; and the constraint section saves the verification results and the thresholds used for subsequent debugging and traceability. The data package employs a hybrid serialization strategy: geometric and constraint information is stored in lightweight JSON format for fast retrieval; while high-dimensional phase field or mesh data is compressed and stored in HDF5 format, balancing efficiency and accuracy. After encapsulation, the system will... The data is passed to the 3D acoustic computation domain construction module in step 2, serving as the sole input source for its mesh generation, material assignment, and boundary setting, thus establishing an end-to-end data flow loop. Its data structure is defined as follows:

[0092]

[0093] In this formula, This represents a concave geometry defined by the radius of curvature R. This represents the complete initialization data package output from step 1; Geometry is a subset of geometry, containing parameterized descriptions of three substructures; Indicates periodicity Feature dimension d, thickness h, and rotation angle Defined acoustic metamaterial array; Represented by the spatial phase distribution function Defined acoustic lens; Material is a subset of material properties; ρ0 and c0 are the density and sound velocity of the surrounding medium, respectively; ρ m and c m These represent the density and sound velocity of structural materials such as photosensitive resin; Constraints is a subset of constraint information; d min The minimum feature size threshold; The upper limit of volume fraction; ConnFlag is the topology connectivity check flag, with a boolean value. .

[0094] Step 2: Based on the finite element or finite difference method, construct a three-dimensional acoustic computing domain, perform unified modeling of the entire physical hardware structure link, complete global mesh generation and multi-scale material mapping, and provide a unified computing foundation for high-fidelity simulation.

[0095] Specifically, such as Figure 2 As shown, the entire physical hardware structure link is divided into three key sub-regions: concave energy focusing domain. Metamaterial equivalent domain With lens-air interface region This is used to achieve differentiated modeling of acoustic behavior at different scales.

[0096] In the concave energy focusing domain, a fine volumetric mesh is used to precisely subdivide the highly curved concave surface to ensure that its scattering wavefront distortion and energy focusing process are accurately captured.

[0097] In the metamaterial equivalent domain, for periodic acoustic metamaterial arrays, a multi-scale homogenization technique is used to transform them into a continuous medium with equivalent refractive index and impedance. A uniform mesh is applied within this region to avoid computational redundancy caused by microstructure. It should be noted that this multi-scale homogenization is only applied when the metamaterial unit response is within the linear range, i.e., the incident sound pressure amplitude is below the material yield threshold. If the local sound pressure is too high, the system marks this region as a nonlinear hotspot and reduces the unit size or replaces it with a high-pressure-resistant material in subsequent iterations. The multi-scale homogenization technique obtains the equivalent parameters by solving a unit cell problem.

[0098] In the lens-air interface domain, a progressive mesh refinement strategy is adopted to ensure that the boundary layer resolution meets the numerical stability requirements for acoustic wave penetration and reflection. Simultaneously, the system explicitly embeds nonlinear parameters and thermoviscous dissipation coefficients into the nonlinear terms of the Helmholtz equation, forming a transient wave equation containing pressure-related source terms, which is used to simulate harmonic generation and waveform distortion effects under high-intensity ultrasound.

[0099] Finally, all sub-regions are solved in a coupled manner under a unified Cartesian grid, forming a physically consistent, cross-scale collaborative high-fidelity sound field simulation platform, providing reliable input for subsequent sensitivity analysis and optimization iterations.

[0100] More specifically, the system design includes a three-dimensional acoustic computational domain construction module, and the module execution details include the following steps:

[0101] Step 2.1: The system first processes the initialization data packet output in Step 1. Extracting the geometry subset: , and These parametric geometries are uniformly converted into boundary representation B, in rep format, and embedded within a cubic computational domain that encloses the entire acoustic link. In this context, the size of the domain is determined by the reference length. Dynamic scaling is used to ensure the focal point is located in the central region of the domain and sufficiently far from the boundary to avoid reflection interference. Subsequently, the system performs Boolean operations to separate each solid structure from the background medium, forming a multi-regional geometric topology. Indicates the air medium area. This represents the solid region of the metamaterial and lens. The resulting closed, multi-connected computational domain contains clear internal interfaces, providing a topological basis for subsequent mesh generation.

[0102] ,

[0103] in, The spatial extent of the external computational domain is represented by a cuboid definition; L x L y L z These represent the half-length or full-length of the computational domain in the x, y, and z directions, respectively; x f y f , z f These are the coordinate components of the target focus; the constants 0.1m and 0.15m are safety margins to ensure the focus is far from the boundary. , , These represent the solid regions occupied by the concave surface, metamaterial, and lens, respectively. The final multi-connected computational domain used for simulation is the background air domain minus all solid structures; the symbol \ represents the set difference operation.

[0104] Step 2.2: The system calls a 3D mesh generator based on Delaunay triangulation to generate a mesh for the multi-connected computational domain. Perform unstructured tetrahedral mesh generation. The initial mesh scale is determined by the operating wavelength. It was decided that the global minimum unit size would be set to λ / 12 to meet the acoustic simulation accuracy requirement of at least 12 units per wavelength.

[0105] Subsequently, adaptive local encryption is performed in three key regions: one is in the target focal neighborhood: radius Within the sphere, the density is increased to λ / 20; secondly, in the gaps between metamaterial units: the width is less than... The region is encrypted to the minimum feature size d. min 1 / 2; thirdly, near the solid-air interface: thickness The shell is encrypted to accurately capture impedance abrupt changes.

[0106] After the mesh is generated, the system calculates the quality indicators of each cell, such as skewness and aspect ratio, and discards cells with quality below a threshold q. min Defective cells with a density of 0.3 are identified and re-meshed until the overall pass rate exceeds 98%. The mesh generation rules are as follows:

[0107]

[0108] in, Indicates the global default grid size; The size of the encrypted grid in the focus area; This refers to the local mesh size at the gaps in the metamaterial microstructure; X represents the mesh size of the densified layer near the solid-gas interface; X represents the coordinates of any point in space. Indicates the distance from point X to the boundary of the solid region. The Euclidean distance; d is the operating wavelength; min This represents the minimum manufacturing feature size; all h values ​​refer to the upper limit of the maximum side length of the tetrahedral element.

[0109] Step 2.3: The system traverses all meshes and automatically assigns corresponding material physics parameters based on the geometric region to which the centroid belongs. Specifically, if the centroid is located in the air domain... Then, density ρ0 and sound velocity c0 are assigned; if it is located in a solid structural domain This assigns the material density ρ m and speed of sound c m For the acoustic lens region, although geometrically it belongs to... However, its equivalent acoustic response is determined by the phase distribution. Implicitly defined, it is therefore treated as a homogeneous solid in this step; its fine-tuning will be reflected in the boundary conditions in step 3. Material mapping is accelerated by spatial indexing, ensuring that the assignment operation of millions of cells is completed in seconds. The assignment rule is expressed as:

[0110] ,

[0111] In this formula, ρ(X) represents the local density at spatial location X; c(X) represents the local sound velocity. The region is an air medium. For the union of all solid structural domains, including concave surfaces, metamaterials, and lenses; ρ0 and c0 are the environmental parameters from step 1, normalized: density and sound velocity; ρ m and c m The preset structural material parameters, such as photosensitive resin, are used; this piecewise function realizes the spatial discretization mapping of material properties.

[0112] Step 2.4: The system applies three types of physical boundary conditions to the boundary of the computational domain: First, a circular region is designated as the ultrasonic incident surface on the back side of the concave energy-concentrating structure, and a plane wave or focusing source excitation with a frequency of f is applied; Second, on the outer surface of the computational domain... Set the PML absorption boundary for a perfect match layer with a thickness of [missing information]. First, to simulate an infinite open space and suppress false reflections; second, on a symmetric plane, such as the y=0 plane, if geometrically symmetric, apply symmetric boundary conditions to reduce computational load.

[0113] After boundary settings are completed, the system performs final mesh quality verification: statistical cell skewness distribution is analyzed to ensure maximum skewness s. max <0.85, mean skewness It also checks if the Jacobian determinant is all positive to prevent mesh flipping. After successful verification, the computational domain is constructed, and the output is in a format readable by the finite element solver, such as .inp or .msh.

[0114] The above-mentioned boundary and quality control are expressed as follows:

[0115] The circular region on the back side of the parabola, with radius r src =0.8R,

[0116] PML physical thickness: ,

[0117] Network quality check: , ,

[0118] In this formula, Indicates the ultrasonic incident surface; The radius of the excitation source; The radius of curvature of the concave surface; To perfectly match the physical thickness of the layer; The skewness of the i-th grid cell, with a value in the range [0,1], is better the smaller it is; s max The maximum skewness; N e This represents the total number of elements; all metrics are used to evaluate the applicability of the mesh to numerical solutions of the sound field.

[0119] In summary, this step creates a three-dimensional mesh structure for the acoustic computational domain based on the geometric and material parameters in the data package generated in step 1, setting corresponding boundary types and material properties for different substructures. To ensure the computational accuracy and stability of the model, refined meshing is used in the concave energy-focusing domain and lens region, while adaptive sparse meshing is used in the far-field or air domain to reduce computational load. Subsequently, a multi-scale homogenization technique is introduced for the metamaterial equivalent domain, transforming the subwavelength-scale periodic microstructures into a continuous medium with effective mass density and effective bulk modulus through the equivalent medium method, thereby significantly reducing the overall computational load without sacrificing the realism of acoustic features. This process is achieved through periodic boundary modeling and feature parameter extraction of local units, allowing the acoustic metamaterial microstructures to participate in the sound field solution in the macroscopic computational domain as equivalent parameters. In this way, a physical mapping from micro-units to the macroscopic sound field is completed in the three-dimensional computational domain, establishing a complete computational foundation for the subsequent numerical solution of nonlinear acoustic equations and sound field distribution analysis.

[0120] Step 3: Perform nonlinear sound field simulation to evaluate the performance of the current configuration in terms of focal energy concentration and sidelobe suppression, and determine whether it meets the preset acoustic indicators to decide whether to proceed to the optimization process.

[0121] Specifically, after completing the construction of the three-dimensional acoustic computational domain, a high-fidelity nonlinear sound field simulation module was designed and launched. Nonlinear terms were introduced into the sound field control equations to accurately characterize the complex physical behavior of high-intensity ultrasound propagation in air. The Westervelt nonlinear acoustic equations were adopted as the core control model, comprehensively considering key factors such as the weak nonlinear effects of the medium, thermoviscous energy dissipation, and gradually varying waveform envelopes. This model is suitable for describing the three-dimensional sound propagation process, including focal convergence, harmonic generation, and wavefront distortion. The Westervelt nonlinear acoustic equations are used to describe the propagation behavior of high-intensity ultrasound, and the equations are as follows:

[0122] ,

[0123] in Sound pressure level, measured in Pascals (Pa), is a spatial unit. and time The function; For the Laplace operator acting on spatial coordinates; ρ is the time variable, in seconds (s); c0 is the velocity of a small signal in the medium, in m / s; ρ0 is the static density of the medium, in kg / m³; β is a dimensionless parameter characterizing the nonlinear properties of the medium. The thermoviscosity dissipation coefficient, expressed in m² / s, is determined by both the shear viscosity and bulk viscosity of the medium, and satisfies the following conditions: This equation is applicable to three-dimensional sound fields with weak nonlinearity, slowly varying envelopes, and thermal viscous energy dissipation.

[0124] In practice, the corresponding excitation source conditions are set according to the input ultrasonic working frequency, and the control equations are implicitly solved in the time domain in combination with the established geometric and material parameters. The dynamic evolution characteristics of the sound pressure field, energy density distribution and focal region are tracked simultaneously. To ensure the accuracy and stability of the calculation, adaptive time step control and residual convergence criteria are adopted in the solution process.

[0125] Furthermore, a local field density output mechanism is enabled in the near-field region of the focal point to accurately capture the main lobe peak, side lobe structure, and energy concentration characteristics at the subwavelength scale.

[0126] The final output of global sound field response data will serve as the physical basis for subsequent sensitivity analysis and closed-loop optimization, fully supporting the system process from modeling real acoustic behavior to performance-driven design.

[0127] More specifically, the execution details of the system's high-fidelity nonlinear sound field simulation module include the following steps:

[0128] Step 3.1: The system calls the frequency domain finite element solver to solve the multi-connected computational domain constructed in step 2. The modified Westervelt governing equations are solved and coupled with a thermoviscous boundary layer model to accurately capture the propagation, focusing, and dissipation behavior of high-amplitude ultrasound in complex media.

[0129] The Westervelt governing equations, in frequency domain form, are expressed as relating complex sound pressure levels. The second-order partial differential equation, which includes the nonlinear term β, the complex equivalent wavenumber k introduced by the thermal conductivity coefficient k and the dynamic viscosity μ. eff The nonlinear parameter β is set according to the medium type: β=1.2 for air, β=3.5 for water, and β=[missing value] for polymer materials. When the sound pressure amplitude The nonlinear equations are automatically enabled when the condition is met; otherwise, they degenerate into linear wave equations.

[0130] The solver employs a direct sparse matrix solution method, such as MUMPS, and sets a convergence tolerance. And enable parallel computing acceleration.

[0131] Excitation source at the ultrasonic incident surface The applied amplitude is The single-frequency harmonic pressure load has a frequency of f. After solving, the full-field complex sound pressure level is output. particle velocity and time-averaged sound intensity vector Its governing equations and output are defined as follows:

[0132] ,

[0133] In this formula, Indicates spatial location Complex sound pressure at a location, unit: Pa; For the Laplace operator; k eff is the complex equivalent wavenumber, used to characterize attenuation and dispersion in a medium; i is the imaginary unit; It is the thermal viscosity absorption coefficient; For dynamic viscosity, by It is derived that; The period of the sound wave; The velocity vector of the particle, in m / s; Angular frequency; The local density is derived from step 2.3; express The complex conjugate; This is the time-averaged sound intensity vector, in W / m². The equations represent taking the real part of a complex number; this system of equations fully describes the propagation of a sound field with losses, negligible nonlinearity, but retaining thermoviscous effects.

[0134] It should be noted that the nonlinearity is mainly reflected in the subsequent heat source term through β. The complete nonlinear form of the Westervelt equation needs to be solved iteratively at high amplitudes. However, for simplification, this system adopts a linear thermoviscosity model in the initial evaluation stage, and the nonlinear thermal deposition is modeled separately in step 3.3.

[0135] Step 3.2: The system is in the target focal neighborhood, sphere ,radius Internal relative sound intensity amplitude Post-processing analysis was performed.

[0136] First, the actual peak sound intensity point is located by searching for local maxima. Subsequently, sound intensity profiles were extracted along three orthogonal directions (x, y, z), and the full width at half maximum (FWHM) was calculated as the lateral and axial resolution indicators. Next, the first local maximum outside the main lobe was identified, and its ratio to the main lobe peak value was calculated to obtain the sidelobe suppression ratio (SSR). Finally, the volume within the focal region where the sound intensity exceeded 50% of the peak value was statistically analyzed. This serves as a proxy indicator for the effective treatment volume. The key feature extraction formula is as follows:

[0137] ,

[0138] in, Indicates the actual peak sound intensity location; Indicates the target focus Centered on, with radius Evaluation sphere; The amplitude of the sound intensity; The full width at half maximum (FWHM) is the height along the d-axis. and These are the two coordinate positions when the sound intensity profile drops to half of its peak value; Peak sound intensity; The maximum local sound intensity at the first sidelobe; SSR is the sidelobe suppression ratio, expressed in decibels (dB). The volume of space where the sound intensity exceeds 50% of the peak value; This is an indicator function that takes the value 1 if the condition is true, and 0 otherwise. This represents a volume element.

[0139] Step 3.3: To assess the biosafety of high-energy focused ultrasound (HIFU) applications, the system calculates the steady-state thermal deposition distribution based on the acoustic field and estimates the local temperature rise.

[0140] Assuming the power of the heat source It is proportional to the sound intensity, and the proportionality coefficient is determined by the absorption characteristics of the medium. A simplified steady-state heat conduction equation is used. Solving the temperature field thermal conductivity Let it be a constant, such as air. The heat source model is as follows:

[0141] ,

[0142] In this formula, This represents the power of a heat source per unit volume, with the unit being W / m³. Local sound absorption coefficient, unit: Np / m; The air absorption coefficient is given by... and Calculated; The absorption coefficient of a solid material is typically much larger than that of a solid material. ; Steady-state temperature field, unit: K; The ambient temperature is 293 K by default. This is the maximum temperature rise; this assessment is used to determine whether there is a risk of overheating, such as... >50K is considered dangerous.

[0143] Step 3.4: The system compares the extracted performance indicators with the preset thresholds item by item, generating a Boolean compliance vector. .

[0144] The setting criteria include: focus offset Lateral focus resolution Sidelobe suppression ratio Maximum temperature rise If all If the design is successful, proceed to step 4; otherwise, the system records the failure dimension and sets the standardized parameters accordingly. With error signal The result is passed to the optimization module, triggering parameter adjustment iterations. The decision logic is as follows:

[0145] ,

[0146] In this formula, This indicates whether the focus positioning accuracy meets the standard; This indicates whether the lateral focusing resolution meets the requirements; This indicates whether the sidelobe suppression capability is sufficient; This indicates whether thermal safety is controllable; This is an indicator function that outputs true or 1 when the condition is met, and false or 0 otherwise. The operating wavelength; and These are the full width and height at half maximum (FWHM) in the x and y directions, respectively. The unit is K; this decision mechanism ensures that only designs that simultaneously meet acoustic performance and safety constraints are accepted.

[0147] Step 4: For different interface types in the physical hardware structure link, automatically match the boundary element method and the transfer matrix method to achieve high-precision hybrid numerical modeling of strong curvature scattering interface and multi-layer planar transmission interface.

[0148] like Figure 3 As shown, the system utilizes a heterogeneous interface hybrid numerical solution module, which mainly consists of a concave energy-concentrating structure solution unit, a metamaterial equivalent domain solution unit, an acoustic lens layer solution unit, an interface matching unit, and an error evaluation unit. These units operate collaboratively within the same computational link, achieving unified modeling of the cross-domain sound field.

[0149] First, for the curved interface between the concave energy-concentrating structure and the air, a solution element for the concave energy-concentrating structure is established in the concave energy-concentrating domain. This element uses the boundary element method (BEM) to solve the frequency domain Helmholtz equation to calculate the scattered sound field and extract the sound pressure boundary data as input to the interface matching element. The Helmholtz equation is solved in the frequency domain using the boundary element method, through discretization of the boundary integral form. The scattered sound field distribution in this region is obtained to accurately reflect the wavefront distortion and energy redistribution caused by surface focusing. For free space Green's function, For the boundary surface, Represents the normal derivative.

[0150] Subsequently, the metamaterial equivalent domain solution element is invoked in the metamaterial equivalent domain. This element adopts the finite element method (FEM) and uses the nonlinear weak form expression of the sound field as the governing equation to solve the spatiotemporal distribution of sound pressure under a multi-scale grid, thereby obtaining high-precision main field results that include nonlinear propagation effects.

[0151] For the approximately planar or layered interface between the acoustic metamaterial array and the acoustic lens, an acoustic lens layer solution unit is set in the lens-air interface domain. This unit uses the Transmission Matrix Method (TMM) to construct a chain product of characteristic matrices based on the equivalent acoustic impedance and thickness of each layer, calculates the transmission and reflection coefficients, and performs layered transmission, phase compensation and wavefront shaping on the transmitted wave to make the energy distribution in the focusing area meet the design target.

[0152] The interface matching unit establishes unified boundary conditions between different physical algorithms to ensure the continuity of sound pressure and normal velocity at the heterogeneous interface, and to ensure seamless connection of calculation results between the concave energy-concentrating structure solution unit, the metamaterial equivalent domain solution unit, and the acoustic lens layer solution unit.

[0153] The error assessment unit performs differential calculations on the integrated sound field solution output by the interface matching unit, extracts the energy comparison information between the focal and side lobe regions, and outputs the error amount δ and the sound field distribution results, providing input for the dual-channel sensitivity analysis module.

[0154] To ensure the conservation of acoustic energy and the continuity of field quantities throughout the entire process, a unified boundary condition coupling strategy is implemented for the boundary element method and the transfer matrix method in the frequency domain. This strategy includes pressure continuity and normal particle velocity matching constraints, thereby achieving high-fidelity and conflict-free acoustic field connection in hybrid interface scenarios. The interface processing results are seamlessly embedded into the global acoustic field solution framework as boundary excitations or internal source terms of the control equations, further enhancing the overall simulation's ability to reproduce real physical processes.

[0155] More specifically, the execution details of the heterogeneous interface hybrid numerical solution module include the following steps:

[0156] Step 4.1: After completing the global geometric modeling, automatically scan all material interface regions in the physical chain, based on the local radius of curvature. Interface type is identified by the rate of change of interface normal: if the following conditions are met... Furthermore, if the curvature continuity is high, such as the interface between a concave energy-concentrating structure and air, it is classified as a highly curvature surface interface, and the Boundary Element Method (BEM) is used. If the interface is approximately planar or layered, such as between a metamaterial array and a gradient refractive index lens, and the adjacent layers have uniform thickness, it is classified as a multilayer planar interface, and the Transfer Matrix Method (TMM) is used. This classification result is written into the interface attribute label for subsequent solver scheduling.

[0157] Step 4.2: For the highly curved interface between the concave energy-concentrating structure and the air, the system constructs the boundary integral form of the Helmholtz equation in the frequency domain and discretizes the surface using quadratic isoparametric elements. The pressure at the interface is obtained by solving the Fredholm second-type integral equation. With normal velocity The distribution is then used to calculate the scattered sound field. This method accurately captures wavefront distortion, focal spot broadening, and sidelobe enhancement effects caused by Gaussian curvature, avoiding numerical dispersion errors in the volume mesh method in thin air domains.

[0158] Step 4.3: For the near-planar interface between the metamaterial and the lens, the Transmission Matrix Method (TMM) is used. The system treats each layer as a homogeneous equivalent medium and extracts its equivalent acoustic impedance. Based on the thickness, construct a single-layer feature matrix for the i-th layer of the medium:

[0159] The total transfer matrix is ​​a chain product of the matrices at each layer. Therefore, the overall transmittance coefficient can be directly calculated analytically. Reflection coefficient and cumulative phase delay Its efficiency is far higher than that of full-wave simulation.

[0160] In this formula, This represents the characteristic transmission matrix of the i-th layer of medium; h is the wave number of that layer. i For layer thickness; Z i The equivalent acoustic impedance is defined by the product of density and sound velocity; this matrix form is the key implementation method of this invention in the rapid modeling of multilayer acoustic interfaces, used to efficiently obtain phase and energy transfer characteristics.

[0161] Step 4.4: To ensure a conflict-free connection between the curved and planar interfaces in the global sound field, the system implements unified coupling between the BEM and TMM regions in the frequency domain: a pressure continuity is forcibly applied at the boundary between the two methods. Matching the velocity of the normal particle Constraints. By introducing Lagrange multipliers or penalty functions, these interface conditions are embedded in the global linear system to ensure that energy conservation and field quantities are connected without conflict.

[0162] Finally, the interface simulation results are injected into the global sound field solver in the form of equivalent source terms or Neumann / Dirichlet boundary conditions to improve the overall simulation fidelity.

[0163] Step 5: Construct a dual-channel sensitivity analysis architecture for the primary and secondary fields to efficiently extract the sensitivity gradient of the design variables with respect to the objective function, providing a precise mathematical basis for parameter updates.

[0164] Specifically, a dual-channel computational architecture consisting of a primary field solver and an adjoint field solver is constructed. The primary field solver, based on the three-dimensional nonlinear acoustic computational domain and heterogeneous interface hybrid model established in steps 2 and 4, performs high-fidelity full-wave simulation to calculate the sound pressure distribution in the spatiotemporal domain. This provides a complete description of the propagation and convergence process of ultrasound waves in a three-stage focusing structure.

[0165] The adjoint field solver solves for sensitivity information under the backpropagation condition of a given objective function, obtains the gradient information of each design variable in the physical link, including the concave surface curvature radius, metamaterial unit geometric parameters and lens refractive index distribution, and evaluates the gradient influence of each design variable in the physical link on the acoustic pressure amplitude or energy density in the focal region.

[0166] The dual-channel architecture achieves efficient sensitivity analysis while ensuring physical consistency by sharing the same computational domain mesh and material parameters, avoiding the huge computational overhead caused by multiple forward simulations in the traditional finite difference method.

[0167] Furthermore, after each full-wave simulation of the main field is completed, the system calculates the true acoustic response vector of each element in the metamaterial array at the operating frequency based on its local incident and transmitted acoustic field data using a parametric inversion method. The gradient results serve as the benchmark for evaluating the accuracy of deep neural network surrogate models; the resulting gradients are directly used to guide subsequent optimization iterations, ensuring that the design update direction always points to performance improvement, and providing a precise mathematical basis for the closed-loop feedback mechanism.

[0168] More specifically, the dual-channel accompanying sensitivity analysis architecture includes the following steps in its execution:

[0169] Step 5.1: Based on the computational domain from Step 2 and the interface processing results from Step 4, the system constructs a frequency domain sound pressure field solver and outputs the complex sound pressure distribution. This solver fully simulates the entire process of ultrasonic wave propagation, diffraction, and focusing in a three-stage focusing structure: concave reflector → metamaterial modulation layer → acoustic lens, providing the fundamental field quantities for calculating the main objective function.

[0170] Step 5.2: To efficiently obtain the sensitivity of design variables to the target performance, the system simultaneously constructs an adjoint field solver. Guided by the objective function defined in Step 7.1, this solver inversely maps the energy distributions of the focal and sidelobe regions to equivalent source terms, solving the adjoint Helmholtz equations, which have the same form but opposite physical meaning. The system defines the objective function. ,in As the focal area, This is the side lobe region. These are the sidelobe suppression weighting coefficients. Based on this, the adjoint Helmholtz equation in the frequency domain is derived, and its expression is:

[0171]

[0172] in To accompany the sound pressure field, For wave number, and These are the indicator functions for the focal region and the sidelobe region, respectively. Accompanying sound pressure field. Physical consistency is ensured by solving in reverse using the same mesh.

[0173] Using the adjoint method theory, any design variable Such as the radius of curvature of a concave surface Metamaterial unit width d, local refractive index of lens The gradient of the objective function can be expressed as relying solely on the inner product of the primary and adjoint fields in the perturbation region, eliminating the need for repeated forward simulations. The system automatically identifies the geometric or material perturbation operators corresponding to each variable and calculates the sensitivity.

[0174] ,

[0175] in To and Associated sensitivity operators, such as shape derivative or material derivative.

[0176] In this formula, This represents the partial derivative of the objective function with respect to the j-th design variable; Home stadium sound pressure; For variables The corresponding perturbation operators are, for example, the boundary deformation derivative with respect to geometric variables and the dielectric / acoustic parameter derivative with respect to material variables; For variables The affected local areas, such as the boundaries of metamaterial units or lens voxels; This represents taking the real part; this gradient expression is the core mathematical tool for achieving efficient optimization in this patent, significantly reducing computational costs.

[0177] Step 5.3: To ensure that step 7 can correctly parse and utilize the sensitivity information, the system calculates the gradient vectors of all design variables. The data is packaged with the current performance metrics and used as input to the optimization feedback layer of the system of this invention. The packaged data structure includes variable names, current values, gradient values, and physical units, ensuring that the subsequent optimizer can correctly parse it and use it for updating direction calculations.

[0178] Step 5.4: For metamaterial microstructures, based on their local incident and transmitted acoustic field data, the parametric inversion method is used to calculate their true acoustic response at the operating frequency.

[0179] Step 6: Introduce a deep neural network model to quickly predict the equivalent acoustic response of the metamaterial array unit, and ensure its physical consistency throughout the optimization process through error monitoring and local correction mechanisms.

[0180] like Figure 4 As shown, the deep neural network model mainly consists of an input layer, a DNN proxy model body, and an output layer. The DNN proxy model body includes a feature extraction block, an activation layer, a pooling layer, and a fully connected layer connected in sequence. Each module works collaboratively in the same data stream to achieve an end-to-end mapping from the manufacturable geometric parameters of metamaterials to the physically equivalent acoustic response.

[0181] The system first receives the set of geometric parameters of the acoustic metamaterial microstructure at the input layer. As the initial input, where d represents the minimum feature size, satisfying h represents the element thickness, a represents the period length, Indicates the material filling rate, which satisfies This parameter set fully characterizes the manufacturable geometry of the unit. The system then feeds this parameter set into the DNN surrogate model for forward propagation. The feature extraction block performs a nonlinear transformation on the input parameters to mine higher-order geometry-acoustic correlation features. The activation layer uses the GeLU function to introduce nonlinear expressive power. The pooling layer retains key features and suppresses redundant information through dimensionality reduction. The fully connected layer fuses all intermediate features to generate a high-dimensional hidden state. Based on this, the system outputs an equivalent acoustic response vector through the output layer. ,in For equivalent refractive index, For equivalent acoustic impedance, The frequency-dependent phase delay is used to characterize the macroscopic acoustic behavior of metamaterials at the target operating frequency.

[0182] Simultaneously, the system integrates an error monitoring and local correction mechanism, which receives prediction results from the output layer. And compare it with the true equivalent response obtained from the high-fidelity full-wave simulation inversion in step 5. Compare and calculate the relative error. ,judge If the value is greater than 5%, the model remains unchanged and the current round of verification ends. If it is, the local parameter correction mechanism is triggered: the system freezes the weight parameters of the first three layers in the main body of the DNN surrogate model to retain the learned general mapping knowledge and avoid catastrophic forgetting, and only performs fine-tuning operations on the weights of the last two layers, correcting their parameters through single-step gradient updates to reduce prediction bias. Finally, the system reconnects the corrected surrogate model to the optimization loop, providing high-precision initial values ​​for the optimizer in step 7, thereby continuously ensuring prediction performance and physical consistency without participating in the calculation of the adjoint field gradient, supporting the efficient and stable operation of the entire "simulation-sensitivity-optimization-feedback" closed loop.

[0183] In one embodiment, the execution details of the deep neural network proxy model include the following steps:

[0184] Step 6.1: To alleviate the computational burden of repeated modeling of metamaterial elements in the high-fidelity full-wave simulation in Step 5, the system introduces a deep neural network as a fast prediction agent, with the input being the normalized geometric parameter vector of the metamaterial element. , representing the feature size ratio, thickness ratio, rotation angle or radians, and fill rate, respectively; the output is the equivalent acoustic response vector. The network activation function is Swish, and the loss function is mean squared error (MSE).

[0185] Step 6.2: The surrogate model is pre-trained using 5000 sets of offline-generated full-wave simulation datasets before optimization begins. At the beginning of each optimization iteration, the system calls this model to quickly predict the macroscopic acoustic response of candidate metamaterial configurations, providing high-quality initial guesses for the optimizer in Step 7 and avoiding blind search.

[0186] Step 6.3: The surrogate model significantly accelerates the design evaluation. Its prediction accuracy needs to be continuously controlled during the optimization process. After each main simulation is completed in step 5, the system extracts the true equivalent parameters of the metamaterial region through inversion or field averaging. , and the surrogate model prediction Comparison. If the relative error of any output component exceeds 5%, i.e. If this occurs, a local correction mechanism is triggered: the sample is added to the training set, and the network undergoes a round of fine-tuning with a learning rate set to [value missing]. Only the weights of the last two layers are updated to prevent catastrophic forgetting. This DNN is explicitly limited to not participating in the adjoint gradient calculation, serving only as an acceleration tool to provide the initial configuration and rapid evaluation. The calculation of actual performance metrics, evaluation of the objective function, and sensitivity analysis all strictly rely on the high-fidelity full-wave simulation results from step 5 to ensure the physical accuracy and reliability of the optimization direction.

[0187] Step 7: Define a weighted objective function to drive closed-loop optimization iteration, combine manufacturability constraint projection and convergence criteria, dynamically update design parameters and output the optimal device configuration that can be used for additive manufacturing.

[0188] Specifically, the weighted objective function is defined with the weighted difference between the energy density of the focal region and the energy of the sidelobe region as the optimization guide, wherein the focal energy weight coefficient and the sidelobe suppression weight coefficient can be dynamically adjusted according to the specific application scenario;

[0189] The optimizer integrates the design variable sensitivity information obtained in step 5 and the metamaterial performance prediction results provided in step 6 to dynamically update the geometric and material parameters in the physical link. After each parameter update, the constraint projector immediately performs manufacturability verification, including: eliminating isolated structures and filling structures smaller than a certain size. The gaps and smoothed sharp edges were eliminated, and morphological closing operations were used to ensure manufacturing feasibility. This ensured that the new configuration met engineering constraints such as a minimum manufacturing feature size of not less than 50 micrometers, an overall integral number of not more than 60%, and structural connectivity.

[0190] The convergence monitor synchronously evaluates the optimization process. When the relative change rate of the objective function between two consecutive iterations is less than 1%, or the cumulative number of iterations reaches the preset upper limit of 20, the optimization process is determined to have converged. If the convergence condition is not met, the process returns to step 2 to rebuild the computational domain and continue iterating. If convergence has been achieved, the final optimal device configuration is output, and a three-dimensional model file that can be used for additive manufacturing is generated. This directly drives the physical link manufacturing process, resulting in a completely passive device structure that requires no external power supply or active control components. It can be directly pre-coupled to the receiving port of a fiber optic Fabry-Perot ultrasonic sensor, fiber optic grating sensor, or MEMS microphone to improve its receiving sensitivity and signal-to-noise ratio for weak ultrasonic signals (frequency range 20 kHz–500 kHz) in the air.

[0191] More specifically, this step includes the following steps:

[0192] Step 7.1: Building upon the sound pressure field and accompanying field sensitivity analysis results provided by the full-wave simulation in Step 5, the system here explicitly defines the objective function used to drive optimization. The objective function is based on the weighted difference between maximizing focal energy and minimizing sidelobe energy, allowing users to dynamically adjust the weight ratios according to the application scenario (e.g., treatment vs. imaging). This flexible definition mechanism enables the same optimization framework to be adapted to various acoustic tasks, forming a guiding benchmark for closed-loop feedback. The system defines the objective function as follows: In this formula, This represents the objective function value for the overall acoustic performance of the current design configuration, used to guide the optimization direction; Indicates the focal area The total acoustic energy within is defined as ,in This refers to local sound intensity. Indicates the side lobe region The total acoustic energy within is defined as ; This is the focus energy weighting coefficient, with a default value of 1, which can be adjusted upwards according to the need for high focusing gain. This is the sidelobe suppression weighting coefficient, with a default value of 0.3, which can be increased when strong sidelobe suppression is needed, for example, to avoid heating of non-target areas. This objective function explicitly balances the two design objectives of "energy concentration" and "spurious emission suppression" through a weighted difference, and its gradient is obtained from step 5.2 based on the main field. With accompanying field Highly efficient computation, ensuring that the optimization process is both physically accurate and computationally efficient, is the core guiding principle of this patent's closed-loop optimization mechanism.

[0193] Step 7.2: Use L-BFGS or Adam optimizer to synthesize the gradient information provided in step 5. Compared with the surrogate model prediction in step 6, for design variables Update:

[0194] In this formula, This represents the design variable vector at the kth iteration, which includes concave curvature, metamaterial geometric parameters, lens refractive index distribution, etc. This represents the updated design variable for the (k+1)th round. Describe the objective function For design variables The gradient vector is provided by the adjoint field solver in step 5; The learning rate represents the optimization step size in the k-th round, which can be determined by line search or an adaptive strategy. This represents the approximate Hessian matrix of the objective function in the k-th round, used to capture the coupling curvature information between variables; Its inverse matrix is ​​efficiently approximated in quasi-Newton methods such as L-BFGS through low-rank updates of historical gradients, avoiding explicit computation. This update formula constitutes the core iterative rule of the closed-loop optimization in this patent, ensuring that the parameters converge efficiently along the performance improvement direction.

[0195] Step 7.3: To prevent the optimization process from generating unmanufacturable structures, the system performs hard constraint projection immediately after each parameter update. The constraint projector performs three hard constraint checks on the new parameters: 1. Minimum feature size II. Total Integral Score 3. Metamaterial unit connectivity (no isolated suspended bodies). If any constraint is violated, the parameters are pulled back to the feasible region by projection operations, for example, truncating the excessively small (d) to 50 μm and scaling other dimensions proportionally to preserve phase properties.

[0196] Step 7.4: The system continuously monitors the improvement rate of the objective function, and the convergence monitor checks the relative rate of change of the objective function in each round: ,like or number of iterations If the condition is met, convergence is determined. The system outputs the final 3D geometric model, material distribution map, and process parameter package, which are directly used to drive additive manufacturing equipment, completing the closed loop from digital design to physical realization. This represents the objective function value in the k-th iteration. The relative rate of change is 1%; the preset convergence threshold is 1%; and the maximum number of iterations is 20. This decision mechanism ensures that the optimization terminates within a finite number of steps, avoiding infinite loops.

[0197] The method of this invention is deployed as an algorithm feedback link on a local workstation or cloud server, supporting linkage with 3D printing equipment. The final output optimal configuration is directly converted into STL or STEP format manufacturing files, realizing an integrated "design-simulation-manufacturing" process.

[0198] In summary, this invention designs a physical hardware structure link and an algorithm coupling feedback collaborative method. The physical hardware structure link sequentially includes a concave energy-concentrating structure, an array of acoustic metamaterial microstructure units, and an acoustic lens, achieving three-level passive focusing and phase correction of ultrasound. The algorithm coupling feedback collaborative method constructs a high-fidelity sound field simulation model, innovatively integrating multi-scale homogenization technology to efficiently characterize the effective acoustic parameters of the composite metamaterial. It introduces a nonlinear pressure acoustic transient equation (including a nonlinear parameter β) to accurately describe harmonic generation and waveform distortion in high-intensity ultrasound propagation, and uses the transfer matrix method or boundary element method to accurately simulate the transmission, reflection, and phase shift behavior at the heterogeneous interface. Furthermore, the system constructs a deep neural network model, using the geometric parameters of the metamaterial units as input and the equivalent acoustic response as output, to achieve real-time performance prediction and dynamic impedance matching optimization. Sensitivity information is obtained through dual-channel solution of the main field and adjoint field. Combined with objective function evaluation and constraint projection mechanisms, a closed loop of "simulation-sensitivity-optimization-feedback" is formed, iteratively outputting the optimal manufacturable configuration. This device can achieve a sound pressure gain of over 15 dB and a focal size of less than half a wavelength in the 20–500 kHz frequency band. Its high-performance focusing effect significantly improves the receiving sensitivity and signal-to-noise ratio of fiber optic ultrasonic sensors. It has advantages such as being passive, low-cost, and highly adaptable, and is suitable for complex industrial application scenarios such as partial discharge monitoring of power equipment and fault diagnosis of rail transit.

[0199] The above description is merely a preferred embodiment of the present invention and is not intended to limit the invention in any way. Although the implementation process of the present invention has been described in detail above, those skilled in the art can still modify the technical solutions described in the foregoing examples or make equivalent substitutions for some of the technical features. All modifications and equivalent substitutions made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A method for modeling and optimizing an ultrasonic enhancement device, characterized in that, Includes the following steps: Step 1: Analyze and process the acoustic field signal of the operating equipment, design the initial geometric configuration of the ultrasonic enhancement device, and generate a parameter data package of the initial geometric configuration. The structure of the initial geometric configuration includes a concave energy-concentrating structure, a metamaterial microstructure, and a lens structure. Step 2: Based on the parameter data package, create a three-dimensional acoustic calculation domain and set the corresponding boundary conditions and material properties for each structure; Step 3: Perform nonlinear sound field simulation to evaluate the performance of the current geometric configuration, in order to determine whether the current geometric configuration should be optimized; Step 4: For different interface types of the current geometric configuration, perform heterogeneous interface hybrid modeling and solve to obtain interface simulation results; Step 5: Based on the boundary conditions, material properties, and interface simulation results, perform full-wave simulation on the current geometric configuration, and conduct sensitivity analysis on the full-wave simulation results according to the preset target results; use the parameter inversion method to obtain the true acoustic response of the metamaterial microstructure from the full-wave simulation results; Step 6: Introduce a deep neural network surrogate model, take the geometric parameters of the metamaterial microstructure as input, quickly predict the equivalent acoustic response, and use the real acoustic response to compare the error, thereby triggering a local correction mechanism to optimize and update the surrogate model; Step 7: Combining the sensitivity analysis results from Step 5 and the predicted equivalent acoustic response results from Step 6, and using an objective function evaluation and constrained projection mechanism, dynamically update the design parameters of the current geometric configuration and output the optimal device configuration that can be used for additive manufacturing.

2. The method for modeling and optimizing an ultrasonic enhancement device as described in claim 1, characterized in that, In step 1, the parameters of the device's acoustic field signal include the ultrasonic operating frequency, the target focal position, and the medium properties; The concave energy-concentrating structure is used to achieve primary energy convergence of ultrasound waves. This structure adopts a parabolic or spherical geometric configuration. The acoustic metamaterial microstructure is used for mesoscale phase compensation and impedance matching. The structure adopts a cross-shaped, H-shaped, labyrinth-shaped or spiral topology. The acoustic lens structure achieves wavefront shaping and focusing through refractive index distribution. This structure adopts the form of a gradient refractive index lens or a planar superlens.

3. The method for modeling and optimizing an ultrasonic enhancement device as described in claim 1, characterized in that, In step 1, the manufacturability constraints of the initial geometric configuration of the design are checked, specifically: the minimum feature size is checked; the volume fraction is checked and the proportion of solid material in the total volume of the metamaterial microstructure is calculated; and the topological connectivity is checked and the physical connections of the metamaterial microstructure are verified by image morphology methods.

4. The method for modeling and optimizing an ultrasonic enhancement device as described in claim 1, characterized in that, Step 2 specifically includes: Step 2.1: Extract the geometric parameter information of the concave energy-concentrating structure, metamaterial microstructure and lens structure from the parameter data package, convert them into boundary representations and embed them into a cubic computational domain. Then perform Boolean operations to separate each solid structure from the background medium, forming an air medium region and a solid region, and finally generate a multi-connected computational domain. Step 2.2: Call the 3D mesh generator based on Delaunay triangulation to perform unstructured tetrahedral meshing on the multi-connected computational domain; among which, adaptive local refinement is performed in three key regions: first, in the neighborhood of the target focus; second, in the gaps between metamaterial elements; and third, near the solid-air interface. Step 2.3: Traverse all meshes and automatically assign corresponding material parameters based on the geometric configuration region to which the centroid of a mesh belongs; Step 2.4: Apply three types of physical boundary conditions to the boundary of the computational domain: First, specify a circular region on the back side of the concave energy-concentrating structure as the ultrasonic incident surface and apply plane wave or focusing source excitation; Second, set a perfectly matched layer absorbing boundary on the outer surface of the cubic computational domain to simulate an infinite open space and suppress false reflections; Third, apply symmetric boundary conditions on the symmetric plane if the geometry is symmetric. After the boundary settings are completed, perform the final mesh quality verification.

5. The method for modeling and optimizing an ultrasonic enhancement device as described in claim 4, characterized in that, The multi-connected computational domain includes a concave energy focusing domain, a metamaterial equivalent domain, and a lens-air interface domain. In the metamaterial equivalent domain, a multi-scale homogenization technique is introduced to transform the metamaterial microstructure and equivalent medium method into a continuous medium with equivalent refractive index and impedance.

6. The method for modeling and optimizing an ultrasonic enhancement device as described in claim 4, characterized in that, Step 3 specifically includes: Step 3.1: Set the corresponding excitation source according to the input ultrasonic working frequency, and use the direct sparse matrix solution method to describe the propagation behavior of high-intensity ultrasonic waves in the multi-connected computational domain using the Westervelt nonlinear acoustic equation; Step 3.2: In the neighborhood of the target focus, firstly, the actual peak sound intensity point is located by searching for local maxima; then, sound intensity profiles are intercepted along the three orthogonal directions (x, y, z), and the full width at half maximum (FWHM) is calculated as the lateral and axial resolution indicators; next, the first local maximum outside the main lobe is identified, and its ratio to the main lobe peak value is calculated to obtain the sidelobe suppression ratio; finally, the volume in the focus region where the sound intensity exceeds 50% of the peak value is counted. Step 3.3: Assume that the heat source power is proportional to the sound intensity, calculate the steady-state heat deposition distribution based on the sound field, and estimate the local temperature rise; Step 3.4: Based on the focus positioning accuracy, lateral focus resolution, sidelobe suppression ratio obtained in Step 3.2 and the maximum temperature rise obtained in Step 3.3, compare them with the preset thresholds one by one. If all comparison results meet the threshold requirements, the geometric configuration design is successful; otherwise, proceed to the optimization process.

7. The method for modeling and optimizing an ultrasonic enhancement device as described in claim 1, characterized in that, Step 4 specifically includes: Step 4.1: Automatically scan all material interface regions. For the curved interface between the concave energy-concentrating structure and the air, use the boundary element method (BEM). For the metamaterial microstructure, use the finite element method (FEM), with the nonlinear weak form expression of the sound field as the governing equation. Solve the spatiotemporal distribution of sound pressure under a multi-scale grid to obtain the main field results including nonlinear propagation effects. For the interface between the metamaterial microstructure and the acoustic lens that is approximately planar or layered, use the transfer matrix method (TMM). Step 4.2: For the highly curved surface interface between the concave energy-concentrating structure and the air, the Boundary Element Method (BEM) is used to solve the frequency domain Helmholtz equation to calculate the scattered sound field and extract the sound pressure boundary data. Step 4.3: For the interface between the metamaterial microstructure and the acoustic lens that is approximately planar or layered, the Transmission Matrix Method (TMM) is used. Each layer is regarded as a homogeneous equivalent medium. The characteristic matrix chain product is extracted and constructed based on the equivalent acoustic impedance and thickness of each layer. The transmission and reflection coefficients and the cumulative phase delay are calculated. The transmitted wave is subjected to layered transmission, phase compensation and wavefront shaping so that the energy distribution in the focusing area meets the design target. Step 4.4: For interfaces with strong curvature surfaces and interfaces with approximate planes or layered stacks, implement unified coupling between the Boundary Element Method (BEM) region and the Transmission Matrix Method (TMM) region in the frequency domain, and establish unified boundary conditions: Force the sound pressure continuity and normal particle velocity matching constraints at the boundary between the regions of the two methods, and embed the above interface conditions in the global linear system by introducing Lagrange multipliers or penalty functions. Finally, differential calculations are performed on the unified boundary conditions to extract the energy comparison information between the focal and sidelobe regions, and the error and sound field distribution results are output.

8. The method for modeling and optimizing an ultrasonic enhancement device as described in claim 7, characterized in that, Step 5 specifically includes: Step 5.1: Construct the master field solver, input the boundary conditions, material properties and interface simulation results, and output the distribution of sound pressure in the spatiotemporal domain; Step 5.2: Simultaneously construct the adjoint field solver, which solves for sensitivity information under the backpropagation condition of the given objective function, and obtains the design variables of the current geometric configuration; Specifically, using the adjoint method theory, the gradient of any design variable with respect to the objective function can be expressed as the inner product of the main field and the adjoint field solvers in the perturbation region, without the need for repeated forward simulations. It automatically identifies the geometric or material perturbation operators corresponding to each variable and calculates the sensitivity. Step 5.3: Package the gradient vectors of all design variables with the current performance metrics for subsequent optimization feedback of the geometric configuration; Step 5.4: For metamaterial microstructures, based on their local incident and transmitted acoustic field data, the parametric inversion method is used to calculate their true acoustic response at the operating frequency.

9. The method for modeling and optimizing an ultrasonic enhancement device as described in claim 8, characterized in that, Step 6 specifically includes: Step 6.1: Introduce a deep neural network model. The input is the normalized geometric parameters of the metamaterial microstructure, and the output is the predicted equivalent acoustic response. The deep neural network surrogate model mainly consists of an input layer, a DNN surrogate model body, an output layer, and a monitoring and correction layer. The DNN surrogate model body includes a feature extraction block, an activation layer, a pooling layer, and a fully connected layer connected in sequence. Step 6.2: The deep neural network model is pre-trained using 5000 sets of offline-generated full-wave simulation datasets before optimization begins. At the beginning of each optimization iteration, the model is called to quickly predict the macroscopic acoustic response of candidate metamaterial microstructures. Step 6.3: Each time the surrogate model compares the real acoustic response with the equivalent acoustic response, if the relative error of any output component exceeds 5%, a local correction mechanism is triggered by monitoring the correction layer: the input sample is added to the training set, fine-tuning is performed only on the weights of the last two layers of the DNN surrogate model, and the DNN is restricted from participating in the adjoint gradient calculation, serving only as an acceleration tool to provide the initial configuration and rapid evaluation; finally, the corrected DNN surrogate model is reconnected.

10. The method for modeling and optimizing an ultrasonic enhancement device as described in claim 7, characterized in that, Step 7 specifically includes: Step 7.1: Define the objective function, using the weighted difference between maximizing the energy in the focal region and minimizing the energy in the sidelobe region as the optimization guide, where the focal energy weight coefficient and the sidelobe suppression weight coefficient can be dynamically adjusted; Step 7.2: The optimizer integrates the gradient information provided in Step 5 with the equivalent acoustic response predicted by the model in Step 6 to update the design variables. Step 7.3: Immediately after each parameter update, perform three hard constraint checks on the new parameters:

1. Minimum feature size; 2. Overall integral; 3. Structural connectivity; If any constraint is violated, pull the parameters back into the feasible region through a projection operation; Step 7.4: Continuously monitor the improvement rate of the objective function. In each round, check the relative rate of change of the objective function to determine convergence. Output the final three-dimensional geometric model, material distribution map and process parameter package for use in additive manufacturing equipment driving.

Images