A dynamic fracture phase field calculation method for a hybrid unit beam structure

Abstract

The application belongs to the technical field of computational fracture mechanics, and particularly relates to a dynamic fracture phase field calculation method for a hybrid unit beam structure. The method comprises the following steps: determining material parameters required for dynamic fracture phase field calculation of the hybrid unit beam structure; establishing a dynamic fracture phase field mathematical model for the hybrid unit beam structure according to the Lagrange action principle, wherein the mathematical model comprises a displacement field control equation and a phase field evolution control equation; performing numerical discretization on the mathematical model to establish a discrete model of the hybrid unit beam structure; performing numerical solution on the numerically-discretized mathematical model by using the Newmark method combined with the staggered algorithm until the displacement field and the phase field both satisfy a convergence condition, and obtaining dynamic fracture phase field evolution results of the hybrid unit beam structure. The application solves the problem that the existing fracture phase field model cannot accurately simulate crack evolution on the thickness of a beam structure under dynamic load, and significantly reduces the calculation cost of simulating dynamic fracture of the beam structure.

Description

Technical Field

[0001] This invention belongs to the field of computational fracture mechanics, and particularly relates to a method for calculating the dynamic fracture phase field of a hybrid unit beam structure. Background Technology

[0002] Accurate and efficient prediction of failure modes of beam structures under dynamic loads is of great significance for the safety design and reliability assessment of large-scale engineering structures. Existing methods for predicting structural fracture often face limitations such as complex equipment and high costs in physical experiments; therefore, numerical simulation has gradually become the primary means of predicting structural fracture behavior. Existing numerical methods for simulating fracture can be broadly classified into discrete and continuous methods. Discrete methods struggle to effectively describe the complex behaviors of cracks, such as branching and merging, during their evolution. As a type of continuous method, the fracture phase-field method has received widespread attention in recent years. This method introduces a phase field to represent cracks in a diffuse manner, naturally simulating crack initiation, propagation, branching, and merging. However, existing dynamic fracture phase-field methods typically use two-dimensional or three-dimensional continuum elements to discretize the structure, resulting in a large number of computational degrees of freedom and high computational costs. Furthermore, for fracture analysis of beam structures, traditional phase-field models based on one-dimensional beam elements often assume that the phase field evolution is consistent along the beam thickness direction. This assumption is difficult to accurately reflect the actual fracture behavior of beam structures under dynamic loading conditions. To overcome the above shortcomings, this invention proposes a method for calculating the dynamic fracture phase field of a hybrid element beam structure. The displacement field of the beam structure and the phase field characterizing fracture are discretized using beam elements and quadrilateral elements, respectively. This method can accurately describe the fracture evolution behavior of the beam structure in the thickness direction under dynamic loads and effectively reduce the computational cost. Summary of the Invention

[0003] To overcome the shortcomings of the prior art, this invention provides a method for calculating the dynamic fracture phase field of a hybrid unit beam structure, comprising the following steps:

[0004] Step 1: Determine the material parameters required for dynamic fracture phase field calculation of the hybrid unit beam structure;

[0005] Step 2: Establish a dynamic fracture phase field mathematical model for the hybrid unit beam structure based on the Lagrange action principle. The mathematical model includes displacement field control equations and phase field evolution control equations.

[0006] Step 3: Numerically discretize the mathematical model to establish a discrete model of the hybrid element beam structure: the displacement field is discretized using beam element mesh, and the phase field is discretized using quadrilateral element mesh;

[0007] The discrete model is as follows:

[0008] ;

[0009] in, This is the quality matrix; Here is the displacement field stiffness matrix; For acceleration; For displacement field; This is the vector of external forces; Here is the phase field stiffness matrix; For phase field; The phase field driving force vector;

[0010] Step 4: Using the Newmark method combined with the staggered algorithm, solve the displacement field and phase field in the discrete model until both the displacement field and phase field satisfy the convergence condition, and obtain the dynamic fracture phase field evolution results of the hybrid unit beam structure.

[0011] The Newmark method is specifically described as follows: In the beam structure at the... The displacement, velocity, and acceleration at time step are determined by the first... The displacement, velocity, and acceleration of the beam structure at each time step are obtained as follows:

[0012] ;

[0013] in, , and The first Displacement, velocity, and acceleration of the beam structure at each time step; , and The first Displacement, velocity, and acceleration of the beam structure at each time step; For time steps; For time step; and Set the parameters used in the Newmark method to ensure the unconditional stability of the numerical integration process. , .

[0014] The beneficial effects of this invention are as follows: This invention realizes the calculation of fracture evolution of beam structures under dynamic load, overcomes the problem that existing methods are difficult to accurately describe the thickness direction crack evolution of beam structures under dynamic load, and effectively reduces the calculation scale and calculation cost in dynamic fracture phase field analysis; in addition, this invention is not only applicable to fracture analysis of simple beam structures, but also to fracture simulation of complex beam structure systems such as lattice structures, and has good applicability and promotion value. Attached Figure Description

[0015] The accompanying drawings, which are included to provide a further understanding of the invention and form part of this application, illustrate the invention and are used to explain it, but do not constitute a limitation on the scope of protection of the invention.

[0016] Figure 1 This is a flowchart of the dynamic fracture phase field calculation method for hybrid unit beam structures according to the present invention;

