CN116486953B - Fracture phase field simulation method containing microstructure effect - Google Patents

Fracture phase field simulation method containing microstructure effect Download PDF

Info

Publication number
CN116486953B
CN116486953B CN202310441617.6A CN202310441617A CN116486953B CN 116486953 B CN116486953 B CN 116486953B CN 202310441617 A CN202310441617 A CN 202310441617A CN 116486953 B CN116486953 B CN 116486953B
Authority
CN
China
Prior art keywords
strain
phase field
fracture
equation
displacement
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202310441617.6A
Other languages
Chinese (zh)
Other versions
CN116486953A (en
Inventor
于红军
黄灿杰
朱帅
郝留磊
张迎宾
刘士状
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Harbin Institute of Technology
Original Assignee
Harbin Institute of Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Harbin Institute of Technology filed Critical Harbin Institute of Technology
Priority to CN202310441617.6A priority Critical patent/CN116486953B/en
Publication of CN116486953A publication Critical patent/CN116486953A/en
Application granted granted Critical
Publication of CN116486953B publication Critical patent/CN116486953B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • GPHYSICS
    • G16INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16CCOMPUTATIONAL CHEMISTRY; CHEMOINFORMATICS; COMPUTATIONAL MATERIALS SCIENCE
    • G16C60/00Computational materials science, i.e. ICT specially adapted for investigating the physical or chemical properties of materials or phenomena associated with their design, synthesis, processing, characterisation or utilisation
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00Computer-aided design [CAD]
    • G06F30/20Design optimisation, verification or simulation
    • G06F30/23Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14Force analysis or force optimisation, e.g. static or dynamic forces
    • YGENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02TCLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Computing Systems (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Investigating Strength Of Materials By Application Of Mechanical Stress (AREA)

Abstract

The invention discloses a fracture phase field simulation method comprising a microstructure effect, which considers the influence of micro-bending deformation and micro-torsion deformation of a microstructure on the expansion behavior of a macroscopic crack. The thermodynamic basic frame and the micropolar elasticity theory based on the fracture phase field method are used for deriving the total potential energy expression. In the total potential energy expression, considering different sources of crack expansion driving force in the mixed fracture mode, the expansion driving force of the mixed mode crack is obtained by decomposing the total elastic strain energy based on the decomposition of the symmetrical strain and the antisymmetric strain of the total strain and the tensile and compressive anisotropic decomposition of the elastic strain energy related to the symmetrical strain. Based on the method, various fracture phase field models containing the microstructure effect can be derived, so that a fracture phase field model with strong universality containing the microstructure effect is established. The method has very important significance for researching scale-related behaviors of damage and fracture failure processes of the microstructure-containing material and identifying parameters in damage identification of the microstructure-containing material.

Description

Fracture phase field simulation method containing microstructure effect
Technical Field
The invention belongs to the technical field of damage and fracture mechanics, relates to a fracture phase field simulation method, and in particular relates to a method for simulating and solving crack propagation behaviors and load response of a microstructure-containing material under the action of complex load.
Background
Many natural and man-made materials, such as bones, rocks with micro-porous structures, foamed metal-like materials, lightweight ceramic materials, lightweight concrete, chiral metamaterials, etc. have complex microstructure features in the macro, fine, and microscopic scale ranges. Such materials generally have light weight, high strength and behavior exhibiting a microstructure effect under the action of multiple fields, are widely used in the fields of science and engineering, such as medical instruments, space shuttles, space stations, large bridges, automobiles, etc., and increasingly exhibit excellent properties thereof.
However, one of the common problems with microstructure-containing materials is that the materials are brittle and have a low fracture toughness. The materials are easy to damage and break failure due to various defects, such as microcracks, micro holes or dislocation, which are inevitably introduced in the preparation, forming, assembly and service processes. Therefore, a reasonable damage and fracture simulation model and a simulation method are established for the materials, and the method is particularly necessary for researching damage and fracture failure processes of the materials.
Damage and failure to fracture of microstructure-containing materials is typically the result of the combined processes of nucleation of microcracks at defects, the aggregation and penetration of numerous microcracks to form a major crack, crack propagation, and bifurcation.
The fracture phase field method is based on the theory of continuous medium damage mechanics and classical fracture mechanics, and is a high-efficiency simulation method for simulating crack initiation, propagation and bifurcation processes. The method has the advantages of clear mechanical concept and attractive mathematical form, does not need to introduce any damage and fracture criteria in the simulation process, and can automatically capture the nucleation, expansion and bifurcation processes of cracks. However, the current fracture phase field method lacks analysis and discussion of microstructure effects, and a more general fracture phase field thermodynamic framework has not yet been developed. Simulation methods based on brittle fracture phase field and micropolar elasticity theory have not been further developed.
Disclosure of Invention
In order to solve the defects existing in the prior art, the invention provides a fracture phase field simulation method containing a microstructure effect aiming at the problems of complex damage and fracture failure modes of a microstructure-containing material. The invention considers that the classical elasticity theory can not represent the microstructure effect of the material, the fracture phase field theory framework based on the micropolar elasticity theory is not completely mature, and provides a fracture phase field simulation method containing the microstructure effect based on a thermodynamic unified framework.
The invention aims at realizing the following technical scheme:
a fracture phase field simulation method comprising a microstructure effect, comprising the steps of:
step one, considering the influence of micro bending freedom and micro torsion freedom in a microstructure, introducing a micro-polar elasticity theory to describe the elasticity mechanism of the material; in order to solve the problem of crack propagation of a mixed fracture mode, decomposing total strain into two parts of symmetrical strain and anti-symmetrical strain, and establishing an elastic strain energy expression corresponding to the symmetrical strain and the anti-symmetrical strain;
step two, considering the crack driving force of tensile and compressive anisotropy in the mixed fracture mode, further carrying out tensile and compressive decomposition on the elastic strain energy related to the symmetrical strain, and obtaining the expression of the elastic strain energy of each part; establishing an expression of a total potential energy functional based on a thermodynamic basic frame of a fracture phase field method and decomposition of elastic strain energy;
thirdly, deducing a control equation in a partial differential form from the total potential energy functional expression based on an energy transformation principle, wherein the control equation comprises a stress balance equation, an even stress balance equation and a fracture phase field driving force equation; aiming at the irreversible evolution process of crack growth, adopting the maximum value of strain energy history until the calculation time as the driving force of crack growth;
step four, a control equation based on a partial differential form is equivalently converted into an integral equation of a weak form by means of a virtual work principle;
fifthly, carrying out numerical discrete on a displacement field, a micro-rotation field and a fracture phase field based on a finite element method, carrying out localized treatment on a weak form integral equation at a cell, and deducing a cell stiffness equation;
and step six, forming a nonlinear total stiffness equation by the cell stiffness equation according to a cell assembly rule, and designing an iterative solution algorithm of the nonlinear total stiffness equation.
Compared with the prior art, the invention has the following advantages:
1. the invention provides a processing method based on a fracture phase field unified thermodynamic frame, and a certain mixed mode fracture phase field simulation method containing a microstructure effect can be formed based on the method according to a specific design of a geometric function and an energy degradation function.
2. The invention can directly characterize the microstructure effect in constitutive relation by considering the influence of micro-structure microbending and micro-torsion degrees of freedom in microstructure-containing materials. The fracture phase field model provided by the invention comprises fracture phase field methods based on the micropolar elasticity theory and the correction even stress theory, can be degenerated into various classical fracture phase field methods, and greatly expands the application range of the fracture phase field method.
3. The invention can analyze the influence of the microstructure effect on the macroscopic crack expansion behavior (comprising crack expansion morphology and load-displacement curve) through the arrangement of the material properties, load boundary conditions and geometric configuration of the microstructure-containing material.
4. The method has good applicability and stability, and can be combined with the existing calculation methods such as a finite element method, an extended finite element method, a grid-free Galerkin method and the like to develop a commercial program so as to flexibly adapt to the change of the required analysis problem.
5. The method has potential scientific research and application values of further development, such as taking the effects of temperature and chemical fields into consideration and introducing the influence of other physical fields. Aiming at the material containing the microstructure, a quasi-static and dynamic fracture phase field simulation method comprehensively considering various physical field coupling effects is formed.
6. The method has very important significance for researching scale-related behaviors of damage and fracture failure processes of the microstructure-containing material and identifying parameters in damage identification of the microstructure-containing material.
Drawings
FIG. 1 is a flow chart of a computational analysis of crack propagation behavior of a microstructure-containing material;
FIG. 2 is a computational analysis process of an alternating solution;
FIG. 3 is a schematic diagram of a displacement field, micro-rotation field and phase field multi-layer element solution based on an alternating solutionRealizing a frame: the first layer is a fracture phase field solving unit, wherein the history field is determined by a displacement field and a micro-rotation field obtained in the last increment step; the second layer is a displacement field and micro-rotation field solving unit, wherein the fracture phase field is obtained by the current increment step; the third layer takes the rigidity value to approach to a smaller value (the elastic modulus is 10 -11 Unit) for result display;
FIG. 4 is a quarter-wave plate model and its geometric boundaries for an infinite microplate with a central circular aperture;
fig. 5 is a graph of the stress component of an infinite plate with a central hole as a function of distance (analytical solution versus numerical solution): (a) Stress component sigma (b) stress component sigma θr
FIG. 6 is a grid, geometric boundary and loading condition pre-run of a prefabricated edge crack plate: (a) a tensile displacement test; (b) shear displacement testing;
FIG. 7 is a graph showing load-displacement curves of a preformed edge crack plate under tensile and shear loads: (a) a tensile displacement test; (b) shear displacement testing;
FIG. 8 is a graph showing crack growth morphology contrast of a preformed edge crack plate under tensile and shear loading: (a) the crack propagation morphology of the stretched plate obtained by the method; (b) the crack propagation morphology of the sheared plate obtained by the method; (c) The stretched plate crack propagation morphology obtained in literature (Patil, r.u., CMAME, 2018); (d) The resulting sheared-sheet crack propagation morphology of literature (Patil, r.u., CMAME, 2018);
FIG. 9 is a graph of load displacement of a preformed edge crack plate containing a microstructure effect under tensile and shear loading: (a) a tensile displacement test; (b) shear displacement testing;
FIG. 10 is a crack propagation morphology of a preformed edge crack plate containing a microstructure effect under tensile and shear loading: (a) - (e) is a tensile test, and (f) - (j) are shear tests.
Detailed Description
The following description of the present invention is provided with reference to the accompanying drawings, but is not limited to the following description, and any modifications or equivalent substitutions of the present invention should be included in the scope of the present invention without departing from the spirit and scope of the present invention.
The invention provides a fracture phase field simulation method comprising a microstructure effect, which considers the influence of micro-bending deformation and micro-torsion deformation of a microstructure on the expansion behavior of a macroscopic crack. The thermodynamic basic frame and the micropolar elasticity theory based on the fracture phase field method are used for deriving the total potential energy expression. In the total potential energy expression, considering different sources of crack expansion driving force in the mixed fracture mode, the expansion driving force of the mixed mode crack is obtained by decomposing the total elastic strain energy based on the decomposition of the symmetrical strain and the antisymmetric strain of the total strain and the tensile and compressive anisotropic decomposition of the elastic strain energy related to the symmetrical strain. By strict deduction, various fracture phase field models containing microstructure effect, such as fracture phase field model based on correction even stress theory, fracture phase field model based on cohesion theory, brittle fracture phase field model and the like, can be derived based on the method, so that a fracture phase field model with strong universality containing microstructure effect is established. As shown in fig. 1, the method comprises the steps of:
step one, considering the influence of micro bending freedom degree and micro torsion freedom degree in the microstructure, a micro pole elasticity theory is introduced to describe the elasticity mechanism of the material. In order to solve the problem of crack propagation in the mixed fracture mode, the total strain is decomposed into two parts of symmetrical strain and anti-symmetrical strain, and strain energy expressions corresponding to the symmetrical strain and the anti-symmetrical strain are established. The method comprises the following specific steps:
the elastic strain energy based on the micropolar elastic theory is expressed as:
wherein W is elastic strain energy density; c (C) ijkl 、D ijkl Elastic coefficient tensors corresponding to the strain and curvature tensors, respectively; epsilon ij As strain tensor, ε ij =u j,i -e ijk θ k ;u j,i Representing the right gradient of displacement, in particular meaning the "j" th displacement component versus the "i" th displacement componentTaking partial derivatives of the space coordinates; e, e ijk Is the third-order Levi-Civita tensor, θ k Is a micro-rotating field; kappa (kappa) ij As a curvature tensor, κ ij =θ j,i ,θ j,i The right gradient of the micro-rotation displacement is expressed, and the specific meaning is that the j-th micro-rotation displacement component takes partial derivative on the i-th space coordinate; i, j, k, l are tensor subscripts, taking integer values 1 to 3; elastic coefficient based on isotropic micro-polar elasticity theory is C ijkl =λδ ij δ kl +(μ+κ)δ ik δ jl +μδ il δ jk ,D ijkl =αδ ij δ kl +γδ ik δ jl +βδ il δ jk ,δ ij ,δ kl ,δ ik ,δ jl ,δ il ,δ jk Different index notations for the second order unit tensor; lambda, mu, kappa, alpha, beta, gamma are the mechanical parameters of the material. The conversion relation between each parameter and the basic mechanical parameter is shown as follows:
where G is the shear modulus, g=e/[ 2 (1+v)]The method comprises the steps of carrying out a first treatment on the surface of the V is poisson's ratio; e is the elastic modulus; chi is the polarizability and has no dimensionless constant coefficient; l (L) t Is the length scale of the torsion feature; l (L) b Is a curved feature length scale; n is the coupling number, and is a dimensionless constant coefficient.
Step two, obtaining the strain energy of the isotropic micro-polar elasticity theory based on the formula (1) as follows:
step one, three, strain tensor ε ij Can be decomposed into symmetrical and antisymmetric parts as follows:
wherein ,classical unpolarized theory as symmetric strain tensor; />Is an antisymmetric strain tensor.
And step four, decomposing strain energy based on the decomposition of strain in the step three, as follows:
wherein ,elastic strain energy related to symmetrical strain; />Is the elastic strain energy related to asymmetric strain and curvature.
Step two, considering the crack driving force of tensile and compressive anisotropy in the mixed fracture mode, carrying out tensile and compressive decomposition on elastic strain energy related to symmetrical strain, and further decomposing the elastic strain energy to obtain an expression of the elastic strain energy of each part; the thermodynamic basic framework based on the fracture phase field method and the decomposition formula of elastic strain energy establish the expression of the total potential energy functional. The method comprises the following specific steps:
step two, considering crack driving force of tensile and compressive anisotropies, carrying out tensile and compressive decomposition on elastic strain energy related to symmetrical strain, and decomposing strain energy corresponding to the symmetrical strain into the following components:
wherein ,the elastic strain energy corresponding to the tensile deformation is represented; />Elastic strain energy corresponding to compression deformation;<·> ± representing the braicket operator-> Is a symmetrical strain tensor in the principal strain space, "a" th principal strain value, n, representing symmetric strain a Representing a feature vector corresponding to the 'a' th main strain; tr (·) represents the apodization operator, +.>
Step two, based on the fracture phase field thermodynamic framework proposed by Miehe, borden et al, the total potential energy of the fracture phase field considering the microstructure effect is generally as follows:
wherein, the II is the total potential energy functional; psi is the total free energy density function;g (d) is an energy degradation function, and the classical energy degradation function, g (d) = (1-d), is selected according to the invention 20 ,κ 0 To ensure non-singular values of stiffness, kappa is generally taken 0 =10 -8 Other forms of energy degradation functions can be studied within the framework of the theory proposed by the present invention; alpha (d) is a geometric function, and the invention selects a classical geometric function alpha (d) =d 2 Other forms of geometric functions may be studied within the theoretical framework of the invention; c 0 Is a geometric coefficient corresponding to a classical geometric function, +.>G c Is the critical energy release rate; l (L) c Is a fracture length scale; />Representing the gradient of the fracture phase field; omega is the body domain; dV is the volume infinitesimal; f (f) i Is the volume stress; u (u) i Is a displacement field; pi i Is the volume couple stress; t (T) i Is the surface stress; Γ -shaped structure σ Is a plane stress distribution domain; dA is an area infinitesimal; m is M i Is the plane couple stress; Γ -shaped structure m Is a plane couple stress distribution domain.
Thirdly, deducing a control equation in a partial differential form from the total potential energy functional expression based on an energy transformation principle, wherein the control equation comprises a stress balance equation, an even stress balance equation and a fracture phase field driving force equation; for the irreversible evolution process of crack growth, the maximum value of the strain energy history up to the calculation time is adopted as the driving force of crack growth. The method comprises the following specific steps:
step three, deducing a corresponding control equation based on a fracture phase field energy variation principle:
for the displacement field part, a stress balance equation and corresponding boundary conditions are derived by a variation principle:
σ ji n j =T i ,X∈Γ σ
wherein δ is a variation symbol;representing an arbitrary value; x is the space coordinate of the object point; e represents belonging to; Γ -shaped structure u A displacement constraint distribution domain; sigma (sigma) ji Representing total stress->σ ji,j Representing the total stress sigma ji Taking partial derivative of the j-th space coordinate, wherein the physical meaning is the left gradient of the stress; n is n j Representing a constraint boundary outer normal vector; u (u) i The displacement is constrained for the boundary.
For the micro-rotating field, an even stress balance equation and corresponding boundary conditions are obtained based on a variation principle:
m ji n j =M i ,X∈Γ m
wherein ,Γθ Is a micro-rotation field constraint domain; m is m ji As even stress, m ji =g(d)·[a(κ kkij +γκ ji +bκ ij ];m ji,j Representing even stress m ji Partial derivative of the j-th space coordinate, the physical meaning is even stress left gradient;the micro-rotational displacement is constrained for boundaries.
For a fracture phase field, obtaining a fracture phase field driving force equation based on a variation principle:
where g '(d) represents the derivative of the energy degradation function with respect to the fracture phase field, g' (d) =d [ g (d) ]]/dd;Taking the historical maximum value of an elastic strain energy driving part for the irreversible evolution process of crack propagation as a historical fieldAs a history field, namely: />t represents the current time, τ represents a time variable; max { } is the maximum operator; a' (d) represents the derivative of the geometric function with respect to the fracture phase field; a' (d) =d [ a (d)]/dd; delta represents the laplace operator, delta= 2 ()/ 2 x i
The stress balance equation, the even stress balance equation and the fracture phase field driving force equation are fracture phase field partial differential form control equations containing microstructure effects.
Substituting the dual stress balance equation into a stress expression to obtain the following formula:
step III, multiplying e by two sides of the pair (3) mni Based onThe method comprises the following steps:
substituting the formula (4) into a stress balance equation to obtain the following components:
based on the formula (5), a fracture phase field method based on the modified even stress theory can be obtained. The method comprises the following steps: selecting an energy degradation function g (d) =d 2 Geometric functionCorrecting curvature tensor to symmetric tensor wherein θi,j The left gradient of the micro-rotation displacement is expressed, and the specific meaning is that the i-th micro-rotation displacement component takes partial derivative on the j-th space coordinate; the relation between the even stress and the curvature tensor is corrected to be m ij =g(d)(γκ ij ) (in planar problems, this formula is obtained based on both the modified even stress theory and the isotropic micropolar elasticity theory, but is different for spatial problems). Furthermore, a fracture phase field method based on the cohesive force model can be obtained based on the formula (5). The method comprises the following steps: the geometric function and the energy degradation function are respectively alpha (d) =d and g (d) = (1-d) 2 /[(1-d) 2 +m·d·(1+p·d)]At this time, geometric coefficientm is a dimensionless coefficient related to material parameters, m=g c /(c 0 ·ψ crit ·l c )≥1,ψ crit The critical strain energy of crack initiation is p, which is a model dimensionless parameter, and p is more than or equal to 1.
It can be seen that the method provided by the invention is general.
And fourthly, the control equation based on the partial differential form is equivalently converted into an integral equation of a weak form by means of the virtual work principle. The method comprises the following specific steps:
step four, one, the stress balance equation, the equal sign two sides are multiplied by the virtual displacement delta u i And integrating the body area to obtain:
Ωji,j +f i )δu i dV=0 (6);
wherein ,δui Representing the variation of the displacement.
Substituting the expression of stress into the formula (6) and deforming based on the Gaussian theorem to obtain:
wherein ,is related to displacement fieldIs an antisymmetric strain tensor of-> A variation representing a symmetrical strain; />A variation representing a displacement-related antisymmetric strain; u (u) j,i Representing the partial derivative of the "j" th displacement component with respect to the "i" th spatial coordinate; u (u) i,j Representing the partial derivative of the "i" th displacement component with respect to the "j" th spatial coordinate.
Step four, two, similar treatment, equal sign two sides of the dual stress balance equation are multiplied by micro rotation virtual displacement delta theta i And integrating the volume domain to obtain:
wherein ,δθi 、δθ j 、δθ k Different index notations representing the micro-rotational displacement variation; delta kappa ij Representing the variation of curvature.
And fourthly, similarly, multiplying virtual displacement delta d on two sides of a fracture phase field driving force equation and integrating a body domain to obtain the three-phase fracture phase field driving force equation:
wherein δd represents the variation of the fracture phase field; δd ,i Representing the variation of the fracture phase field gradient.
Equation (7), equation (8), equation (9) constitute a weak form of integral equation.
From the weak form of the integral equations (7) to (9)), the invention will degenerate into classical fracture phase field models without microstructural effects when the correlation coefficients α, β, γ, κ of the microstructure are all zero. This conclusion applies to both planar and three-dimensional models.
And fifthly, carrying out numerical value dispersion on each field variable (comprising a displacement field, a micro-rotation field and a fracture phase field) based on a finite element method, carrying out localized processing on a weak form integral equation at a cell, and deducing a cell stiffness equation. The method comprises the following specific steps:
step five, the present invention gives mainly a finite element discrete form of the planar problem, and the spatial problem can be deduced according to the following similar steps based on the formulas (7) to (9) in step four. Based on the finite element method, finite element interpolation expressions of a displacement field, a micro-rotation field and a fracture phase field are as follows:
wherein ,{ui The displacement column vector, { u }, is represented by i }=[u x u y ] TFor the "I" th node displacement, i=x, y,respectively representing the displacement of the 'I' node along the x-axis and the y-axis; i is an angle sign, and if no special description exists, the angle sign 'I' in each variable takes positive integer values of 1 to N; n is the node number of one unit; />Micro-rotational displacement of the 'I' th node around the z axis; d, d I Is the 'I' th nodeA fracture phase field; />For the displacement of the column vector for the node,N I interpolation shape function for the "I" th node; />For the node micro-rotation field column vector, +.>{d I The node fracture phase field column vector is { d } I }=[d 1 ... d N ] T ;N u A node displacement shape function matrix; n (N) θ A node micro-rotation displacement shape function matrix; n (N) d A matrix of phase field shape functions is broken for the nodes.
N θ =N d =[N 1 ...N i ...N N ]。
And step five, substituting the displacement, micro-rotation displacement and finite element interpolation expression of the fracture phase field in the step five into the weak form integral equation in the step four to respectively obtain a node force balance equation, a node even stress balance equation and a node fracture phase field driving force balance equation as follows:
node force balance equation:
node even stress balance equation:
node fracture phase field driving force balance equation: k (K) dd {d I }=F d
wherein ,Kuu A displacement dependent stiffness matrix; k (K) And K is equal to θu Is a displacement and micro-rotation displacement coupling stiffness matrix; k (K) θθ A micro-rotational displacement dependent stiffness matrix; k (K) dd Is a fracture phase field stiffness matrix; f (F) u Is the node force; f (F) θ Coupling stress for the node; f (F) d A fracture phase field driving force for the node; each physical quantity expression is as follows.
K θθ =∫ Ω g(d)·[γ(B θ ) T B θ +2κ·(N θ ) T N θ ]dV;
wherein ,a shape function gradient matrix for symmetrical strain; />A shape function gradient matrix for asymmetric strain; b (B) d Gradient matrix for fracture phase field shape function; b (B) θ A shape function gradient matrix for the micro-rotating field; { f } is the column vector of the volume stress, { f = [ f ] x f y ] T ,f x And f y Representing the volumetric stress along the x-axis and y-axis directions, respectively; { T } is the plane stress column vector, { T = [ T ] x T y ] T ,T x And T is y Respectively representing the plane stress along the x-axis and the y-axis; { pi } is the even stress column vector of the volume, { pi = [ pi ] xz π yz ] T ,π xz Representing the volumetric even stress in the x-direction in the plane x-z, pi yz Representing the volumetric even stress in the y-direction in plane y-z; { M } is the plane even stress column vector, { M = [ M = xz M yz ] T ,M xz Representing plane couple stress in the x-direction in the plane x-z, M yz Representing the plane even stress in the y-direction in plane y-z.
Step six, forming a nonlinear total stiffness equation (a highly nonlinear equation set) according to a cell assembly rule by using a cell stiffness equation; designing an iterative solution algorithm of the nonlinear total stiffness equation; the numerical analysis scheme can combine pre-processing software and post-processing software compatible with ABAQUS, such as finite element pre-processing software Hypermesh, post-processing software Hyperview and the like, so as to further popularize the proposed method to engineering application.
The unit stiffness equation in the fifth step is a highly nonlinear non-convex equation set, and the stable solving method of the highly nonlinear non-convex equation set is to use an alternate solving method to decompose the highly nonlinear non-convex equation set into two convex equation sets for solving. The basic idea of the alternating solution is shown in fig. 2 and 3. As shown in fig. 2, for each load increment step, the fracture phase field solved in the previous increment step is input into a node force balance equation and a node even stress balance equation as known values. At this time, since the fracture phase field is known, the displacement and micro-rotation field are solved and converted into a linear equation set; after the displacement field and the micro-rotation field are obtained, the current strain energy and the history field can be determined, then the fracture phase field of the current increment step is solved through Newton iteration based on the node fracture phase field driving force balance equation, and all field variables are updated. As shown in FIG. 3, each node has four degrees of freedom based on an alternate solution, including node displacement, node micro-rotation displacement and node fracture phase field parameters; the algorithm is provided with three layers of units in total, wherein the first layer of units only comprises a unit stiffness equation of a fracture phase field, and the overall stiffness equation of the fracture phase field can be formed based on the assembly rule of the units and the nodes, and is used for solving the fracture phase field. The second layer unit only comprises a unit stiffness equation of displacement and micro-rotation displacement, the displacement can be formed based on the assembly rule of the units and the nodes, and the micro-rotation displacement coupling total stiffness equation is used for solving the displacement and the micro-rotation displacement. The third layer unit was used for image display of all the solving field variables, and in order for the third layer unit to have no influence on the solving of the overall stiffness equation of each field variable, the stiffness of the layer unit was set to a minimum value (elastic modulus is 10 -11 Units).
The correctness and characteristics of the method provided by the invention will be described below in connection with three examples.
Calculation example one: the infinite microplates with the central circular aperture are subjected to remote tensile stress. The purpose of this calculation is to verify the correctness of the programmed program by means of an analytical solution. As shown in fig. 4, the infinite microplates with a central circular hole were analyzed for symmetry by taking only one quarter of the plate, with the dimensions of 5d'5d.This example is considered as a plane strain problem, with each material point having a degree of freedom u x ,u y θ and d. The material parameters are respectively as follows: elastic modulus e=130 GPa, poisson ratio v=0.3, bending length scale L b =d, α=β=0. The coupling number N is set to 0.1,0.5,0.9. The cells of this example used 8-node reduced integral (CPEQ 8R) to refine the cell size near the hole and relatively larger cell size far from the hole, with a total number of cells of 35500. The analytical solution referred to is the result of the analysis by Ariman et al. As can be seen from FIG. 5, the numerical solution of the microplates with the center circular holes is very consistent with the theoretical analytical solution, and the correctness of the programmed elastic portion is verified.
Calculating example II: prefabricated edge crack plate under stretching and shearing action. The present example is a verification of the degradation of the proposed fracture phase field method into a classical fracture phase field method. As shown in FIG. 6, the problem area is a square plate with a size of 1X 1mm 2 . The middle part of the plate is provided with a preset crack with the length of 0.5 mm. Tensile and shear load testing was performed by applying displacements of the top of the vertical plate and the top of the parallel plate, respectively, on top of the plates. The bottom of the plate was held stationary during the test. The meshing of the area along the path of possible crack propagation is relatively fine, with a cell size of about 0.002mm, about one quarter of the fracture length scale. In the whole analysis process, according to the existing research, the elastic modulus E=210 GPa, the Poisson ratio v=0.3 and the critical energy release rate G are selected according to the material parameters c =0.0075 GPa mm, breaking length scale l c =0.0075 mm, bending length scale l b =0.0 mm, the coupling number n=0.0. Referring to fig. 7 and 8, it can be known that the load displacement curve and crack propagation morphology under the action of the tensile and shear displacement loads obtained by the method of the present invention are highly consistent with the analysis results of the classical fracture phase field method, which proves that the method of the present invention can be degraded into the fracture phase field method based on the linear elastic theory. Calculating example III: this example is used to illustrate the microstructure effects of a microstructure-containing material. The material parameters are selected as elastic modulus E=210 GPa, poisson ratio v=0.3 and critical energy release rate G c =0.0075 GPa mm, breaking length scale l c =0.0075 mm, bending length scale l b The coupling numbers N were 0.0,0.2,0.4,0.6,0.8 and 0.5mm, respectively. The mesh size of the refinement region was chosen to be h=0.002 mm. As can be seen from fig. 9 and 10, for the stretched prefabricated edge crack plate, the coupling numbers take different values and have no influence on the crack growth morphology and the load displacement curve; for a sheared edge crack plate, as the coupling number increases, the elastic segment stiffness increases, the load peak increases, the critical deformation value for material failure decreases, the crack initiation angle (the angle between the crack propagation direction and the horizontal direction) of the crack decreases, and the propagation morphology of the crack tends to be more toward the horizontal direction (fig. 10). This is because crack propagation of microstructure-containing materials can be divided into two phases: firstly, micro-rotation in the opposite direction is generated by the micro-structure at the tip of the crack, secondly, the crack is further expanded to cut off the connection of the micro-structure, and compared with the material without the micro-structure effect, the elastic energy stored by the connection of the micro-structure is released, so that the crack expansion is driven more easily, and the critical deformation value of the material failure is reduced. And based on the isotropic micro-polar elasticity theory, the coupling between the stretching deformation and the micro-rotation deformation does not exist. In contrast, there is a strong coupling effect between shear deformation and micro-rotational deformation. The increase in the coupling number directly reflects the enhancement of the coupling effect of the rotational deformation and the shear deformation of the microstructure. The above simulation results appear by combining the reasons of the above aspects.

Claims (10)

1. A method of fracture phase field simulation comprising a microstructure effect, the method comprising the steps of:
step one, considering the influence of micro bending freedom and micro torsion freedom in a microstructure, introducing a micro-polar elasticity theory to describe the elasticity mechanism of the material; in order to solve the problem of crack propagation of a mixed fracture mode, decomposing total strain into two parts of symmetrical strain and anti-symmetrical strain, and establishing an elastic strain energy expression corresponding to the symmetrical strain and the anti-symmetrical strain;
step two, considering the crack driving force of tensile and compressive anisotropy in the mixed fracture mode, further carrying out tensile and compressive decomposition on the elastic strain energy related to the symmetrical strain, and obtaining the expression of the elastic strain energy of each part; establishing an expression of a total potential energy functional based on a thermodynamic basic frame of a fracture phase field method and decomposition of elastic strain energy;
thirdly, deducing a control equation in a partial differential form from the total potential energy functional expression based on an energy transformation principle, wherein the control equation comprises a stress balance equation, an even stress balance equation and a fracture phase field driving force equation; aiming at the irreversible evolution process of crack growth, adopting the maximum value of strain energy history until the calculation time as the driving force of crack growth;
step four, a control equation based on a partial differential form is equivalently converted into an integral equation of a weak form by means of a virtual work principle;
fifthly, carrying out numerical discrete on a displacement field, a micro-rotation field and a fracture phase field based on a finite element method, carrying out localized treatment on a weak form integral equation at a cell, and deducing a cell stiffness equation;
and step six, forming a nonlinear total stiffness equation by the cell stiffness equation according to a cell assembly rule, and designing an iterative solution algorithm of the nonlinear total stiffness equation.
2. The method for simulating a fracture phase field containing a microstructure effect according to claim 1, wherein the specific steps of the step one are as follows:
the elastic strain energy based on the micropolar elastic theory is expressed as:
wherein W is elastic strain energy density; c (C) ijkl 、D ijkl Elastic coefficient tensors corresponding to the strain and curvature tensors, respectively; epsilon ij Is the strain tensor; kappa (kappa) ij For curvature tensors i, j, k and l, tensor subscripts are taken as integer values 1 to 3;
step two, obtaining the strain energy of the isotropic micro-polar elasticity theory based on the formula (1) as follows:
step one, three, strain tensor ε ij Can be decomposed into symmetrical and antisymmetric parts as follows:
wherein ,is a symmetrical strain tensor; />Is an antisymmetric strain tensor;
and step four, decomposing strain energy based on the decomposition of strain in the step three, as follows:
wherein ,elastic strain energy related to symmetrical strain; />Is the elastic strain energy related to asymmetric strain and curvature.
3. The method of claim 1, wherein the C is a phase field simulation method ijkl =λδ ij δ kl +(μ+κ)δ ik δ jl +μδ il δ jk ,D ijkl =αδ ij δ kl +γδ ik δ jl +βδ il δ jk ,δ ij ,δ kl ,δ ik ,δ jl ,δ il ,δ jk Different index notations for the second order unit tensor; lambda, mu, kappa, alpha, beta and gamma are mechanical parameters of the material, and the conversion relation between each parameter and the basic mechanical parameter is shown as the following formula:
wherein G is the shear modulus; v is poisson's ratio; e is the elastic modulus; chi is the polarizability and has no dimensionless constant coefficient; l (L) t Is the length scale of the torsion feature; l (L) b Is a curved feature length scale; n is the coupling number, and is a dimensionless constant coefficient.
4. The method for simulating a fracture phase field containing a microstructure effect according to claim 1, wherein the specific steps of the second step are as follows:
step two, considering crack driving force of tensile and compressive anisotropies, carrying out tensile and compressive decomposition on elastic strain energy related to symmetrical strain, and decomposing strain energy corresponding to the symmetrical strain into the following components:
wherein ,the elastic strain energy corresponding to the tensile deformation is represented; />Elastic strain energy corresponding to compression deformation;<·> ± representing a braicket operator; />Is a symmetrical strain tensor under the main strain space;
based on a fracture phase field thermodynamic frame, the total potential energy form of the fracture phase field considering the microstructure effect is as follows:
wherein, pi is the total potential energy functional; psi is the total free energy density function; g (d) is an energy degradation function; α (d) is a geometric function; c 0 Is a geometric coefficient; g c Is the critical energy release rate; l (L) c Is a fracture length scale;representing the gradient of the fracture phase field; omega is the body domain; dV is the volume infinitesimal; f (f) i Is the volume stress; u (u) i Is a displacement field; pi i Is the volume couple stress; t (T) i Is the surface stress; Γ -shaped structure σ Is a plane stress distribution domain; dA is an area infinitesimal; m is M i Is the plane couple stress; Γ -shaped structure m Is a plane couple stress distribution domain.
5. The fracture phase field simulation method comprising the microstructure effect according to claim 1, wherein the specific steps of the third step are as follows:
step three, deducing a corresponding control equation based on a fracture phase field energy variation principle:
for the displacement field part, a stress balance equation and corresponding boundary conditions are derived by a variation principle:
σ ji n j =T i ,X∈Γ σ
wherein δ is a variation symbol;representing an arbitrary value; x is the space coordinate of the object point; e tableShowing the belongings; Γ -shaped structure u A displacement constraint distribution domain; sigma (sigma) ji Representing the total stress; n is n j Representing a constraint boundary outer normal vector; />Constraining the displacement for the boundary;
for the micro-rotating field, an even stress balance equation and corresponding boundary conditions are obtained based on a variation principle:
m ji n j =M i ,X∈Γ m
wherein ,Γθ Is a micro-rotation field constraint domain; m is m ji Is even stress; m is m ji,j Representing even stress m ji Partial derivatives to the "j" th spatial coordinate;constraint of micro-rotational displacement for boundaries;
for a fracture phase field, obtaining a fracture phase field driving force equation based on a variation principle:
where g' (d) represents the derivative of the energy degradation function with respect to the fracture phase field;taking a historical maximum value of an elastic strain energy driving part as a historical field aiming at an irreversible evolution process of crack propagation; alpha' (d) represents the derivative of the geometric function with respect to the fracture phase fieldThe method comprises the steps of carrying out a first treatment on the surface of the Delta represents the laplace operator;
substituting the dual stress balance equation into a stress expression to obtain the following formula:
step III, multiplying e by two sides of the pair (3) mni Based on e kli e mni =δ km δ lnkn δ lm A kind of electronic device with high-pressure air-conditioning systemThe method comprises the following steps:
substituting the formula (4) into a stress balance equation to obtain the following components:
6. the fracture phase field simulation method comprising the microstructure effect according to claim 1, wherein the specific steps of the fourth step are as follows:
step four, one, the stress balance equation, the equal sign two sides are multiplied by the virtual displacement delta u i And integrating the body area to obtain:
Ωji,j +f i )δu i dV=0 (6);
wherein ,δui A variation representing the displacement;
substituting the expression of stress into the formula (6) and deforming based on the Gaussian theorem to obtain:
wherein ,is an antisymmetric strain tensor related to the displacement field; />A variation representing a symmetrical strain; />A variation representing a displacement-related antisymmetric strain;
step four, two, equal sign two sides of dual stress balance equation are multiplied by micro rotation virtual displacement delta theta i And integrating the volume domain to obtain:
wherein ,δθi 、δθ j 、δθ k Different index notations representing the micro-rotational displacement variation; delta kappa ij Representing a variation of curvature;
and step four, multiplying virtual displacement delta d on two sides of a fracture phase field driving force equation and integrating a body domain to obtain the three-phase fracture phase field driving force equation:
wherein δd represents the variation of the fracture phase field; δd ,i A variation representing a fracture phase field gradient;
equation (7), equation (8), equation (9) constitute a weak form of integral equation.
7. The fracture phase field simulation method comprising the microstructure effect according to claim 1, wherein the specific steps of the fifth step are as follows:
step five, for the finite element discrete form of the plane problem, based on the finite element method, the finite element interpolation expressions of the displacement field, the micro-rotation field and the fracture phase field are as follows:
wherein ,{ui -representing a displacement column vector;for the "I" th node displacement, i=x, y; i is an angle sign, and the angle sign 'I' in each variable takes positive integer values of 1 to N; n is the node number of one unit; />Micro-rotational displacement of the 'I' th node around the z axis; d, d I Breaking the phase field for the "I" th node; />A column vector is displaced for the node; n (N) I Interpolation shape function for the "I" th node; />A micro-rotation field column vector is used for the node; { d I And the phase field column vector is the node fracture; n (N) u A node displacement shape function matrix; n (N) θ A node micro-rotation displacement shape function matrix; n (N) d A node fracture phase field shape function matrix;
and step five, substituting the displacement, micro-rotation displacement and finite element interpolation expression of the fracture phase field in the step five into the weak form integral equation in the step four to respectively obtain a node force balance equation, a node even stress balance equation and a node fracture phase field driving force balance equation as follows:
node force balance equation:
node even stress balance equation:
node fracture phase field driving force balance equation: k (K) dd {d I }=F d
wherein ,Kuu A displacement dependent stiffness matrix; k (K) And K is equal to θu Is a displacement and micro-rotation displacement coupling stiffness matrix; k (K) θθ A micro-rotational displacement dependent stiffness matrix; k (K) dd Is a fracture phase field stiffness matrix; f (F) u Is the node force; f (F) θ Coupling stress for the node; f (F) d Is the fracture phase field driving force of the node.
8. The method for fault phase field simulation including a microstructure effect as claimed in claim 7, wherein the K is uu 、K 、K θu 、K θθ 、K dd 、F u 、F θ 、F d The expression is as follows:
K θθ =∫ Ω g(d)·[γ(B θ ) T B θ +2κ·(N θ ) T N θ ]dV;
wherein ,a shape function gradient matrix for symmetrical strain; />A shape function gradient matrix for asymmetric strain; b (B) d Gradient matrix for fracture phase field shape function; b (B) θ A shape function gradient matrix for the micro-rotating field; { f } is the volume stress column vector; { T } is the plane stress column vector; { pi } is the volume even stress column vector; { M } is the plane even stress column vector.
9. The method of claim 7, wherein in the fifth step, for the spatial problem, the corresponding finite element discrete form can be derived in the plane problem step based on the integration equation of the weak form in the fourth step.
10. The method of claim 1, wherein in the sixth step, the iterative solution algorithm of the nonlinear total stiffness equation uses an alternate solution method.
CN202310441617.6A 2023-04-23 2023-04-23 Fracture phase field simulation method containing microstructure effect Active CN116486953B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202310441617.6A CN116486953B (en) 2023-04-23 2023-04-23 Fracture phase field simulation method containing microstructure effect

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202310441617.6A CN116486953B (en) 2023-04-23 2023-04-23 Fracture phase field simulation method containing microstructure effect

Publications (2)

Publication Number Publication Date
CN116486953A CN116486953A (en) 2023-07-25
CN116486953B true CN116486953B (en) 2023-09-05

Family

ID=87215114

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202310441617.6A Active CN116486953B (en) 2023-04-23 2023-04-23 Fracture phase field simulation method containing microstructure effect

Country Status (1)

Country Link
CN (1) CN116486953B (en)

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2005060631A2 (en) * 2003-12-11 2005-07-07 Ohio University Titanium alloy microstructural refinement method and high temperature, high strain rate superplastic forming of titanium alloys
WO2021248850A1 (en) * 2020-06-11 2021-12-16 大连理工大学 Method for predicting structural damage by using strength criterion-driven near-field dynamic model
CN113820474A (en) * 2021-11-24 2021-12-21 中南大学 Phase field method for simulating composite crack propagation of brittle rock
CN114255825A (en) * 2021-12-22 2022-03-29 中北大学 Method for predicting influence of strain loading on microcrack expansion based on crystal phase field method
CN115410663A (en) * 2022-08-16 2022-11-29 大连理工大学 Dynamic impact/contact elastoplasticity large deformation fracture analysis explicit phase field material point method

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2005060631A2 (en) * 2003-12-11 2005-07-07 Ohio University Titanium alloy microstructural refinement method and high temperature, high strain rate superplastic forming of titanium alloys
WO2021248850A1 (en) * 2020-06-11 2021-12-16 大连理工大学 Method for predicting structural damage by using strength criterion-driven near-field dynamic model
CN113820474A (en) * 2021-11-24 2021-12-21 中南大学 Phase field method for simulating composite crack propagation of brittle rock
CN114255825A (en) * 2021-12-22 2022-03-29 中北大学 Method for predicting influence of strain loading on microcrack expansion based on crystal phase field method
CN115410663A (en) * 2022-08-16 2022-11-29 大连理工大学 Dynamic impact/contact elastoplasticity large deformation fracture analysis explicit phase field material point method

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
晶体相场模型及其在材料微结构演化中的应用;高英俊;卢昱江;孔令一;邓芊芊;黄礼琳;罗志荣;;金属学报(第02期);全文 *

Also Published As

Publication number Publication date
CN116486953A (en) 2023-07-25

Similar Documents

Publication Publication Date Title
Placidi et al. Simulation results for damage with evolving microstructure and growing strain gradient moduli
Fish et al. Computational damage mechanics for composite materials based on mathematical homogenization
Amendola et al. Experimental response of additively manufactured metallic pentamode materials confined between stiffening plates
Teotia et al. Applications of finite element modelling in failure analysis of laminated glass composites: A review
Ye et al. Failure analysis of fiber-reinforced composites subjected to coupled thermo-mechanical loading
Fish et al. Computational mechanics of fatigue and life predictions for composite materials and structures
Thamburaja Length scale effects on the shear localization process in metallic glasses: a theoretical and computational study
Zhang et al. Phenomenological crystal plasticity modeling and detailed micromechanical investigations of pure magnesium
Forest et al. Estimating the overall properties of heterogeneous Cosserat materials
US11193267B2 (en) Tensegrity structures and methods of constructing tensegrity structures
Chen et al. A nonlocal lattice particle model for fracture simulation of anisotropic materials
Ren et al. Modeling and simulation of large-scale ductile fracture in plates and shells
Mbiakop et al. On void shape effects of periodic elasto-plastic materials subjected to cyclic loading
Brünig et al. Nonlocal continuum theory of anisotropically damaged metals
Miller et al. Modeling large strain multiaxial effects in FCC polycrystals
Rousselier et al. A novel approach for anisotropic hardening modeling. Part I: Theory and its application to finite element analysis of deep drawing
Haj-Ali et al. Nonlinear micromechanical formulation of the high fidelity generalized method of cells
Naumenko et al. Damage patterns in float glass plates: Experiments and peridynamics analysis
Germain et al. Composite layered materials: Anisotropic nonlocal damage models
Wanji et al. A study of scale effect of composite laminated plates based on new modified couple stress theory by finite-element method
Tjahjanto et al. A novel grain cluster-based homogenization scheme
CN116486953B (en) Fracture phase field simulation method containing microstructure effect
Hommel et al. A hybrid modeling concept for ultra low cycle fatigue of metallic structures based on micropore damage and unit cell models
Muliana et al. Nested nonlinear viscoelastic and micromechanical models for the analysis of pultruded composite materials and structures
Joshan et al. A couple stress model in non-polynomial framework to examine structural responses of laminated composite micro-plates: An analytical solution

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant