Parallel simulation method for guided wave focusing imaging of composite structure piezoelectric array

Abstract

The invention discloses a parallel simulation method for guided wave focusing imaging of a piezoelectric array of a composite structure, belongs to the technical field of nondestructive testing, and aims at solving the problems that the sound field control precision of a current multi-array-element piezoelectric structure in a composite material plate is insufficient, and a traditional finite element method is large in calculation time consumption, low in convergence speed and high in reliability in the guided wave propagation simulation process. In order to solve the problem that effective convergence cannot be realized under the excitation condition of a high-frequency multi-channel array and even under the excitation condition of a high-frequency multi-channel array, the invention provides an efficient guided wave focusing simulation strategy combining spectral element high-order precision and GPU parallel computing, and the method comprises the following steps: establishing a system control equation considering a piezoelectric-structure coupling effect on the basis of a delay superposition principle; energy focusing and sidelobe suppression of a guided wave signal in a specific direction are realized by adjusting spatial arrangement and excitation time sequence of a piezoelectric array sensor, Gauss-Labatto-Legendre integral nodes are introduced to construct a high-order spectral element unit, unit-level parallel solution is realized under a CUDA architecture, and dynamic simulation capability and overall calculation efficiency of focusing field distribution are remarkably improved.

Description

TECHNICAL FIELD

[0001] The present application relates to a parallel simulation method for composite structure piezoelectric array guided wave focusing imaging, belonging to the technical field of nondestructive testing. BACKGROUND

[0002] The spectral element method has the flexibility of finite element in dealing with complex boundaries and structures, and the rapid convergence characteristics of spectral method. Under the premise of achieving the same accuracy, it can adopt sparser unit division than traditional finite element. In the spectral element method, the global matrix is aggregated from the element matrix according to a specific structure, so most of the work can be done on the local element in the calculation process, and load balancing between elements is easy to achieve. Through the Compute Unified Device Architecture (CUDA), the calculation task can be distributed to multiple computing nodes, and each node independently calculates the allocated element matrix and stores it in its own memory. In this way, not only the calculation efficiency is improved, but also the memory demand is greatly reduced. The phased array technology is an ultrasonic detection technology that uses multiple piezoelectric sensors to form a piezoelectric array and controls the excitation timing of each sensor to realize the deflection and focusing of the sound beam. It can make the ultrasonic beam focus at different angles and different distances.

[0003] In some complex industrial detection scenarios, such as nondestructive testing of large structural parts, large-scale computing tasks and high-precision sound field focusing problems need to be handled simultaneously. The time-domain spectral element method can be used for parallel computing to quickly simulate the propagation and scattering of ultrasonic guided waves in structures, and combined with the phased array focusing method, high-precision sound field control and defect detection in specific areas can be realized. The combination of the two can fully leverage their respective advantages and improve detection efficiency and accuracy. This combination has broad application prospects in the field of nondestructive testing and can promote the research and application development in related fields. SUMMARY

[0004] The present application proposes a parallel simulation method for composite structure piezoelectric array guided wave focusing imaging, aiming to realize fast detection of large-area composite structures and improve the calculation efficiency of piezoelectric array ultrasonic guided wave focusing simulation using GPU parallel computing.

[0005] In order to achieve the above technical purpose, the technical scheme of the present application is as follows:

[0006] A parallel simulation method for composite structure piezoelectric array guided wave focusing imaging, characterized in that the method comprises:

[0007] Step 1: Use ABAQUS software to divide the composite material plate into 4-node quadrilateral element mesh, calculate the one-dimensional n-node Gauss-Lobatto-Legendre integral point coordinates, and in an n-order spectral element, the node coordinates are in the range of [-1, 1], and the normalized coordinates of the nodes can be represented as:

[0008] ξ1 = -1, ξ2 = r1, ξ3 = r2, ξ n-1 = r n-2 , ξ n = 1

[0009] where r1, r2,..., r n-2 represent the roots of (n-2) order Lobatto polynomial L n-2 (ξ), the coordinates of the nodes ξ i of a n order spectral element are obtained by the following formula:

[0010] (1-ξ 2 )L n-2 (ξ) = 0

[0011] where L n-2 (ξ) is Lobatto polynomial:

[0012]

[0013] where P n (ξ) is the first derivative of (n+1) order Legendre polynomial:

[0014]

[0015] The 4-node quadrilateral element is converted into n 2 node spectral element by Gauss-Lobatto-Legendre integral point coordinates, and the new node coordinates are obtained by the following formula:

[0016]

[0017] ξ i is the integral point coordinate, x is the coordinate of the midpoint of the element edge when calculating the node on the element edge, and l is the length of the element edge, x is the coordinate of the midpoint of the connecting line between the two nodes on the element edge when calculating the middle node in the element, and l is the distance between the two nodes;

[0018] Step two: define the parameters required for solving, including material properties, geometric dimensions of the model, and excitation signal s(t) and time step Δt:

[0019]

[0020] where f c is the center frequency of the excitation signal, t is the time, and n is the number of modulation periods;

[0021] Step three: according to Hamilton's principle, the system control equation considering piezoelectric-structure coupling is established:

[0022] σ = c E ε - d T E

[0023] U D = dε + g ε E

[0024] where U D is the electric displacement matrix, σ is the stress vector, ε is the strain vector, c E is the elastic constant matrix at constant electric field intensity (superscript E indicates measured at constant electric field intensity), d (superscript T indicates transpose) is the piezoelectric constant matrix, g ε is the dielectric constant tensor at constant stress (superscript ε indicates measured at constant strain), and E is the electric field intensity. The relationship between the electric field intensity and the electric potential in the electric field is:

[0025]

[0026] where φ represents the electric potential field, and the electric field-node potential matrix is defined as:

[0027]

[0028] where the subscript represents the physical quantity of the electric field, and S(ξ, η, ζ) is the shape function. According to the Hamilton principle, the system control equation of three-dimensional piezoelectric coupling can be obtained as:

[0029]

[0030] where M e is the mass matrix, is the stiffness matrix of the mechanical field, is the dielectric constant matrix, and are the piezoelectric coupling matrices, φ e is the electric potential, f e is the node external force vector, u e , are the node displacement and acceleration vectors, respectively, and Q is the externally applied charge; and are defined as:

[0031]

[0032] where the subscript u represents the physical quantity of the displacement field, w is the integral weight, is the electric field-node potential matrix, B u is the strain-displacement matrix (superscript T indicates transpose), and J eis the Jacobian matrix, which represents the mapping from the local coordinate system to the global coordinate system, for typical piezoelectric materials, The value of 8 is about 10 The value of -11 is about 10 s This huge difference in magnitude will cause a non-negligible error in calculation, in order to overcome this problem, the static condensation method is used to condense the term related to , which is expressed in the form of displacement field vector:

[0033]

[0034] In the formula, K s is the stiffness matrix caused by piezoelectric coupling, f A is the equivalent node force vector caused by the voltage applied to the upper surface of the piezoelectric sensor, which can be expressed as:

[0035]

[0036] When the piezoelectric sensor is used as a driver, the excitation voltage V is applied to the upper surface, and the potential of the lower surface is 0, then the unknown induced potential in the middle layer of the piezoelectric element can be expressed as:

[0037]

[0038] In the formula, the upper surface i represents the middle layer of the piezoelectric element, and n represents the upper surface of the piezoelectric element.

[0039] Step four: considering a one-dimensional linear wave phased array composed of N piezoelectric sensor elements, the array is numbered from left to right as 1 to N, assuming that there is a point P (|d|, θ) in the plane, d is the vector from the coordinate origin O to the point P, there is a sensor i in the column, the vector from the sensor to the origin is z i , and the vector of point P is d i , where θ is the angle between vector d and x axis, and has:

[0040]

[0041] d i = d - z i

[0042] In the formula, l is the arrangement distance of the piezoelectric sensor, let the unit vector of the vector from sensor i to point P be τ i , and τ i is:

[0043]

[0044] Then the wave number vector of sensor i is:

[0045]

[0046] When the guided wave excited by the sensor at the origin of the plate propagates to the focus point:

[0047]

[0048] Where d is the distance from the focus point to the origin of the coordinate axis, k is the wave number, A is the amplitude coefficient, e is the base of the natural logarithm, and ω is the angular frequency. When the guided wave excited by the i-th sensor propagates to the focus point, there is:

[0049]

[0050] Where d i is the distance from the focus point to the i-th sensor, and the synthesized beam of the guided wave received at the focus point generated by the N sensor excitation signals is:

[0051]

[0052] Where w i is the weight coefficient of each sensor, and when the amplitude of the sensor excitation signal is consistent, w i = 1:

[0053]

[0054] The latter part is the directivity function of the guided wave beam, which plays a key role in controlling the deflection and focusing of the beam. Therefore, for a one-dimensional linear guided wave phased array composed of N piezoelectric sensor elements, the final beam synthesis formula after applying the delay is:

[0055]

[0056] The delay time of the excitation signal of the i-th sensor is:

[0057]

[0058] Where θ is the target angle of the deflection and focusing of the synthesized beam, and v is the phase velocity corresponding to the propagation angle.

[0059] Step five: decompose the common nodes of the elements using the grid numbering algorithm. Considering four spectral elements with common nodes, if each spectral element has 36 nodes, the correspondence between the global node number and the local node number is:

[0060]

[0061] Where W is the global node number when the spectral elements are connected, L is the local node number after domain decomposition, and the superscript is the element number, indicating that the node belongs to which element.

[0062] Strain is calculated for the separated spectral elements Stress Internal force The mapping connecting global index W and local index L is created, and the index pair is used for parallel computation to assemble global force vector F W And displacement vector

[0063]

[0064] Wherein And The element node force matrix and the element node displacement matrix obtained by parallel computation, F W And U W The assembled global force matrix and global displacement matrix;

[0065] Finally, the displacement at time step t+Delta t is calculated according to the step-by-step time integration algorithm, wherein t represents the current time, and Delta t is the time step length;

[0066]

[0067] In the formula, U t-Δt , U t , U t+Δt Respectively represent the displacement at the last time, the displacement at the current time and the displacement to be solved at the next time, Delta t is the time step length, F t Is the excitation at the current time t, M is the global mass matrix, and C is the global damping matrix.

[0068] Step six: the full wave field data is visualized, and the wave field is approximated by using node displacement and shape function.

[0069] The application has the advantages and beneficial effects:

[0070] The application discloses a parallel simulation method for focusing wave imaging of a piezoelectric array of a composite structure, and establishes system control equations considering piezoelectric-structure coupling; based on a classical delay and superposition principle, subwaves generated by each sensor are superposed by changing arrangement intervals between the piezoelectric sensors and excitation delay time, so that spatial selectivity and amplification of a transmission signal in a required direction are realized, and suppression in other directions is realized, defects with different angles in a composite material plate are detected in a full-angle scanning mode, influence of a phase velocity difference between different propagation angles caused by an anisotropic material is considered, and multi-angle defect positioning of the composite material is realized. When a dynamic propagation process of the coupled wave is solved, a grid is divided into subblocks (domain decomposition), parallelization is realized through calculation of a unified device architecture (CUDA), unit-level parallel computation for solving a wave equation is realized on a GPU, and calculation efficiency is improved, so that the method has wide application prospects in ultrasonic wave phased array detection of a large composite material plate-shaped structure. BRIEF DESCRIPTION OF DRAWINGS

[0071] In order to more clearly illustrate the technical solutions in the embodiments of the present application or the prior art, the following will briefly introduce the drawings needed to be used in the embodiments or the prior art description. Obviously, the drawings in the following description are only some embodiments of the present application, and other drawings can be obtained by those skilled in the art without any creative effort on the basis of the drawings shown.

[0072] Figure 1 is a flowchart of the method of the present application;

[0073] Fig. 2 is a schematic diagram of a model grid;

[0074] Figure 3 is a 90-degree direction focusing wave field diagram of a [0\90\45\-45\45\-45\90\0]s laminated plate;

[0075] Figure 4 is a 90-degree direction focusing defect wave field diagram of a [0\90\45\-45\45\-45\90\0]s laminated plate;

[0076] Figure 5 is a 60-degree direction focusing wave field diagram of a [0\90\45\-45\45\-45\90\0]s laminated plate;

[0077] Figure 6 is a 60-degree direction focusing defect wave field diagram of a [0\90\45\-45\45\-45\90\0]s laminated plate;

[0078] Figure 7 is a 90-degree direction focusing directivity diagram;

[0079] Figure 8is a 60-degree directional wave focusing directivity diagram. DETAILED DESCRIPTION

[0080] In order to make the purposes, technical solutions and advantages of the embodiments of the present application clearer, the technical solutions in the embodiments of the present application will be described clearly and completely below with reference to the drawings in the embodiments of the present application. Obviously, the described embodiments are part of the embodiments of the present application, rather than all the embodiments of the present application. Based on the embodiments in the present application, all other embodiments obtained by a person of ordinary skill in the art without creative labor fall within the protection scope of the present application.

[0081] ABAQUS is a powerful engineering simulation finite element software, and the meshing function of ABAQUS is a crucial part of finite element analysis, which decomposes large structures into small and simple elements. By setting local mesh density, defining characteristic lines, local refinement areas and other ways, it ensures sufficient mesh density in important areas to improve solving accuracy and efficiency. Gauss-Lobatto-Legendre is a numerical integration method used to calculate definite integrals within a given interval. It combines Gauss-Legendre integration points and endpoints (interval boundaries) to provide higher accuracy in calculations. GPU (Graphics Processing Unit) is a processor specifically designed for processing graphics and images, widely used in parallel computing and other fields. Compared with CPU, the parallel processing capability of GPU makes it more efficient in processing large amounts of data, and can handle up to thousands of threads in parallel. Hamilton's principle is a principle based on the variational method, which describes the motion law of a physical system through the extremum condition of the action. The action is a scalar quantity, and for a piezoelectric-structure coupled system, the action is the integral of the Lagrangian over time. Through Hamilton's principle, the control equations of the system considering piezoelectric-structure coupling can be established. These equations include the motion equations of the structure, the motion equations of the piezoelectric material and the electric field equations, as well as the coupling conditions. The role of Hamilton's principle is to unify the kinetic energy, potential energy and coupling terms of the system through the extremum condition of the action, so as to obtain the complete dynamics description of the system.

[0082] Example: As shown in Figure 1 The present application provides a parallel simulation method for composite structure piezoelectric array wave focusing imaging:

[0083] 1) Step one of the method needs to calculate one-dimensional 6-node Gauss-Lobatto-Legendre integral point coordinates. In an n-order spectral element, the node coordinates are in the range of [-1, 1], and the normalized coordinates of the nodes can be expressed as:

[0084] ξ1=-1,ξ2=r1,ξ3=r2,ξn-1 = r n-2 , ξ n = 1

[0085] where r1, r2,..., r n-2 denote the roots of the (n-2)th order Lobatto polynomial L n-2 (ξ), the coordinates of the nodes ξ i of a one n-th order spectral element are obtained by the following formula:

[0086] (1-ξ 2 )L n-2 (ξ) = 0

[0087] where L n-2 (ξ) is the Lobatto polynomial:

[0088]

[0089] where P n (ξ) is the set of first order derivatives of the (n+1)th order Legendre polynomial:

[0090]

[0091] The 4-node quadrilateral element is converted into an n 2 node spectral element by the Gauss-Lobatto-Legendre integration point coordinates, and the new node coordinates are obtained by the following formula:

[0092] X i = x + (l / 2)*ξ i

[0093] ξ i is the integration point coordinate, x is the coordinate of the midpoint of the element edge when calculating the node on the element edge, and l is the length of the element edge. When calculating the intermediate node of the element, x is the coordinate of the midpoint of the connecting line between the two nodes on the element edge, and l is the distance of the connecting line between the two nodes.

[0094] 2) The material used in this example is a 16-layer [0\90\45\-45\45\-45\90\0] s symmetric structure carbon fiber reinforced polymer (CFRP) composite plate, the single layer thickness is 0.15mm, the total thickness is 2.4mm, the geometric size is 1000mm×1000mm×2.4mm, as shown in Figure 2, a one-dimensional piezoelectric sensor array is arranged at the center of the plate, and two delamination defects with an angle of 60 degrees and 90 degrees and a diameter of 10mm are arranged at a distance of 200mm from the center of the plate between the 4th layer and the 5th layer of the composite material plate;

[0095] The overall stiffness matrix and mass matrix are calculated by the defined material parameters. The piezoelectric layer is also regarded as a composite layer to calculate the overall stiffness of the model. The formula for calculating the mass matrix and stiffness matrix is:

[0096]

[0097] H (c) (ξ,η) is the coefficient matrix related to the inertia term:

[0098]

[0099] where m is the number of composite layers, ρ r is the mass density of the rth layer, l r is the distance from the neutral surface of the composite plate to the upper surface of the rth layer, l r-1 is the distance from the neutral surface of the composite plate to the lower surface of the rth layer.

[0100] D (c) in the formula is the elastic matrix:

[0101]

[0102] According to the stacking sequence, the off-axis stiffness matrix of the rotating layer is calculated from the on-axis stiffness matrix

[0103]

[0104] where m = cosθ, n = sinθ, and θ is the layup angle. After obtaining the off-axis stiffness matrix of each rotating layer, the integral of the layer thickness in the elastic matrix is analyzed and calculated to obtain the stiffness coefficients:

[0105]

[0106] where n is the number of layers of the CFRP plate, ρ is the density, l is the total thickness of the shell, z r is the mid-surface coordinate of the rth layer of the layup. a ij is the stiffness coefficient related to the internal force and the neutral surface strain, which is independent of the material layup sequence, b ij , d ij are related to the material layup sequence, where b ij represents the coupling relationship between torsion, bending, and tension, and d ij is the stiffness coefficient related to the internal force moment and curvature and torsion curvature, which is the bending stiffness coefficient; h ij is the shear stiffness coefficient.

[0107] 3) According to Hamilton's principle, the system control equation considering piezoelectric-structure coupling is established:​

[0108] σ = c E ε - d T E

[0109] U D = dε + g ε E

[0110] where U D is the electric displacement matrix, p is the mass density, p e is the free charge density, f is the body force, σ is the stress vector, ε is the strain vector, c E is the elastic constant matrix at constant electric field intensity (superscript E indicates measured at constant electric field intensity), d is the piezoelectric constant matrix, g ε is the dielectric constant tensor at constant stress (superscript ε indicates measured at constant strain), E is the electric field intensity. In the electric field, the relationship between the electric field intensity and the electric potential is:

[0111]

[0112] where φ represents the electric potential field, the electric field-node potential matrix is defined as:

[0113]

[0114] where the subscript represents the physical quantity of the electric field, S(ξ, η, ζ) is the shape function, and according to the Hamilton principle, the system control equation of three-dimensional piezoelectric coupling can be obtained as:

[0115]

[0116] where M e is the mass matrix, is the stiffness matrix of the mechanical field, is the dielectric constant matrix, and are the piezoelectric coupling matrices, φ e is the electric potential, f e is the node external force vector, u e , are the node displacement and acceleration vectors, respectively, and Q is the externally applied charge; and are defined as:

[0117]

[0118] where the subscript u represents the physical quantity of the displacement field, w is the integral weight, is the electric field-node potential matrix, B uis the strain-displacement matrix (superscript T denotes transpose), J e is the Jacobian matrix, representing the mapping of the element from the local to the global coordinate system, for a typical piezoelectric material, is of the order of 10 8 while is of the order of 10 -11 This huge difference in magnitude can cause non-negligible errors in the calculation, in order to overcome this problem, the static condensation method is used to condense the terms related to in the form of displacement field vector:

[0119]

[0120] where K s is the stiffness matrix caused by piezoelectric coupling, f A is the equivalent node force vector caused by the voltage applied on the upper surface of the piezoelectric sensor, which can be expressed as:

[0121]

[0122] When the piezoelectric sensor is used as a driver, the excitation voltage V is applied on the upper surface, and the potential on the lower surface is 0, then the unknown induced potential in the middle layer of the piezoelectric element can be expressed as:

[0123]

[0124] where the upper i represents the middle layer of the piezoelectric element, and n represents the upper surface of the piezoelectric element;

[0125] 4) Consider a one-dimensional linear wave phased array composed of N piezoelectric sensor elements, the array is numbered from left to right as 1 to N, assuming that there is a point P (|d|, θ) in the plane, d is the vector from the coordinate origin O to the point P, there is a sensor i in the column, the vector from the sensor to the origin is z i , the vector of point P is d i , where θ is the angle between vector d and x-axis, and has:

[0126]

[0127] d i = d - z i

[0128] where l is the arrangement distance of the piezoelectric sensor, let the unit vector of the vector from sensor i to point P be τ i , τ i is:

[0129]

[0130] The wave number vector of sensor i is:

[0131]

[0132] When the guided wave excited by the sensor at the origin of the plate propagates to the focus point, we have:

[0133]

[0134] where d is the distance from the focus point to the origin of the coordinate axis, k is the wave number, A is the amplitude coefficient, e is the base of the natural logarithm, and ω is the angular frequency. When the guided wave excited by the i-th sensor propagates to the focus point, we have:

[0135]

[0136] where d i is the distance from the focus point to the i-th sensor, and the synthesized beam at the focus point is:

[0137]

[0138] where w i is the weight coefficient of each sensor, and w i = 1 when the amplitudes of the sensor excitation signals are consistent:

[0139]

[0140] The latter part is the directivity function of the guided wave beam, which plays a key role in controlling the deflection and focusing of the beam. Therefore, for a one-dimensional linear guided wave phased array composed of N piezoelectric sensor elements, the final beam synthesis formula after applying the delay is:

[0141]

[0142] As shown in Figure 7 and 8 , the synthesis effect of wave velocity in 60-degree and 90-degree directions is calculated and compared with the numerical simulation results. The delay time of the excitation signal of the i-th sensor is:

[0143]

[0144] where θ is the target angle of the deflection and focusing of the synthesized beam, and v is the phase velocity corresponding to the propagation angle.

[0145] 5) The common nodes of the elements are decomposed using the grid numbering algorithm. Considering four spectral elements with common nodes, if each spectral element has 36 nodes, the correspondence between the global node number and the local node number is:

[0146] ​

[0147] where W is the global node number of the spectral element when it is connected, L is the local node number after domain decomposition, the superscript is the element number, and it indicates that the node belongs to which element;

[0148] Calculate the strain for the separated spectral element Stress Internal force A mapping connecting the global index W and the local index L is created, and the index pair is used to assemble the global force vector F after parallel computing W and displacement vector

[0149]

[0150] where and are the element node force matrix and the element node displacement matrix obtained by parallel computing, F W and U W are the assembled global force matrix and the global displacement matrix;

[0151] 6) Finally, the displacement at time step t+Δt is calculated according to the step-by-step time integration algorithm, where t represents the current time, and Δt is the time step;

[0152]

[0153] where U t-Δt , U t , U t+Δt represent the displacement at the previous time, the current time, and the next time to be solved, respectively, Δt is the time step, F t is the excitation at the current time t, M is the global mass matrix, and C is the global damping matrix.

[0154] 7) After obtaining all the node displacement data, a uniform grid point interpolation technique is used to visualize the full wave field data, as shown in Figure 3 , 4 , 5 and 6, the 60-degree and 90-degree directional layered defect and defect focusing wave field diagrams are calculated.

