Phase field control equation numerical solution method
By employing the PDE module of COMSOL Multiphysics finite element software for numerical solution in the phase field method, the weak form of the phase field control equations is derived and discretized, thus solving the computational difficulty of the phase field method in multi-crack simulation and achieving efficient and stable numerical solution and multi-physics coupling.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- LANZHOU JIAOTONG UNIV
- Filing Date
- 2026-02-05
- Publication Date
- 2026-05-01
AI Technical Summary
The phase-field method suffers from computational convergence difficulties and high costs in simulating the initiation and propagation of multiple cracks, which limits its application in engineering.
The phase field problem was solved numerically using the PDE module in the COMSOL Multiphysics finite element software. The weak form of the phase field control equation was derived, and the displacement field and phase field were discretized using shape functions. The discretization of vector functions was achieved by combining the Helmholtz equation module. Boundary conditions and initialization parameters were strictly applied to complete the solution of the numerical model.
It significantly improves computational efficiency and numerical stability, enables multi-physics coupled simulation, and provides the advantage of dynamically and concurrently accessing the calculation results of other physics fields.
Smart Images
Figure CN121958725A_ABST
Abstract
Description
A numerical solution method for phase field control equations Technical Field
[0001] This invention belongs to the field of geological disaster prevention and control technology, and specifically relates to a numerical solution method for phase field control equations. Background Technology
[0002] The phase-field method, by introducing an order parameter to construct a system of differential equations and employing a diffuse interface description method, achieves continuous characterization of crack boundaries, avoiding the explicit tracking of crack surfaces required by traditional methods. This method automatically acquires crack propagation paths and spatial distributions through the spontaneous evolution of the order parameter, demonstrating significant advantages in simulating complex fracture behaviors such as crack initiation, propagation, bifurcation, and multi-crack penetration, providing an important technical means for studying the evolution mechanism of complex fracture systems. In the construction of the phase-field model, the phase-field variable 'd' is introduced to characterize the damage evolution process of materials or structures. Its theoretical framework originates from the variational principle of brittle fracture developed by Francfort and Marigo based on the Griffith fracture theory framework, while its numerical implementation relies on a regularized variational model. Compared to traditional damage mechanics theories, the phase-field model can more accurately predict crack propagation paths and penetration modes in materials. Compared to classical Griffith theory, this method addresses the singularity problem at the crack tip through regularization techniques, transforming the discrete crack interface into a continuous diffuse field description. This modeling approach can completely reproduce the entire process of a material from its initial intact state to crack initiation, propagation, merging, and eventual failure without requiring pre-defined crack propagation paths or additional fracture criteria. Although the phase-field method has become one of the mainstream methods for simulating crack behavior in recent years due to its unique advantages in simulating the initiation and propagation of multiple cracks, its numerical implementation faces significant challenges: the high nonlinearity of the governing equations, the strong coupling between phase-field variables and displacement fields, and mesh sensitivity issues lead to difficulties in computational convergence and high computational costs, which severely restricts the engineering application of this method. Summary of the Invention
[0003] To address the aforementioned problems, embodiments of the present invention propose a numerical solution method for the phase field control equations.
[0004] The numerical solution method for the phase field control equations of the present invention includes the following steps:
[0005] S1. Establish the phase-field control equations characterizing the fundamental physical processes;
[0006] S2. Derive the weak form of the phase field control equations;
[0007] S3. The displacement field and phase field in the phase field control equation are discretized using shape functions to obtain the weak domain forms of the displacement field equation and the phase field equation, respectively.
[0008] S4. Create a user phase-field model through the interactive interface of the PDE module builder;
[0009] S5. After discretizing the displacement field equation and phase field equation into vector functions using the PDE module, specify the associated coefficient field in the module;
[0010] S6. Solve the numerical model by strictly applying boundary conditions and initialization parameters.
[0011] S1 includes the following steps:
[0012] S201. Perform a weighted integral on both sides of the control equation;
[0013] S202. Multiply by the experimental function;
[0014] S203. Perform the integration over the entire computational domain and apply the divergence theorem.
[0015] The phase field control equation in S1 is:
[0016]
[0017] In the formula, Let be the spatial rate of change of the stress tensor. For stress tensor, Let d be the volume force, d be the phase field variable, H be the energy history function, and G be the energy history function. c The critical energy release rate. Parameters for controlling the degree of crack surface morphology propagation, It is the square of the spatial rate of change of the phase field.
[0018] The weak form of the phase field control equation obtained by S2 is:
[0019]
[0020] In the formula, For stress tensor, Let d be the strain tensor and d be the phase field variable. For displacement field test functions, For phase field test functions, The surface force acting on the boundary. Let dS be the volume force, and dS be the differential of the area where the boundary force acts. To solve for the differential of the solution domain, H is the energy history function, and G... c The critical energy release rate. Parameters for controlling the degree of crack surface morphology propagation, The spatial rate of change of the phase field. The spatial rate of change of the phase field experimental function. For the strain tensor test function, To solve for the domain, For surface force The boundaries of its function.
[0021] The weak domain form of the displacement field equation in S3 is:
[0022]
[0023] In the formula, This is the weak domain form of the displacement field equation, where d is the phase field variable. and Let Lamé constant be . and These represent the components of the lateral displacement along the x and y axes. and This represents the vertical displacement components along the x and y axes. and These are the components of the volume force along the x and y axes. and For experimental functions of lateral and vertical displacements, Let x be the experimental function for lateral displacement in the x-axis direction. Let be the experimental function for lateral displacement in the y-axis direction. Let x be the experimental function of the vertical displacement in the x-axis component. This is the experimental function for the vertical displacement in the y-axis component.
[0024] The weak domain form of the phase field equation in S3 is:
[0025]
[0026] In the formula, Let be the weak domain form of the phase-field equations, where d is the phase-field variable. G is the experimental function of the phase field variables. c The critical energy release rate. H is an energy history function used to control the degree of crack surface morphology diffusion.
[0027] The coefficients in S5 are the diffusion coefficient c and the source term f in the Helmholtz equation.
[0028] The initialization parameter in S6 is the displacement field. The phase field variable d, and the boundary conditions are the specified displacements in the x and y directions.
[0029] The beneficial effects of this invention are that the numerical solution method for the phase field control equations of this invention achieves a significant improvement in computational efficiency and numerical stability through innovative algorithm design; the use of the PDE module in the COMSOL Multiphysics finite element software to numerically solve the phase field problem has significant advantages in dynamically and concurrently accessing the calculation results of other physical fields, thus providing convenience for multi-physics coupled simulation. Attached Figure Description
[0030] Figure 1 is a configuration diagram of the three-point bending test of the semi-circular specimen of Embodiment 1 of this application: (a) model dimensions, (b) model mesh division.
[0031] Figure 2 is a comparative analysis of the load-displacement curves of the experimental measurement results and numerical simulation results of Embodiment 1 of this application.
[0032] Figure 3 is a comparative analysis of the failure modes of experimental observation and numerical simulation in Embodiment 1 of this application: (a) experimental results of Hou et al. (2021), (b) numerical calculation results of the model in this paper.
[0033] Figure 4 shows the crack propagation behavior and corresponding load curve of the sample of Example 1 of this application under load: (a) failure evolution characteristics, (b) load curve.
[0034] Figure 5 shows the evolution of the maximum shear stress, maximum principal stress, and displacement of the specimen in Example 1 of this application: (a) maximum principal stress, (b) maximum shear stress, and (c) displacement.
[0035] Figure 6 is a configuration diagram of specimens with different pre-fabricated cracks in Embodiment 2 of this application.
[0036] Figure 7 shows the failure evolution process and axial stress comparison of specimens with different pre-fabricated cracks in Example 2 of this application: (a) failure evolution of horizontal crack specimens, (b) failure evolution of 135° angle crack specimens, (c) failure evolution of vertical crack specimens, (d) failure evolution of 45° angle crack specimens, and (e) numerically calculated axial stress-strain curves.
[0037] Figure 8 shows the evolution of maximum shear stress in Example 2 of this application, including specimens with different pre-fabricated cracks: (a) horizontal crack specimen, (b) 135° angle crack specimen, (c) vertical crack specimen, and (d) 45° angle crack specimen. Detailed Implementation
[0038] Embodiments of the present invention are described in detail below, examples of which are illustrated in the accompanying drawings. The embodiments described below with reference to the accompanying drawings are exemplary and intended to explain the present invention, and should not be construed as limiting the present invention.
[0039] This application uses the PDE module in the COMSOL Multiphysics finite element software to numerically solve the phase field problem. This module has significant advantages in dynamically and concurrently accessing the calculation results of other physics fields, thus facilitating multiphysics coupled simulation.
[0040] The numerical solution method for the phase field control equations in this application includes the following steps:
[0041] S1. Establish the phase-field control equations characterizing the fundamental physical processes;
[0042] In the variational principle-based phase-field fracture theory framework, the total potential energy functional of the system is defined as a function of two independent field variables: the displacement field *u* and the phase field variable *d*. The displacement field describes the mechanical deformation of the material, while the phase field characterizes the damage evolution and crack propagation within the material; together, they constitute a complete physical picture describing the fracture problem. According to the variational principle of energy, the true physical state of the system (i.e., the equilibrium state) corresponds to the case where the total potential energy functional takes a stationary value. To obtain this equilibrium state, the displacement field must be considered separately. The phase field variable d is subjected to variational operations, i.e., a first-order variation is performed and its value is set to zero. This crucial mathematical process is physically equivalent to simultaneously satisfying the mechanical equilibrium condition and the thermodynamic driving force criterion for damage evolution. Through rigorous variational derivation, starting from a single energy functional, two sets of governing equations governing the displacement field distribution and phase field evolution are systematically derived, thus fully constructing the mathematical expression of the phase field governing equations:
[0043] In the formula, Let be the spatial rate of change of the stress tensor. For stress tensor, For volume force, d is the phase field variable, and G is the volume force. c The critical energy release rate. Parameters for controlling the degree of crack surface morphology propagation, Let be the square of the spatial rate of change of the phase field. It is the elastic strain energy.
[0044] For brittle materials, cracks during unloading cannot heal spontaneously, and their cumulative damage exhibits typical irreversible characteristics. This physical phenomenon fundamentally determines the mechanical response path of the material and constitutes a key thermodynamic constraint that must be considered in the constitutive model. To rigorously describe this irreversibility within the framework of phase-field fracture theory and avoid non-physical crack healing or damage recovery during the unloading phase, an energy history function H is introduced. The core function of this function is to record and lock the maximum tensile energy state reached by the material during the historical loading process. The expression for H is:
[0045] In the formula, H is the energy history function. It is the elastic strain energy. For time variables, For time intervals.
[0046] In the phase-field fracture model, the energy history function H is defined as the function H over the time interval [0, ... The maximum effective tensile strain energy experienced by the material is defined as H. This key definition mathematically captures the material's "memory effect," meaning it always remembers the highest energy state reached during the loading process. Since H records this unsurpassable peak energy potential, it naturally replaces the instantaneous elastic strain energy, becoming the true thermodynamic driving force for phase field evolution (i.e., damage propagation). Therefore, when updating the phase field evolution equations, the strain energy term driving the evolution of the phase field variable d is formally replaced by H. Based on this more physically reasonable correction, the governing equations of the phase-field model can be reformulated. The core difference lies in that the mechanical equilibrium equations are still derived from the variational application of the total potential energy to the displacement field, but the phase-field evolution equations are instead obtained by variationally applying the phase-field variables to a functional containing the historical maximum energy state. This reconstruction ensures that the model strictly adheres to the thermodynamic constraint of irreversible damage, thus accurately simulating the real physical behavior of materials such as rock where cracks only propagate and do not heal under cyclic loading. The revised phase-field governing equations are:
[0047] In the formula, Let be the spatial rate of change of the stress tensor. For stress tensor, Let d be the volume force, d be the phase field variable, H be the energy history function, and G be the energy history function. c The critical energy release rate. Parameters for controlling the degree of crack surface morphology propagation, It is the square of the spatial rate of change of the phase field.
[0048] S2. Derive the weak form of the phase field control equations;
[0049] The core of this numerical calculation framework includes two coupled field models: the stress field model and the phase field model. To successfully implement finite element calculations, their strong-form governing equations must be transformed into weak forms that are more suitable for numerical solutions. The essence of this transformation is to reduce the continuity requirements of the solution function and to transform the complex differential equations into computable integral forms. To derive the weak forms of the stress field model and the phase field model, the following variational operations need to be performed: (1) Weighted integration of both sides of the governing equations: The essence of this is to adopt the idea of the weighted residual method. The strong-form governing equations require that they be strictly true at every point in the computational domain, which is almost impossible to achieve numerically. By performing weighted integration, the strict "point-to-point" satisfaction is weakened to satisfaction in the sense of integral average, which is the fundamental first step from strong form to weak form. (2) Multiply by the experimental function (weight function): The experimental function is the shape function in the Galerkin framework. The mathematical essence of this operation is to construct a function space and check whether the residuals of the governing equations are orthogonal in this space. Physically, it can be understood as introducing virtual displacement (corresponding to the mechanical equilibrium equation) or virtual phase field variable (corresponding to the phase field evolution equation) into the principle of virtual work, thereby expressing the physical equilibrium condition as an energy orthogonal form. (3) Complete the domain integration over the entire computational domain Ω and apply the divergence theorem: This is the most crucial step in the derivation process. By integrating by parts (i.e. applying the divergence theorem), the higher-order derivatives acting on the trial function in the equation are transferred to the trial function. The mathematical essence of this operation is to reduce the smoothness requirement of the trial function, thereby allowing the use of simple piecewise polynomial functions as shape functions. Its physical essence is to naturally introduce natural boundary conditions in the weak form, so that physical conditions such as surface force and boundary flux can be automatically satisfied and become part of the equation.
[0050] First, the phase field control equations are all multiplied by a unified experimental function. We can obtain:
[0051] In the formula, To unify the form of the test function, Let be the spatial rate of change of the stress tensor. For stress tensor, Let d be the volume force, d be the phase field variable, H be the energy history function, and G be the energy history function. c The critical energy release rate. Parameters for controlling the degree of crack surface morphology propagation, Let be the square of the spatial rate of change of the phase field. To solve for the differential of the domain, To find the solution domain.
[0052] Then, by applying the divergence theorem and based on the principles of virtual work and minimum potential energy, the weak form of the phase field control equations is derived as follows:
[0053] In the formula, For stress tensor, Let d be the strain tensor and d be the phase field variable. For displacement field test functions, For phase field test functions, For surface forces acting on the boundary, Let dS be the volume force, and dS be the differential of the area where the boundary force acts. To solve for the differential of the solution domain, H is the energy history function, and G... c The critical energy release rate. Parameters for controlling the degree of crack surface morphology propagation, The spatial rate of change of the phase field. The spatial rate of change of the phase field experimental function. For the strain tensor test function, To solve for the domain, The boundary is where the surface force t acts.
[0054] S3. The displacement field and phase field in the phase field control equation are discretized using shape functions to obtain the weak domain forms of the displacement field equation and the phase field equation, respectively.
[0055] Equation (5) has the exact same solution as the original governing equations. Transforming the differential equations into a weak form means that the requirements for the solution are relaxed—the solution does not need to satisfy the equations at every point, but only needs to satisfy equilibrium in the integral sense over the finite field. This fundamental property allows the solution to exist in discrete form, constituting one of the basic principles of the finite element method. Shape functions are used to represent the displacement field in equation (5). Discretizing the phase field variable d, its expression is:
[0056] In the formula, and Let these represent the element shape function matrices corresponding to the displacement field and phase field, respectively. and These represent the nodal value vectors of the displacement field and phase field at the element nodes, respectively.
[0057] For plane problems, the expression for the shape function is:
[0058]
[0059] In the formula, n is the number of nodes in each finite element element.
[0060] The gradients of the displacement field and phase field variables can be defined as:
[0061]
[0062] In the formula, and These represent the nodal value vectors of the displacement field and phase field at the element nodes, respectively. and The derivatives of the shape functions representing the displacement field and phase field variables are expressed mathematically as follows:
[0063] In the formula, n is the number of nodes in each finite element element.
[0064] Using the Galerkin finite element method, displacement field The expressions for the phase field variable d and its gradient are as follows:
[0065] In the formula, and Let represent the nodal value vector variation of the displacement field and phase field at the element nodes, respectively. and Let these represent the element shape function matrices corresponding to the displacement field and phase field, respectively. and The derivatives of the shape functions representing the displacement field and phase field variables. It is an experimental function of the spatial rate of change of the phase field.
[0066] Substituting equation (10) into equation (5), we can obtain the discrete forms of the displacement field equation and the phase field equation as follows:
[0067] (11)
[0068] In the formula, This represents the degenerate elastic modulus tensor. and Let represent the nodal value vector variation of the displacement field and phase field at the element nodes, respectively. and Let these represent the element shape function matrices corresponding to the displacement field and phase field, respectively. and The derivatives of the shape functions representing the displacement field and phase field variables. and These represent the nodal value vectors of the displacement field and phase field at the element nodes, respectively, where d is the phase field variable. The surface force acting on the boundary. Let dS be the volume force, and dS be the differential of the area where the boundary force acts. To solve for the differential of the solution domain, H is the energy history function, and G... c The critical energy release rate. Parameters for controlling the degree of crack surface morphology propagation, To solve for the domain, Let T be the boundary where the surface force t acts, and let T denote the transpose of the matrix.
[0069] This application employs the theoretical framework of the hybrid phase-field method. The displacement field is governed by isotropic stress-strain constitutive relations. It can be calculated using the following formula:
[0070] In the formula, Let d represent the initial elastic tensor of the complete material, and d be the phase field variable.
[0071] Since equation (11) applies to any test function, the discrete forms of the displacement field and phase field equations can be derived as follows:
[0072]
[0073] In the formula, This represents the degenerate elastic modulus tensor. and Let these represent the element shape function matrices corresponding to the displacement field and phase field, respectively. and The derivatives of the shape functions representing the displacement field and phase field variables. The vector represents the nodal values of the displacement field at the element nodes, where d is the phase field variable. The surface force acting on the boundary. Let dS be the volume force, and dS be the differential of the area where the boundary force acts. To solve for the differential of the solution domain, H is the energy history function, and G... c The critical energy release rate. Parameters for controlling the degree of crack surface morphology propagation, To solve for the domain, Let T be the boundary where the surface force t acts, and let T denote the transpose of the matrix.
[0074] The discrete forms of the displacement field and phase field equations in equation (13) need to be input in a weak form in the PDE module of the COMSOL numerical analysis platform. Then, the global parameter definition function and the transient iterative calculation program are called to realize the coupling of each physical field.
[0075] The displacement field equations are implemented using a solid mechanics module, and the constitutive equations of the model are as follows:
[0076]
[0077] In the formula, For stress tensor, Let d represent the initial elastic tensor of the complete material, and d be the phase field variable. For strain tensor.
[0078] For displacement field Define COMSOL test functions and Then the domain weak form of the displacement field equation is:
[0079]
[0080] In the formula, This is the weak domain form of the displacement field equation, where d is the phase field variable. and Let Lamé constant be . and These represent the components of the lateral displacement along the x and y axes. and These are the components of the vertical displacement along the x and y axes. and These are the components of the volume force along the x and y axes. and For experimental functions of lateral and vertical displacements, Let x be the experimental function for lateral displacement in the x-axis direction. Let be the experimental function for lateral displacement in the y-axis direction. Let x be the experimental function of the vertical displacement in the x-axis component. This is the experimental function for the vertical displacement in the y-axis component.
[0081] The phase-field equations are implemented using the Helmholtz equation module, and the constitutive equations of the model are:
[0082]
[0083] In the formula, c is the diffusion coefficient of the Helmholtz equation, and f is the source term of the Helmholtz equation.
[0084] The corresponding parameters in the model are:
[0085]
[0086] For the phase field variable d, define the COMSOL test function. Then the weak domain form of the phase field equation is:
[0087] In the formula, Let be the weak domain form of the phase-field equations, where d is the phase-field variable. G is the experimental function of the phase field variables.c The critical energy release rate. H is an energy history function used to control the degree of crack surface morphology diffusion.
[0088] S4. Create a user phase-field model through the interactive interface of the PDE module builder;
[0089] By creating a user phase field model for equations (14), (15), (16), (17), and (18) through the interactive interface of the PDE module builder, and filling in the specified associated coefficient field in the module according to the discretization function of the PDE module, the phase field can be solved.
[0090] S5. After discretizing the displacement field and phase field equations into vector functions using the PDE module, specify the associated coefficient field in the module; the coefficients are mainly the diffusion coefficient c and the source term f in the Helmholtz equation.
[0091] S6. The numerical model is solved by strictly applying boundary conditions and initialization parameters. The initialization parameters are the displacement field. The phase field variable d, and the boundary conditions are the specified displacements in the x and y directions.
[0092] Example 1
[0093] To verify the rationality of the numerical solution method, this application uses three-point bending test data from Hou et al. for comparison and verification. The test configuration is shown in Figure 1. The semi-circular specimen (radius 50 mm) contains a pre-fabricated crack located at the center. The crack's geometric dimensions are as follows: width 0.5 mm, length 13.5 mm, and inclination angle 40°. The calculation model parameters are: the material is isotropic, and the initial elastic tensor of the intact material is... The axial stiffness coefficient is 11.41 GPa, the Poisson's ratio is 0.35, and the critical energy release rate G c =1200 J / m², a parameter controlling the degree of crack surface morphology diffusion. .
[0094] Figure 2 shows a comparative analysis of the load-displacement curves obtained from experimental measurements and numerical simulations of this application, demonstrating good agreement between the two sets of data. During the initial loading stage, both the experimental and simulated curves exhibit linear elastic behavior with nearly identical slopes. When the load exceeds the 800 kN threshold, the experimental curve shows significant nonlinear hardening characteristics, while the yield load predicted by the numerical simulation is 760 kN (underestimated by 5%). The experimentally measured hardening modulus is 6.5 kN / mm², while the simulated value is 5.2 kN / mm² (hardening slope reduced by 20%). Energy dissipation analysis indicates that the area enclosed by the experimental curve is 15% larger than the numerical simulation prediction, suggesting an additional energy dissipation mechanism in the actual specimen that the computational model failed to fully capture.
[0095] Figure 3 illustrates a comparative analysis of the failure modes observed experimentally and those simulated in this application. Experimental results show that the specimen failure primarily propagates along the direction of the maximum principal stress, but exhibits path fluctuations due to material inhomogeneity (maximum deflection angle reaching 15°), and significant path deflection occurs upon interaction with pre-existing defects. The fracture surface displays obvious roughness characteristics, accompanied by clear local crushing zones. In contrast, the failure mode predicted by the numerical simulation is dominated by a single principal crack, with a smoother propagation path that strictly adheres to the maximum circumferential stress criterion, failing to reproduce the microstructural complexity observed in the experiment.
[0096] Figure 4 illustrates the fundamental relationship between crack propagation behavior and the corresponding stress-strain response under load. In stage A, microcracks nucleate at inherent defects in the material, with limited propagation length and slow speed, during which the material maintains elastic deformation. The stable propagation stage (B→C) is characterized by crack extension along the direction of maximum principal stress or weak interface, with the formation of a plastic zone at the crack tip accompanied by energy dissipation. Stage D (instability fracture) is characterized by accelerated crack propagation, formation of macroscopic fracture surfaces, and eventual instability.
[0097] Figure 5 illustrates the evolution of maximum shear stress, maximum principal stress, and displacement in a semi-circular bending specimen. Based on the characteristic points (AD) of the load-displacement curve of the semi-circular bending specimen, the evolution of its maximum principal stress, maximum shear stress, and displacement follows a clear mechanical process. Elastic loading stage (A→B): The maximum principal stress and maximum shear stress increase linearly with displacement, dominated by the stress concentration effect at the notch tip; the displacement exhibits proportional elastic deformation characteristics. Crack initiation stage (B): When the tensile strength is reached, microcrack formation leads to a deceleration in the growth of the maximum principal stress; while the development of the shear band promotes a nonlinear acceleration in the growth of the maximum shear stress; the displacement deviates from the linear trajectory, marking the onset of plastic deformation. Peak state stage (C): This stage corresponds to the critical displacement: under the combined effect of tensile crack penetration and shear localization, the maximum principal stress and maximum shear stress reach their peak values simultaneously, indicating the occurrence of fracture instability. Post-failure stage (D): The maximum principal stress drops sharply to near zero due to the loss of tensile strength; the maximum shear stress stabilizes at the residual stress level along the frictional slip along the fracture surface.
[0098] Example 2
[0099] The computational method of this application was used to study the crack propagation behavior of four multi-crack numerical specimens. The specimens were 100 mm high and 50 mm wide. Each model contained four pre-fabricated crack configurations with identical geometric dimensions (12 mm long, 0.5 mm wide) but different inclination angles (0°, 135°, 90°, and 45°). The cracks were uniformly distributed within each specimen, and the specific crack spacing configuration is shown in Figure 6. Uniaxial compression numerical simulations were performed on the four specimens. The model parameters were: isotropic material, initial elastic tensor of intact material. The axial stiffness coefficient is 35.95 GPa, the Poisson's ratio is 0.33, and the critical energy release rate G c =1320 J / m², a parameter controlling the degree of crack surface morphology diffusion. .
[0100] Figure 7 illustrates the failure evolution process and corresponding stress-strain curves of specimens with different pre-existing cracks. Analysis shows that specimen 1 exhibits five typical evolution stages: Stage A is characterized by initial elastic deformation, with a linear stress-strain relationship and no macroscopic damage; Stage B begins with the direct initiation of a main crack upon reaching the yield strength, which propagates stably from the stress concentration zone and sustains a slow increase in load; Stage C marks the transition to unstable crack propagation, where rapid bifurcation or zigzag propagation occurs due to the energy release rate exceeding a critical value, accompanied by the formation of secondary cracks; Stage D achieves full-section crack penetration, forming a macroscopic fracture surface and dividing the specimen into multiple segments; Stage E represents complete instability and the complete loss of load-bearing capacity. This staged evolution system systematically elucidates the entire mechanical degradation process from damage initiation to catastrophic failure.
[0101] The evolution characteristics of specimen 2 are as follows: Stage A is the elastic deformation phase, with a linear stress-strain relationship and no macroscopic damage; Stage B is the stable propagation of the main crack from the stress concentration zone, with the load-bearing capacity maintained or slightly increased; Stage C is the transition to unstable propagation, manifested by the bifurcation of the main crack and / or the initiation of multiple secondary cracks; Stage D is the full-section penetration of the main crack to form a macroscopic fracture surface, dividing the specimen into independent blocks; Stage E is the complete loss of load-bearing capacity.
[0102] The failure process of specimen 3 includes: stage A is the initial elastic deformation phase accompanied by micro-damage nucleation; stage B enters the yielding stage, characterized by the initiation of the main crack, significant local plastic deformation and crack propagation along the direction of the maximum principal stress; stage C achieves stable and accelerated propagation; stage D exhibits rapid propagation and bifurcation; stage E completes the global crack penetration of the specimen and the final fracture surface formation.
[0103] The mechanical response of specimen 4 is divided into the following stages: Stage A is the micro-damage incubation period, during which elastic strain energy continues to accumulate; Stage B initiates nonlinear deformation and is accompanied by the nucleation of the main crack; Stage C is the damage localization phase, during which the main crack achieves breakthrough propagation by crossing the process zone; Stage D reflects instability-driven failure, dominated by Type I cracks, with the participation of potential Type II / III cracks and significant dynamic bifurcation effects; Stage E achieves complete disintegration of the specimen through the establishment of the final fracture surface.
[0104] Figure 8 illustrates the evolution of maximum shear stress in specimens with different pre-fabricated fractures. It can be seen that at point A, the shear stress level is controlled by the angle between the fracture and the principal stress. The results show that the 45° inclined fracture exhibits the highest value (15 MPa) due to the maximum shear stress concentration effect, while the horizontal fracture shows the lowest stress (3 MPa) due to the stress shielding effect. At point E, the residual shear stress is determined by both the failure mode and the friction mechanism. The results show that the 45° inclined specimen maintains the highest value (12 MPa) through coordinated frictional slip along the shear fracture surface, while the specimen with vertical fractures relies on the interlocking effect of the uneven body generated by tensile-dominated failure, resulting in a lower residual stress (8 MPa). This evolution pattern indicates that rock masses with near-45° fractures pose a dual risk: high initial shear stress sensitivity and continuous slip potential. Therefore, it is crucial to focus on controlling this risk in engineering stability assessments to mitigate progressive failure.
[0105] In this invention, the terms "one embodiment," "some embodiments," "example," "specific example," or "some examples," etc., refer to a specific feature, structure, material, or characteristic described in connection with that embodiment or example, which is included in at least one embodiment or example of the invention. In this specification, the illustrative expressions of the above terms do not necessarily refer to the same embodiment or example. Furthermore, the specific features, structures, materials, or characteristics described may be combined in any suitable manner in one or more embodiments or examples. Moreover, without contradiction, those skilled in the art can combine and integrate the different embodiments or examples described in this specification, as well as the features of different embodiments or examples.
[0106] Although the above embodiments have been shown and described, it is understood that the above embodiments are exemplary and should not be construed as limiting the present invention. Any changes, modifications, substitutions and variations made to the above embodiments by those skilled in the art are within the protection scope of the present invention.
Claims
1. A numerical solution method for phase-field control equations, characterized in that, Includes the following steps: S1. Establish the phase-field control equations characterizing the fundamental physical processes; S2. Derive the weak form of the phase field control equations; S3. The displacement field and phase field in the phase field governing equations are discretized using shape functions to obtain the domain weak forms of the displacement field equations and phase field equations, respectively; the domain weak form of the displacement field equation is: In the formula, This is the weak domain form of the displacement field equation, where d is the phase field variable. and Let Lamé constant be . and These represent the components of the lateral displacement along the x and y axes. and This represents the vertical displacement components along the x and y axes. and These are the components of the volume force along the x and y axes. and For experimental functions of lateral and vertical displacements, Let x be the experimental function for lateral displacement in the x-axis direction. Let be the experimental function for lateral displacement in the y-axis direction. Let x be the experimental function of the vertical displacement in the x-axis component. Let be the experimental function of the vertical displacement in the y-axis component; the domain weak form of the phase field equation is: In the formula, This is the field weak form of the phase field equation. G is the experimental function of the phase field variables. c The critical energy release rate. H is an energy history function used to control the degree of crack surface morphology diffusion. S4. Create a user phase-field model through the interactive interface of the PDE module builder; S5. After discretizing the displacement field equation and phase-field equation into vector functions using the PDE module, specify the associated coefficient field in the module; the coefficients include the diffusion coefficient in the Helmholtz equation, source term: In the formula, c is the diffusion coefficient of the Helmholtz equation, and f is the source term of the Helmholtz equation; S6. The numerical model is solved by strictly applying boundary conditions and initialization parameters.
2. The numerical solution method for the phase-field control equations according to claim 1, characterized in that, S1 includes the following steps: S201. Perform a weighted integral on both sides of the control equation; S202. Multiply by the test function; S203. Perform the integral over the entire computational domain and apply the divergence theorem.
3. The numerical solution method for the phase-field control equations according to claim 2, characterized in that, The phase field control equation in S1 is: In the formula, Let be the spatial rate of change of the stress tensor. For stress tensor, Let d be the volume force, d be the phase field variable, H be the energy history function, and G be the energy history function. c The critical energy release rate. Parameters for controlling the degree of crack surface morphology propagation, It is the square of the spatial rate of change of the phase field.
4. The numerical solution method for the phase-field control equations according to claim 2, characterized in that, The weak form of the phase field control equation obtained by S2 is: In the formula, For stress tensor, Let d be the strain tensor and d be the phase field variable. For displacement field test functions, For phase field test functions, For surface forces acting on the boundary, Let dS be the volume force, and dS be the differential of the area where the boundary force acts. To solve for the differential of the solution domain, H is the energy history function, and G... c The critical energy release rate. Parameters for controlling the degree of crack surface morphology propagation, The spatial rate of change of the phase field. The spatial rate of change of the phase field experimental function. For the strain tensor test function, To solve for the domain, For surface force The boundaries of its function.
5. The numerical solution method for the phase-field control equations according to claim 1, characterized in that, The initialization parameter in S6 is the displacement field. The phase field variable d, and the boundary conditions are the specified displacements in the x and y directions.