[0017] Figure 2 This is a schematic diagram of the boundary conditions for the cantilever beam of the present invention;

[0018] Figure 3 This is a schematic diagram of the specific mesh structure of the discrete model of the hybrid unit beam structure of the present invention;

[0019] Figure 4 This is a discrete schematic diagram of the hybrid unit beam structure of the present invention;

[0020] Figure 5 This is a comparison diagram of the calculation results of the hybrid unit beam structure model of the present invention and the comparative two-dimensional fracture phase field model;

[0021] Figure 6 This is a comparison diagram of the energy response of the hybrid unit beam structure model of the present invention and the comparative two-dimensional fracture phase field model;

[0022] Figure 7 This is a comparison chart of the calculation time of the hybrid unit beam structure model of the present invention and the comparative two-dimensional fracture phase field model. Detailed Implementation

[0023] The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention. To make the above-mentioned objects, features and advantages of the present invention more apparent and understandable, the present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments.

[0024] like Figure 1 As shown, the method for calculating the dynamic fracture phase field of a hybrid unit beam structure according to the present invention includes the following steps:

[0025] Step 1: Determine the material parameters required for dynamic fracture phase field calculation of the hybrid unit beam structure. This embodiment uses a cantilever beam structure subjected to dynamic impact load as an example for illustration. Figure 2 As shown, the left end of the beam structure is fixed, forming a cantilever boundary condition; the right end, which is free, is subjected to a downward dynamic impact with an impact velocity of 16.5 m / s. The geometric dimensions and material properties used in the calculation are as follows:

[0026] Young's modulus Mass density Critical energy release rate ; Length dimension parameters ; beam width ; length of beam ; Beam thickness .

[0027] Step 2: Establish a dynamic fracture phase field mathematical model of the beam structure based on the Lagrange action principle. The mathematical model includes displacement field control equations and phase field evolution control equations.

[0028] Further details are as follows: The Lagrange action mentioned in step two is defined as the total kinetic energy of the system minus the total potential energy. The total potential energy includes strain energy, surface energy, and work done by external forces. The specific expression of the Lagrange action is as follows:

[0029] (1)

[0030] (2)

[0031] (3)

[0032] (4)

[0033] (5)

[0034] in, Let Lagrange's action function be used. The kinetic energy of the system; The strain energy of the system; The surface energy of the system; The work done by external forces is represented by u, which represents the displacement field. For speed; For phase field; Mass density; Parameters to ensure stability during numerical calculations; The tensile strain energy density; The compressive strain energy density; The critical energy release rate; For length dimensions; The phase field gradient; This indicates the distributed load applied along the axial direction of the beam structure; This indicates the distributed load applied laterally to the beam structure; and For the axial and lateral displacements of points on the neutral axis of the beam; For the computation of the domain volume integral; For the computational domain; It is a volume fractional element. In this embodiment, ; ; ; .

[0035] According to Hamilton's principle, for a conservative or non-conservative dynamical system, its true trajectory over any time interval should result in zero first variation of the system's action functional, i.e. Substituting formula (1) into the above equation, we get:

[0036] (6)

[0037] in, For variational operators; It is a time variable; Integral operators for time variables; It is the Lagrange action function; and These are the start and end times of the dynamic problem in the variational principle, respectively. Based on the Lagrange variational principle, it is necessary to take first-order variations for the displacement field and the phase field, respectively. When taking first-order variations for the displacement field, the kinetic energy term, strain energy term, and external force work term are respectively varied with respect to the displacement field. When taking first-order variations for the phase field, the strain energy term and surface energy term are respectively varied with respect to the phase field.

[0038] For a displacement field u, the variational form of its kinetic energy term is:

[0039] (7)

[0040] in, Let be the imaginary velocity. Integrating the above equation by parts yields:

[0041] (8)

[0042] in, This is a virtual displacement; Virtual velocity, For acceleration. According to Hamilton's principle, the variation at the time endpoint is zero, so the first term on the right-hand side of equation (8) is zero, and the final result of the kinetic energy variation is:

[0043] (9)

[0044] For a displacement field u, the variational form of its strain energy term is:

[0045] (10)

[0046] Since strain energy is a function of strain, according to the chain rule... Formula (10) can be written as:

[0047] (11)

[0048] in, For strain tensor; This is virtual strain. Based on the stress tensor... Definition, Formula (11) can be written as:

[0049] (12)

[0050] in, and These are tensile stress and compressive stress, respectively. Since the displacement field of this invention is an Euler-Bernoulli beam model considering axial deformation, the strain tensor... It has only one non-zero component, which is:

[0051] (13)

[0052] in, Let be the axial strain at any point within the beam. Its expression is:

[0053] (14)

[0054] in, and These represent the axial and lateral displacements of the beam's neutral axis, respectively. and Let be the axial and lateral coordinates of any point within the beam. In the Euler-Bernoulli beam model considering axial deformation, the stress tensor is... It also has only one non-zero component, written as:

[0055] (15)

[0056] in, Let be the axial stress at any point within the beam. The relationship between axial stress and axial strain is:

[0057] (16)

[0058] in, In this embodiment, the modulus is Young's modulus. According to formulas (13) and (15), the double dot product of the stress tensor and strain tensor in formula (11) can be simplified to:

[0059] = (17)

[0060] = (18)

[0061] in, and The axial direction of the beam is respectively The tensile stress component and compressive stress component in the direction; substituting formulas (17) and (18) into formula (12), formula (12) can be simplified to:

[0062] (19)

[0063] According to the stress expression of tension-compression decomposition, we have:

[0064] (20)

[0065] Substituting equation (20) into equation (19), the final form of the strain energy variation is:

[0066] (twenty one)

[0067] For the displacement field u, the variational form of its external force work term is:

[0068] (twenty two)

[0069] in, and These represent the axial and lateral virtual displacements of points on the neutral axis of the beam.

[0070] Phase field The variational form of its surface energy term is:

[0071] (twenty three)

[0072] in, It is a virtual phase field; This is the gradient of the virtual phase field.

[0073] Phase field The variational form of its strain energy term is:

[0074] (twenty four)

[0075] Substituting the kinetic energy variational formula (9), the displacement field strain energy variational formula (21), the external force work variational formula (22), the phase field surface energy variational formula (23), and the phase field strain energy variational formula (24) into Hamilton's principle formula (6), we obtain:

[0076] (25)

[0077] The variational factor for the displacement field in the above equation ( , , , The terms and the phase field variation ( , Grouping and rearranging the terms, and then inverting the full expression sign, we get:

[0078] (26)

[0079] According to the basic principles of variational methods, due to virtual displacement With virtual phase field They are independent and arbitrary, so that the above expression is valid in any time interval. For all conditions to hold true, the integral terms within the curly braces of equation (26) regarding both must be identically zero. Thus, the single equation (26) above is decoupled into two independent control equations, namely the equivalent weak integral forms of the displacement field control equation and the phase field evolution control equation:

[0080] (27)

[0081] To ensure the irreversibility of phase field evolution, a history variable function H is introduced, whose physical meaning is the maximum tensile strain energy density experienced by the material throughout the entire loading history, specifically expressed as:

[0082] (28)

[0083] in, The position vector is denoted by 'max'; 'max' is the maximum value operator. Time range At any time between.

[0084] The tensile strain energy density in formula (27) Replacing the historical variable function H in formula (28), the final forms of the displacement field control equation and the phase field evolution control equation are as follows:

[0085] (29)

[0086] The above governing equations represent the general form of the dynamic fracture phase-field model for the hybrid unit beam structure of this invention. In this embodiment, since the cantilever beam structure is loaded using end dynamic impact loading and no distributed load is applied to the beam, the following formula is used in the subsequent numerical implementation process: The displacement field control equation and the phase field evolution control equation (29) are simplified as follows:

[0087] (30)

[0088] Step 3: Numerically discretize the mathematical model to establish a hybrid element beam structure discretization model. Specifically, this includes: discretizing the displacement field and phase field using beam elements and quadrilateral elements respectively, based on the finite element method. For the displacement field, beam element meshes are generated along the beam structure's axial direction. The defined beam element node degrees of freedom include axial displacement, lateral displacement, and cross-sectional rotation, used to calculate the macroscopic deformation of the structure. For the phase field, two-dimensional planar four-node quadrilateral elements are used for discretization. In the thickness direction of each beam element's cross-section, the beam thickness region is divided into an upper region and a lower region, with the beam's neutral axis as the boundary, and upper and lower quadrilateral elements are respectively deployed. Both the upper and lower quadrilateral elements are configured to have only phase field degrees of freedom, excluding displacement degrees of freedom, thus refining the crack description in the thickness direction. In terms of coupling, the upper and lower quadrilateral elements are connected to the corresponding beam elements through shared Gaussian integration points. The specific coupling mechanism is as follows: The beam element calculates the strain energy density of the upper and lower regions at the integration points. The strain energy density of the upper region is transferred to the corresponding integration point of the upper quadrilateral element as the driving force for calculating the phase field of the upper region; the strain energy density of the lower region is transferred to the corresponding integration point of the lower quadrilateral element as the driving force for calculating the phase field of the lower region. The phase field calculated by the upper quadrilateral element is fed back to the beam element to reduce the material stiffness of the beam element at the integration point of the upper region; the phase field calculated by the lower quadrilateral element is fed back to the beam element to reduce the material stiffness of the beam element at the integration point of the lower region, thereby achieving data exchange between the displacement field and the phase field.

[0089] The specific mesh structure and nodal degrees of freedom definition of the discrete model of the hybrid element beam structure are as follows: Figure 3 As shown. Figure 3 The left side shows an overall schematic diagram of the hybrid unit mesh. The thick solid line in the middle represents the beam unit mesh arranged at the beam's neutral axis. Above and below this beam unit mesh, upper quadrilateral unit meshes and lower quadrilateral unit meshes are regularly arranged to cover the thickness area of ​​the beam. Figure 3 The right side shows the node parameter definitions for a single element; the lower part shows a beam element, which contains two nodes, left and right. These represent the axial displacements at the two end nodes of the beam element, respectively. Lateral displacement and the degree of freedom of the cross section rotation angle The upper part shows a quadrilateral element, which contains four nodes. These represent the phase field degrees of freedom at the four nodes, used to characterize the material fracture state in that region.

[0090] In this embodiment, the cantilever beam structure is discretized using a hybrid element method, such as... Figure 4As shown. In the axial direction of the beam, the neutral axis displacement field of the cantilever beam structure is discretized into 200 beam elements. Figure 4 The bold black line segment in the middle shows 200 interconnected beam elements. Adjacent beam elements are connected by nodes. These beam elements are used to describe the axial displacement, lateral displacement, and section rotation of the structure. Figure 4 The upper and lower quadrilateral element meshes are regularly arranged on the upper and lower sides of the central bold black line segment, respectively. Along the beam thickness direction, two quadrilateral element meshes are arranged corresponding to each beam element to discretize the phase field, thus forming a total of 400 quadrilateral elements in the entire structure. These quadrilateral elements maintain the same projection range as their corresponding beam elements in the beam axial direction and are located on the upper and lower sides of the neutral axis in the thickness direction.

[0091] To clarify the geometrical mapping relationship between beam elements and quadrilateral elements in the hybrid element discrete model, this invention adopts the following geometrical mapping relationship:

[0092] Let the i-th beam element be... Consisting of two adjacent beam nodes and The beam element corresponds to an upper quadrilateral element in the beam thickness direction. and a lower quadrilateral unit The corresponding quadrilateral element node numbering rules are defined as follows:

[0093] (31)

[0094] in, The upper quadrilateral element is the beam element, with its four nodes located at the top left, middle left, middle right, and top right positions respectively. The lower quadrilateral element corresponding to the beam element has its four nodes located at the left-middle, left-bottom, right-bottom, and right-middle positions, respectively. For each beam element, its corresponding upper and lower quadrilateral elements contain a total of eight nodes. Since the bottom edge of the upper quadrilateral element coincides with the top edge of the lower quadrilateral element, they share the node numbered on the neutral axis. and There are two nodes, so each beam element actually corresponds to six quadrilateral element nodes.

[0095] Let the i-th beam node be... Located on the neutral axis, its coordinates are This is then mapped to three quadrilateral element nodes along the beam thickness direction:

[0096] (32)

[0097] in, Let be the horizontal coordinate of the i-th beam node; Let be the vertical coordinate of the i-th beam node; This is the top edge node of the upper quadrilateral element; This is the bottom edge node of the lower quadrilateral unit; It is a common node located at the neutral axis, which simultaneously serves as the bottom edge node of the upper quadrilateral element and the top edge node of the lower quadrilateral element. The thickness of the beam.

[0098] In this embodiment, the first beam element is used. For example, the two beam nodes of this beam element are numbered as follows:

[0099] (33)

[0100] The corresponding coordinates of the beam neutral axis nodes are as follows:

[0101] (34)

[0102] Substituting the beam node numbers of the beam element in formula (33) into formula (31), the node numbers of the upper quadrilateral element are:

[0103] (35)

[0104] The node numbers of the lower quadrilateral elements are:

[0105] (36)

[0106] Substituting the coordinates of the beam's neutral axis nodes (34) into formula (32), the beam nodes... and The coordinates of the quadrilateral element nodes obtained by mapping along the thickness direction are as follows:

[0107] (37)

[0108] (38)

[0109] Based on the above calculation results, by substituting the calculated node coordinate formulas (37) and (38) into the element topology connections determined in formulas (35) and (36), the final geometric definitions of the upper and lower quadrilateral elements corresponding to the first beam element can be obtained:

[0110] Upper quadrilateral unit The physical coordinates of the four nodes (in counter-clockwise order: top left, middle left, middle right, top right) are:

[0111] (39)

[0112] Lower quadrilateral unit The physical coordinates of the four nodes (in counter-clockwise order: left-middle, left-bottom, right-bottom, right-middle) are:

[0113] (40)

[0114] Similarly, the remaining beam elements and their corresponding quadrilateral elements are generated according to the above mapping rules, thereby completing the hybrid element mesh discretization of the entire cantilever beam structure.

[0115] Based on the aforementioned discrete mesh model, the nodal degrees of freedom and field variable interpolation schemes for each element are further defined. Figure 3 The definition of nodal degrees of freedom is shown, and the nodal displacements of the beam element are also shown. Includes the axial displacement of each node i Lateral deflection and rotation angle The specific expression is:

[0116] (41)

[0117] The displacement field u at any point on the beam element can be expressed as:

[0118] (42)

[0119] in, Let be the shape function matrix of the displacement field. In the discrete model of the hybrid element beam structure of this invention, the physical region of the beam element in the thickness direction is divided into an upper region and a lower region, with the geometric neutral axis of the beam section as the boundary. In terms of spatial correspondence, each beam element is respectively arranged with an upper quadrilateral element in the upper region and a lower quadrilateral element in the lower region. Based on the displacement field of the Euler-Bernoulli beam considering axial deformation, the axial strain of the upper and lower regions of the beam is expressed as:

[0120] (43)

[0121] (44)

[0122] in, For the axial strain in the upper region of the beam; The axial strain in the lower layer region of the beam; Let be the derivative matrix of the displacement shape functions of the upper region of the beam; and Let U be the derivative matrix of the displacement shape function of the lower region of the beam; the superscripts U and D are used to characterize the upper and lower regions of the beam in the thickness direction, respectively.

[0123] The phase field is discretized using quadrilateral elements. Upper and lower quadrilateral elements are arranged along the thickness of the beam. The phase field at the nodes of the upper and lower quadrilateral elements are defined as follows:

[0124] (45)

[0125] (46)

[0126] in, The phase field of the upper quadrilateral element is the nodal phase field. The phase field of the lower quadrilateral element is the nodal phase field. The phase field of the i-th node of the upper quadrilateral unit; This represents the phase field of the i-th node in the lower quadrilateral element. The phase field of any point within the upper quadrilateral element is also shown. Phase field with any point within the lower quadrilateral unit The following can be obtained by interpolating the nodal phase field and the phase field shape function:

[0127] (47)

[0128] (48)

[0129] in, Let be the phase field shape function matrix. The phase field gradient can be expressed as:

[0130] (49)

[0131] (50)

[0132] in, The phase field gradient within the upper quadrilateral unit; The phase field gradient within the lower quadrilateral unit; The derivative matrix of the phase field shape function of the upper quadrilateral unit; Let be the derivative matrix of the phase field shape function of the lower quadrilateral unit.

[0133] Substituting formulas (41) to (44) into the control equation (30) established in this embodiment, the acceleration of the displacement field control equation is... Virtual displacement , adapt to change Axial stress The finite element discretization scheme is as follows:

[0134] (51)

[0135] in, For nodal acceleration; These are the displacements of virtual nodes; For the degenerate function of the upper region of the beam; These are the degeneracy functions for the lower region of the beam. Their expressions are:

[0136] (52)

[0137] (53)

[0138] (54)

[0139] (55)

[0140] in, , and These represent the axial displacements at the two nodes of the beam element, respectively. Lateral displacement and the degree of freedom of the cross section rotation angle The second derivative with respect to time. Substituting equation (51) into the displacement field governing equation in equation (30), we have:

[0141] (56)

[0142] Using the rules of vector dot product ( The scalar terms in equation (56) are then matrix-reorganized. The displacement field shape function interpolation terms are then... Convert to This will result in the virtual displacement of the nodes. Moving the integral sign to the outside, formula (56) can be rearranged as follows:

[0143] (57)

[0144] because It is arbitrary. For equation (57) to hold, the expression within the square brackets must be zero. The expression within the square brackets is defined as the element residual of the displacement field, and its specific expression is as follows:

[0145] (58)

[0146] Substituting formulas (45) to (50) into the phase field evolution control equation in formula (30) of this embodiment, the phase field in the phase field evolution control equation... Virtual Field Phase field gradient Virtual phase field gradient The finite element discretization schemes are as follows:

[0147] (59)

[0148] Since the phase field is discretized using a double-layer quadrilateral element mesh, the solution domain is... Geometrically decomposed into upper quadrilateral unit regions With the lower quadrilateral unit region Therefore, the integral term of the phase field evolution control equation in formula (30) can be split into two independent integrals for the upper quadrilateral unit and the lower quadrilateral unit. Substituting formula (59) into these integrals, we obtain the discrete phase field evolution control equations for the upper and lower quadrilateral units as formulas (60) and (61), respectively:

[0149] (60)

[0150] (61)

[0151] Using the rules of vector dot product ( The scalar terms in equations (60) and (61) are then matrix-reorganized. The phase field shape function interpolation terms are then... Convert to , Convert to This will enable the virtual phase field of the nodes. and Moving the integral sign to the outside, equations (60) and (61) are rearranged as follows:

[0152] (62)

[0153] Due to the virtual phase field of the node and It is arbitrary. For the above equation to hold, the expression within the curly braces must be zero. The expression within the curly braces is defined as the element residual of the phase field, and its specific expression is:

[0154] (63)

[0155] Based on formulas (58) and (63), the element residuals of the displacement field and the element residuals of the phase field are rearranged as follows:

[0156] (64)

[0157] in, For the element residuals of the displacement field; The displacement field shape function matrix; Mass density; For nodal acceleration; Young's modulus; For the degeneracy function of the lower region of the beam; Let be the derivative matrix of the displacement shape function of the lower region of the beam; Degeneracy function of the upper region of the beam; Let be the derivative matrix of the displacement shape function of the upper region of the beam; For nodal displacement; The phase field residual of the upper quadrilateral unit; The phase field shape function matrix; For historical variables; The critical energy release rate; For length dimensions; The derivative matrix of the phase field shape function of the upper quadrilateral unit; The phase field of the upper quadrilateral element is the nodal phase field. The phase field residual of the lower quadrilateral unit; The derivative matrix of the phase field shape function of the lower quadrilateral unit; The phase field of the lower quadrilateral element is represented by the nodal phase field; the superscripts U and D are used to characterize the upper and lower regions of the beam structure in the thickness direction, respectively; T is the matrix transpose symbol. For the computation of the domain volume integral; For the global computational domain; and These are the computational domains for the upper and lower regions of the beam, respectively. For the integral infinitesimal element of the global computational domain; For the integral infinitesimal element of the upper region; This is the integral element of the lower region.

[0158] The element residuals of the displacement field in formula (64) nodal displacement Find the partial derivative; phase field residual of the upper quadrilateral element. Phase field of the upper quadrilateral element nodes Find the partial derivative; phase field residual of the lower quadrilateral element. Phase field of the lower quadrilateral element nodes Find the partial derivatives. The element stiffness matrices for the displacement field and phase field are respectively:

[0159] (65)

[0160] in, Here is the element stiffness matrix of the displacement field; This is the phase field stiffness matrix of the upper quadrilateral element; This is the phase field stiffness matrix of the lower quadrilateral element.

[0161] Step 4: Using the Newmark method combined with the staggered algorithm, solve the displacement field and phase field in the discrete model until both the displacement field and phase field satisfy the convergence condition, and obtain the dynamic fracture phase field evolution results of the hybrid unit beam structure.

[0162] The staggered algorithm is specifically described as follows: (1) Solve the displacement field under the condition of a fixed phase field and update the historical variables to meet the irreversible condition of crack evolution; (2) Solve the phase field under the condition of a fixed displacement field; (1) and (2) are executed alternately until both the displacement field and the phase field meet the convergence condition.

[0163] In further detail, neglecting structural damping, the basic equations for the structural dynamic time history problem are as follows:

[0164] (66)

[0165] in, and These are the system's mass and stiffness matrices, respectively. and These are the displacement and acceleration responses of the system at any given time; It is the time-varying external load vector acting on the system; It is a time variable.

[0166] According to the Newmark method, the system in displacement of time step ,speed It can be obtained from the following formula:

[0167] (67)

[0168] in, , and Represented as the first Displacement, velocity, and acceleration of the system at each time step; , and Represented as the first Displacement, velocity, and acceleration of the system at each time step; This is the iteration time step; and These are the parameters used in the Newmark algorithm to ensure the unconditional stability of the numerical integration process.

[0169] Formula (67) is only used for initial value prediction at the beginning of the time step. During Newton's iteration process, when the displacement increment is obtained... Then, the acceleration and velocity need to be corrected according to the Newmark method. The correction formula is defined as follows:

[0170] (68)

[0171] in, for The acceleration of the system after the k-th iteration at each time step; for The first time step The system speed after the next iteration; for The first time step The acceleration of the system after the next iteration; for The first time step The system speed after the next iteration; For the first The displacement increment of the next iteration. Substituting equation (67) into equation (66) yields:

[0172] (69)

[0173] in This is the equivalent stiffness matrix; for The equivalent load vector at each time step. , The specific form is:

[0174] (70)

[0175] (71)

[0176] in, This is the quality matrix; for The external load vector at the time step; These are the parameters used in the Newmark algorithm to ensure the unconditional stability of the numerical integration process. This is the iteration time step. In this embodiment, it is set to... , .

[0177] The displacement field u and phase field in the governing equation (30) are solved using the Newton-Raphson iteration method. First, the residual is set at the current time step. Performing a first-order Taylor expansion at this point, ignoring higher-order terms, yields the linear equation:

[0178] (72)

[0179] in, Here is the tangent stiffness matrix; , Let the increment of the variable be the variable to be determined. Solve for the increment. Substitute it into the update formula In the middle, we get:

[0180] (73)

[0181] Define the unknown quantity as The corresponding global residual vector is Because of the use of an interleaved solution algorithm, It presents a segmented diagonal form, that is The defined vector of the unknown variable. and global residual vector and the block diagonal form of the tangent stiffness matrix Substituting into the general linearization formula (73) above, we can expand to obtain the incremental iteration format of the governing equation (30):

[0182] (74)

[0183] Based on the above incremental iteration format, displacement field at time step and phase field The staggered solution process is as follows: In the first... Each time step, assuming displacement ,speed acceleration Phase field d n and historical variables All are known. The first prediction is calculated using Newmark's prediction formula (67). Initial displacement prediction value of the step Then it enters an interleaved solution loop, in the first... In the first iteration, the first... Second phase field By utilizing the block diagonal properties of the stiffness matrix in formula (74), the linear equation corresponding to the displacement field is extracted. Solving this linear equation yields the first... displacement increment of the next iteration Using the update formula Calculate the first The displacement of the next iteration Subsequently, based on displacement increments Update the first one using formula (68) Speed ​​of each iteration and acceleration Based on the updated displacement Calculate the first using formulas (43) and (44). strain field of the next iteration Then, the historical variables are updated using formula (28). To ensure the irreversibility of crack evolution; then fix the displacement of the k-th iteration. Again, by utilizing the block diagonal characteristics of the stiffness matrix in formula (74), the linear equation corresponding to the phase field is extracted. Solving this linear equation yields the first... Phase field increment of the next iteration Using the update formula Calculate the first Phase field of the next iteration Continue this iterative process until the displacement field increment obtained in the current iteration step is reached. and phase field increment The Euclidean norms are all less than the preset tolerance. That is, when the convergence criterion is satisfied, the iteration stops, and the first iteration is obtained. displacement field at time step and phase field .

[0184] The convergence criterion for stopping the iteration is specifically stated as follows:

[0185] (75)

[0186] (76)

[0187] in, It is the Euclidean norm; This is the i-th component of the displacement increment; The i-th component of the phase field increment; This represents the total number of components.

[0188] Using the hybrid unit discrete model and numerical solution algorithm constructed above, the dynamic fracture process of the cantilever beam structure in this embodiment is simulated to obtain its fracture evolution results under dynamic impact load.

[0189] To verify the accuracy of the method of this invention, the calculation results of this embodiment are compared and analyzed with the comparative example of a traditional two-dimensional fracture phase-field model. In the comparative example, both the displacement field and the phase field are discretized using 400 quadrilateral elements, while the remaining material parameters and loading conditions are the same as in this embodiment.

[0190] The fracture evolution results of the cantilever beam under impact load are as follows: Figure 5 As shown, Figure 5 The gradual transition of color from blue to red indicates the progression of the structural material from an intact state to a fractured state. Figure 5 (a)-(d) are the calculation results of the method of the present invention; Figure 5(e)-(h) show the calculation results of the comparative traditional two-dimensional fracture phase-field model. The results show that the cantilever beam of the present invention fractures near the fixed constraint end, and the crack evolves along the beam thickness direction. The crack evolution trend of the present invention is consistent with the prediction results of the comparative traditional two-dimensional fracture phase-field model, indicating that the hybrid unit beam structure model proposed in this invention can accurately describe the fracture evolution behavior in the beam thickness direction. Furthermore, a comparative analysis of the energy response of the method of the present invention and the comparative traditional two-dimensional fracture phase-field model is conducted, such as... Figure 6 As shown. The kinetic energy, strain energy, and surface energy of the embodiments and comparative examples of the present invention are calculated using formulas (2), (3), and (4). The evolution results of the kinetic energy, strain energy, and surface energy obtained by the calculation method of the present invention are basically consistent with the calculation results of the traditional two-dimensional fracture phase field model of the comparative example, indicating that the method of the present invention has good accuracy in predicting the dynamic fracture energy evolution. Furthermore, the computational efficiency of the two methods of the present invention and the comparative example is compared, and their computation time is as follows: Figure 7 As shown, the computation time of the method of the present invention is 76 s, while the computation time of the traditional two-dimensional fracture phase field model in comparison is 147 s. The computation time of the present invention is reduced by about 48%, which shows that the present invention significantly improves computational efficiency while ensuring computational accuracy.

[0191] The above description is merely a preferred embodiment of the present invention and is not intended to limit the invention. Various modifications and variations can be made to the invention by those skilled in the art. Any modifications, equivalent substitutions, or improvements made to the present invention should be included within the scope of protection of the present invention.

Claims

1. A method for calculating the dynamic fracture phase field of a hybrid unit beam structure, characterized in that, Includes the following steps: Step 1: Determine the material parameters required for dynamic fracture phase field calculation of the hybrid unit beam structure; Step 2: Establish a dynamic fracture phase field mathematical model for the hybrid unit beam structure based on the Lagrange action principle. The mathematical model includes displacement field control equations and phase field evolution control equations. Step 3: Numerically discretize the mathematical model to establish a discrete model of the hybrid element beam structure: the displacement field is discretized using beam element mesh, and the phase field is discretized using quadrilateral element mesh; The discrete model is as follows: ； in, This is the quality matrix; Here is the displacement field stiffness matrix; For acceleration; For displacement field; This is the vector of external forces; Here is the phase field stiffness matrix; For phase field; The phase field driving force vector; Step 4: Using the Newmark method combined with the staggered algorithm, solve the displacement field and phase field in the discrete model until both the displacement field and phase field satisfy the convergence condition, and obtain the dynamic fracture phase field evolution results of the hybrid unit beam structure. The Newmark method is specifically described as follows: In the beam structure at the... The displacement, velocity, and acceleration at time step are determined by the first... The displacement, velocity, and acceleration of the beam structure at each time step are obtained as follows: ； in, , and The first Displacement, velocity, and acceleration of the beam structure at each time step; , and The first Displacement, velocity, and acceleration of the beam structure at each time step; For time steps; For time step; and Set the parameters used in the Newmark method to ensure the unconditional stability of the numerical integration process. , .

2. The method for calculating the dynamic fracture phase field of a hybrid unit beam structure according to claim 1, characterized in that, In step one, the material parameters required for the dynamic fracture phase field calculation of the hybrid unit beam structure are determined, specifically including: Young's modulus. Mass density Critical energy release rate Length dimension parameters , beam width Length of beam and beam thickness .

3. The method for calculating the dynamic fracture phase field of a hybrid unit beam structure according to claim 1, characterized in that, In step two, the Lagrange action is composed of kinetic energy, surface energy, strain energy, and work done by external forces, as specifically described below: ； in, Let Lagrange's action function be used. The kinetic energy of the system; The strain energy of the system; The surface energy of the system; This is work done by external force.

4. The method for calculating the dynamic fracture phase field of a hybrid unit beam structure according to claim 1, characterized in that, In step two, the mathematical model includes the displacement field control equation and the phase field evolution control equation, which are specifically described as follows: ； Where u is the displacement field; For acceleration; Mass density; Let be the axial stress at any point within the beam; Let be the axial strain at any point within the beam; and These represent the distributed loads applied to the beam in the axial and transverse directions, respectively. and These represent the axial and lateral displacements of points on the neutral axis of the beam, respectively. For phase field; For historical variables; The critical energy release rate; For length dimensions; The phase field gradient; For variational operators; This is a virtual displacement; For axial strain Arbitrary virtual strain; and These represent the virtual displacement components of the beam's neutral axis in the axial and transverse directions, respectively. It is a virtual phase field; The gradient of the virtual phase field; For the computational domain; It is a volume integral element; To compute the volume integral of the domain.

5. The method for calculating the dynamic fracture phase field of a hybrid unit beam structure according to claim 1, characterized in that, In step three, establishing the discretization model of the hybrid element beam structure includes discretizing the displacement field using beam elements and the phase field using quadrilateral elements; and arranging two layers of quadrilateral elements within the projection area of ​​each beam element in the thickness direction. Beam elements and quadrilateral elements exchange information between displacement fields and phase fields by establishing a geometric mapping relationship.

6. The method for calculating the dynamic fracture phase field of a hybrid unit beam structure according to claim 5, characterized in that, In step three, the specific steps for establishing the discrete model of the hybrid element beam structure are as follows: A one-dimensional beam element mesh is created along the axial direction of the beam structure. The nodal degrees of freedom of the beam element include axial displacement, lateral displacement, and cross-sectional rotation, to discretize the displacement field. For each beam element, upper and lower quadrilateral element meshes are sequentially established along its cross-sectional thickness direction. Both the upper and lower quadrilateral elements are set as two-dimensional planar quadrilateral elements and have only phase field degrees of freedom, to discretize the phase field. The upper and lower quadrilateral elements and the corresponding beam element are coupled by sharing the same Gaussian integration point to establish the coupling relationship between the displacement field and the phase field, thus obtaining the discrete model of the hybrid element beam structure.

7. The method for calculating the dynamic fracture phase field of a hybrid unit beam structure according to claim 6, characterized in that, The geometric mapping relationship between the beam element and the corresponding quadrilateral element can be expressed as follows: Let the i-th beam element be... Consisting of two adjacent beam nodes and The beam element corresponds to an upper quadrilateral element in the thickness direction. and a lower quadrilateral unit The corresponding quadrilateral element node numbering rules are defined as follows: ； in, Number the beam nodes; The upper quadrilateral element is the beam element, with its four nodes located at the top left, middle left, middle right, and top right positions respectively. The lower quadrilateral element corresponding to the beam element has four nodes located at the left center, lower left, lower right, and right center positions, respectively. For each beam element, its corresponding upper and lower quadrilateral elements together contain eight nodes. The bottom edge of the upper quadrilateral element coincides with the top edge of the lower quadrilateral element, and they share the same node number on the neutral axis. and Each beam element actually corresponds to six quadrilateral element nodes, with two nodes in total. Let the i-th beam node be... Located on the neutral axis, its coordinates are This is then mapped to three quadrilateral element nodes along the beam thickness direction: ； in, Let be the horizontal coordinate of the i-th beam node; Let be the vertical coordinate of the i-th beam node; This is the top edge node of the upper quadrilateral element; This is the bottom edge node of the lower quadrilateral unit; It is a common node located at the neutral axis, which simultaneously serves as the bottom edge node of the upper quadrilateral element and the top edge node of the lower quadrilateral element. The thickness of the beam.

8. The method for calculating the dynamic fracture phase field of a hybrid unit beam structure according to claim 1, characterized in that, In step four, the Newmark method combined with the interleaved algorithm solution process adopts a nested loop structure, that is, the Newmark method is used as the outer loop, and the interleaved iteration of the displacement field and the phase field is used as the inner loop. At the beginning of the outer loop at each time step, the predicted values ​​of displacement, velocity, and acceleration are first calculated using the Newmark method; then, the inner loop is entered, and the following steps are performed: (1) With the phase field fixed, solve the displacement field control equations to obtain the displacement field increment. And update the values ​​of displacement, velocity, and acceleration; (2) Update the historical variables based on the updated displacement field; (3) With the displacement field fixed, solve the phase field evolution control equation to obtain the phase field increment. ; Repeat steps (1) to (3) until the displacement field increment at that time step. and phase field increment If the Euclidean norm satisfies the convergence condition, the outer loop proceeds to the next time step; the inner loop is used to solve for the displacement field increment. and phase field increment The mathematical model is specifically expressed as follows: ； in, The stiffness matrix is ​​related to the displacement field; The stiffness matrix is ​​related to the phase field; This is the residual vector after discretizing the displacement field control equations; This is the residual vector after discretization of the phase field evolution control equations; This represents the displacement field increment; This is the phase field increment.

9. The method for calculating the dynamic fracture phase field of a hybrid unit beam structure according to claim 8, characterized in that, The historical variables are updated based on the updated displacement field, where the historical variables are the maximum tensile strain energy density experienced by the material during the loading history, i.e.: ； in, It is a position vector; For phase field; max is the maximum value operator; In response to the situation; The tensile strain energy density; For time variables, Time range At any time between.

10. The method for calculating the dynamic fracture phase field of a hybrid unit beam structure according to claim 1, characterized in that, The convergence condition is the displacement field increment in the current iteration step. and phase field increment The Euclidean norms are all less than the preset tolerance. ,Right now: ； ； in, It is the Euclidean norm; Let i be the i-th component of the displacement field increment; The i-th component of the phase field increment; This represents the total number of components.

Images