[0155] First, create uniform data points (xp, yp), find the spectral element containing the data points, find the spectral node (x0, y0) closest to the data point (xp, yp) in the spectral element as the initial point, and perform parallel computing for all data points in the regular grid:

[0156]

[0157] where k=0 in the first iteration. Then the shape functions and Jacobian matrix are calculated at the point (ξ1,η1), and the global coordinates (x1,y1) corresponding to the interpolation point of the local coordinates (ξ1,η1) are calculated using the following formula:

[0158]

[0159] After 3 iterations, the new interpolation local coordinates (ξ p ,η p ) are obtained, and the spectral element shape functions are calculated at the interpolation local coordinates (ξ p ,η p ):

[0160] S (c) =(Q) T *k

[0161] where k(i,j)=ξ i-1 *η j-1 (i=1,...,6,j=1,...,6), and Q is the coefficient of the Vandermonde matrix;

[0162] Finally, the displacement data of all nodes of the spectral element containing the data point (xp,yp) are extracted, and the wave field is approximated by the node displacement and the shape function:

[0163]

[0164] The purpose of the present application is to solve the problems of low efficiency of multi-pressure structure coupling calculation and high-precision sound field control and arbitrary position defect detection of large composite plate structure, and to realize rapid detection of large-area composite structure and improve the efficiency of piezoelectric array ultrasonic guided wave focusing simulation calculation using GPU parallel calculation. The method of the present application combines high-order time-domain spectral element method parallel calculation with piezoelectric array ultrasonic guided wave deflection focusing technology to realize large composite material defect detection, improve detection efficiency and accuracy.

[0165] The above is a further detailed description of the present application in combination with specific preferred embodiments, and the specific implementation of the present application should not be limited to these descriptions. For ordinary skilled persons in the technical field to which the present application belongs, some simple deductions or substitutions can be made without departing from the concept of the present application, and all of them should be regarded as falling within the protection scope of the present application.

Claims

1. A parallelized simulation method for piezoelectric array guided wave focusing imaging of composite structures, characterized in that, The method comprises the following steps: Step 1: using ABAQUS software to divide the composite plate into a quadrilateral element grid, and converting the quadrilateral element model into a spectral element model through Gauss-Lobatto-Legendre integral point coordinates; Step 2: defining parameters required for solving, including material properties, geometric dimensions of the model, and excitation signal S(t) and time step Δt; Step 3: according to the Hamilton principle, establishing a system control equation considering piezoelectric-structure coupling, and solving to obtain equivalent node forces caused by voltage; Step 4: taking the center of the plate as the origin, determining the distance of the focal point to the origin and the deflection angle, and based on the classic delay and superposition principle, establishing a general formula of anisotropic composite plate phased array beam forming to calculate the arrangement spacing and excitation delay time between piezoelectric sensors; Step 5: splitting the spectral element model, and using a parallel computing method to solve the wave equation at the element level; Step 6: visualizing the full wave field data, and using node displacement and shape function to approximate the wave field.

2. A parallelized simulation method for piezoelectric array guided wave focusing imaging of composite structures according to claim 1, characterized in that, In step three, the constitutive equation considering piezoelectric-structure coupling is: σ = c E ε - d T E U D = dε + g ε E Among them U D It is the electric displacement matrix, ρ is the mass density, ρ e Let f be the free charge density, f be the body force, σ be the stress vector, ε be the strain vector, and c be the free charge density. E The elastic constant matrix under constant electric field strength (superscript E indicates measurement at a fixed electric field strength), d is the piezoelectric constant matrix, g ε Let ε be the dielectric constant tensor under constant stress (the superscript ε indicates measurement at constant strain), and E be the electric field strength. Within the electric field, the relationship between electric field strength and electric potential is: where φ represents the electric potential field, the electric field-nodal potential matrix is defined as: where the subscript The physical quantity representing the electric field is S(ξ,η,ζ), which is a shape function. According to Hamilton's principle, the control equation of the three-dimensional piezoelectric coupling system is obtained as where M e is the mass matrix, is the stiffness matrix of the mechanical field, is the permittivity matrix, and is the piezoelectric coupling matrix, φ e is the electric potential, f e is the nodal force vector, u e , are the nodal displacement, acceleration vectors, respectively, and Q is the externally applied charge; and are defined as: where subscript u denotes the physical quantity of the displacement field, w is the integral weight, is the electric field-node potential matrix, B u is the strain-displacement matrix (superscript T denotes transposition), J e is the Jacobi matrix, which represents the mapping of the element from the local coordinate system to the global coordinate system, for typical piezoelectric materials, The value of is about 10 8 and the value of is about 10 The value of is about 10 -11 This huge difference in magnitude will cause a non-negligible error in calculation, in order to overcome this problem, the static condensation method is used to condense the term related to in the form of displacement field vector: where K s is the stiffness matrix due to piezoelectric coupling, f A is the equivalent nodal force vector due to the voltage applied on the top surface of the piezoelectric sensor, which can be expressed as: When the piezoelectric sensor is used as a driver, the excitation voltage V is applied to the upper surface, and the potential of the lower surface is 0, then the unknown induced potential of the middle layer of the piezoelectric element can be expressed as: In the formula, the upper surface i represents the middle layer of the piezoelectric element, and n represents the upper surface of the piezoelectric element.

3. The parallelized simulation method for piezoelectric array guided wave focusing imaging of composite structures according to claim 1, wherein, In step four, considering a one-dimensional linear guided wave phased array consisting of N piezoelectric sensor elements, the array is numbered from left to right as 1 to N, and assuming that there is a point P (|d|, θ) in the plane, d is the vector from the coordinate origin O to the point P, and there is a sensor i in the column, the vector from the sensor to the origin is z i , the vector of the point P is d i , where θ is the angle between the vector d and the x-axis, and there is: d i = d - z i where l is the piezoelectric sensor array pitch, and let τ be the unit vector of the vector from sensor i to point P i , τ i is given by: Then the wave number vector of the sensor i is: When the sensor at the origin of the plate excites the guided wave to the focal point: In the formula, d is the distance from the focal position to the coordinate axis origin, k is the wave number, A is the amplitude coefficient, e is the base of natural logarithm, ω is the angular frequency, then when the i-th sensor excites the guided wave to the focal point: where d i is the distance from the focus position to the ith sensor, and the received wave at the focus point is the sum of the N guided waves: where w i is the weight coefficient of each sensor, and w i = 1 when the sensor excitation signal amplitudes are consistent. In the formula, the latter part is the directivity function of the guided wave beam, which plays a key role in controlling the deflection and focusing of the beam, therefore, for a one-dimensional linear guided wave phased array composed of N piezoelectric sensor elements, the final beam synthesis formula after applying the delay is: The excitation signal delay time of the i-th sensor is: In the formula, θ is the target angle of the deflection and focusing of the synthesized beam, and v is the phase velocity corresponding to the propagation angle.

4. The parallelized simulation method for piezoelectric array guided wave focusing imaging of composite structures according to claim 1, wherein, In step five, four spectral elements with a common node are considered, if each spectral element has 36 nodes, the corresponding relationship between the global node number and the local node number is: Wherein W is the global node number when the spectral elements are connected, L is the local node number after domain decomposition, the superscript is the element number, and it is indicated that the node belongs to which element; Strain is computed for the individual spectral elements Stress Internal forces A map is created that connects the global indices W and the local indices L, which are used to assemble the global force vector F after parallel computation W And displacement vectors wherein and are the element node force matrix and element node displacement matrix obtained by parallel computation, F W and U W are the assembled global force matrix and global displacement matrix; Finally, the displacement at time step t+Δt is calculated according to the step-by-step time integration algorithm, wherein t represents the current time, and Δt is the time step. In the formula, U t-Δt , U t , U t+Δt respectively represent the displacement at the previous time, the displacement at the current time and the displacement at the next time to be solved, Δt is the time step, F t is the excitation at the current time t, M is the global mass matrix, and C is the global damping matrix.

Images