CN112861305A - Crack propagation direction prediction method and device and storage medium - Google Patents

Crack propagation direction prediction method and device and storage medium Download PDF

Info

Publication number
CN112861305A
CN112861305A CN202011523316.0A CN202011523316A CN112861305A CN 112861305 A CN112861305 A CN 112861305A CN 202011523316 A CN202011523316 A CN 202011523316A CN 112861305 A CN112861305 A CN 112861305A
Authority
CN
China
Prior art keywords
time
displacement
strain
stress
initial
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.)
Pending
Application number
CN202011523316.0A
Other languages
Chinese (zh)
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.)
Tsinghua University
Institute of Flexible Electronics Technology of THU Zhejiang
Original Assignee
Tsinghua University
Institute of Flexible Electronics Technology of THU Zhejiang
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 Tsinghua University, Institute of Flexible Electronics Technology of THU Zhejiang filed Critical Tsinghua University
Priority to CN202011523316.0A priority Critical patent/CN112861305A/en
Publication of CN112861305A publication Critical patent/CN112861305A/en
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/15—Correlation function computation including computation of convolution operations
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Data Mining & Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Pure & Applied Mathematics (AREA)
  • Computational Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • Software Systems (AREA)
  • Databases & Information Systems (AREA)
  • Algebra (AREA)
  • Computing Systems (AREA)
  • Operations Research (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • Investigating Strength Of Materials By Application Of Mechanical Stress (AREA)

Abstract

The application relates to a method, a device and a storage medium for predicting crack propagation direction, wherein the prediction method comprises the following steps: calculating the displacement field of the material according to the strain-rotation and decomposition theorem and a gridless method; calculating a stress-strain field according to the calculation result of the displacement field; and calculating the crack propagation direction according to the calculation result of the stress-strain field and by combining a crack propagation criterion. On the basis of a gridless method, the method provides the strain-rotation and decomposition theorem, forms a new crack propagation direction prediction method, can improve the accuracy of numerical analysis of the nonlinear deformation fracture problem, and provides new ideas and means for the numerical analysis of the fracture problem.

Description

Crack propagation direction prediction method and device and storage medium
Technical Field
The application relates to the technical field of fracture problem numerical analysis, in particular to a crack propagation direction prediction method and device and a storage medium.
Background
The problem of fracturing has long been one of the most important and complex problems in various engineering fields. Fracture analysis becomes particularly important, especially with the rapid development of advanced materials and structures. For the fracture problem, the analytical solution can be found only in a very simple case, and the approximate calculation of the fracture problem can be usually realized only by a numerical method.
Finite element is one of the numerical calculation methods commonly used in engineering calculation, and the finite element based expansion finite element method (XFEM) is a very common method for dealing with the crack problem at present. However, the finite element method needs to be processed by means of grid cells, and has the problems of cell distortion, complex preprocessing process, time consumption, large calculation amount and the like, so that the problem of complex crack propagation caused by the computational efficiency and accuracy of the discontinuous problem (such as the crack propagation problem) is greatly reduced; in addition, because the expansion path and the original unit interface can not be matched, the finite element is difficult to realize accurate calculation, and if the large deformation of materials and structures is accompanied before and after the fracture, the distortion of the unit can seriously influence the solving precision of the finite element.
The emerging meshless-Free Methods in the last 30 years can partially or completely eliminate meshes and directly utilize discrete nodes to analyze problems, and the unique characteristic can overcome partial problems existing in finite elements and ensure the precision of results, so the meshless method has great advantages for processing crack propagation problems. However, all the fracture problem numerical methods based on the meshless method are based on the small deformation theory basis at present, and although the linear elastic fracture mechanics is concise and perfect, the linear elastic fracture mechanics is strictly not suitable for fracture cusp field analysis due to the limitation of the linear elastic theory in the aspect of describing larger deformation, so that the accuracy of the numerical analysis of the nonlinear deformation fracture problem is limited from the theory itself.
Disclosure of Invention
In view of the above technical problems, the present application provides a method, an apparatus, and a storage medium for predicting a crack propagation direction, so as to improve accuracy of numerical analysis of a fracture problem due to nonlinear deformation, and provide a new idea and means for numerical analysis of the fracture problem.
In order to solve the above technical problem, the present application provides a method for predicting a crack propagation direction, including: calculating the displacement field of the material according to the strain-rotation and decomposition theorem and a gridless method; calculating a stress-strain field according to the calculation result of the displacement field; and calculating the crack propagation direction according to the calculation result of the stress-strain field and by combining a crack propagation criterion.
In one embodiment, the strain-rotation and resolution theorem includes:
F＝S+R
Figure BDA0002849582520000021
Figure BDA0002849582520000022
wherein,
Figure BDA0002849582520000023
Figure BDA0002849582520000024
wherein F is a reversible linear differential transformation function; s, S,
Figure BDA0002849582520000025
Is the strain tensor; r, R,
Figure BDA0002849582520000026
Is the tensor of rotation; u. ofi|jIs a displacement component uiTo coordinate xjA derivative of (a);
Figure BDA0002849582520000027
is ui|jThe transposed matrix of (2);
Figure BDA0002849582520000028
is a kronecker function;
Figure BDA0002849582520000029
is a unit vector of the rotating shaft; theta is the average angle of full rotation.
In one embodiment, the step of calculating the displacement field of the material according to the strain-rotation and decomposition theorem and the gridless method comprises:
obtaining an incremental variational equation according to the strain-rotation and decomposition theorem:
Figure BDA00028495825200000210
wherein the prime symbol "-" represents the physical quantity under the initial tow at time t in the tow coordinate system,
Figure BDA00028495825200000211
is a stress tensor,
Figure BDA00028495825200000212
Is an elastic matrix; Δ t is the time increment;
Figure BDA00028495825200000213
is the strain rate;
Figure BDA00028495825200000214
is the rate of rotation;t+ΔtWineandt+ΔtWextrespectively is inertia force power and external force virtual power; delta is a variation sign; Ω is the area of the material;
according to the gridless method, a moving least square interpolation function is adopted to obtain the displacement increment and the speed at the time t:
Figure BDA00028495825200000215
Figure BDA00028495825200000216
wherein,tψk(x) A moving least squares function for node k at time t;tΔukandtVkrespectively the displacement increment and the speed of the node k at the time t;tΔ u is the displacement increment at the time t,tAnd V is the speed at the time t.
In one embodiment, the step of calculating the displacement field of the material according to the strain-rotation and decomposition theorem and the gridless method comprises:
substituting the displacement increment of the time t and the speed of the time t into the increment variational equation to obtain a system discrete equation:
t+ΔtKe tΔu＝t+ΔtFe
wherein,
Figure BDA0002849582520000031
Figure BDA0002849582520000032
wherein,t+ΔtKe、t+ΔtFerespectively an equivalent stiffness matrix and an equivalent force vector at the moment of t + delta t;tKL、tKN、tM、tf is a linear stiffness matrix, a nonlinear stiffness matrix, a mass matrix and a stress integral vector at the moment t respectively; kαFor the penalty function stiffness matrix, FαIs a penalty function force vector;t+Δtr is an external force vector at the moment of t + delta t;ta is the acceleration at the time t; the parameter β is determined according to the accuracy and stability requirements of the integration.
In one embodiment, the step of calculating the displacement field of the material according to the strain-rotation and decomposition theorem and the gridless method further comprises:
according to the Newmark method, the relationship among the displacement increment, the speed and the acceleration at different moments is obtained:
Figure BDA0002849582520000033
wherein,t+Δta is the acceleration at time t + Deltat,tA is the acceleration at time t,t+ΔtV is the speed at time t + Deltat,tV is the speed at time t,tAnd delta u is the displacement increment at the moment t, delta t is the time increment, and the parameters beta and gamma are determined according to the requirements of the integral precision and stability.
In one embodiment, the step of calculating the displacement field of the material according to the strain-rotation and decomposition theorem and the gridless method comprises:
obtaining an initial equivalent stiffness matrix, an initial mass matrix, an initial equivalent force vector and an initial displacement, and calculating an initial acceleration according to the following formula:
A＝M-1(Fe-Keu0)
wherein A is the initial acceleration, KeIs the initial equivalent stiffness matrix, M is the initial mass matrix, FeIs the initial equivalent force vector, u0Is the initial displacement.
In one embodiment, the step of calculating the displacement field of the material according to the strain-rotation and decomposition theorem and the gridless method further comprises:
obtaining a correction increment variational equation according to the strain-rotation and decomposition theorem:
Figure BDA0002849582520000041
wherein the prime symbol "-" represents the physical quantity under the initial tow at time t in the tow coordinate system,
Figure BDA0002849582520000042
is a stress tensor,
Figure BDA0002849582520000043
Is an elastic matrix; Δ t is the time increment;
Figure BDA0002849582520000044
is the strain rate;
Figure BDA0002849582520000045
a rate of rotation;t+Δ tWineandt+ΔtWextrespectively is inertia force power and external force virtual power; delta is a variation sign; Ω is the area of the material;
Figure BDA0002849582520000046
is a penalty factor; s is a boundary; vi、VjIs the speed;
Figure BDA0002849582520000047
is the velocity determined by the boundary conditions.
In one embodiment, the step of calculating a stress-strain field from the calculation of the displacement field includes:
and acquiring a stress tensor according to the strain tensor by adopting the following formula:
Figure BDA0002849582520000048
wherein,
Figure BDA00028495825200000415
is the stress tensor;
Figure BDA0002849582520000049
is an elastic matrix;
Figure BDA00028495825200000410
is the strain tensor.
In one embodiment, the step of calculating the crack propagation direction according to the calculation result of the stress-strain field and the crack propagation criterion includes:
and obtaining the main stress and the main direction thereof by adopting the following formula according to the weighted stress of the crack tip:
Figure BDA00028495825200000411
wherein σ1＞σ2，
Figure BDA00028495825200000412
Wherein σtipA weighted stress of the crack tip; sigma1、σ2Is the principal stress, n1，n2Is the primary direction.
In one embodiment, the step of calculating the crack propagation direction according to the calculation result of the stress field and the crack propagation criterion includes:
the weighted stress of the crack tip is calculated according to the following formula:
Figure BDA00028495825200000413
wherein,
Figure BDA00028495825200000414
wherein σ1Selecting the stress of a node near the tip of the crack; n isqThe number of integration points near the crack tip; w is aIIs the weight of the integral point; dmaxIs the crack tip andmaximum value of distance between selected nodes in the vicinity thereof, dIFor the corresponding integration point xiA distance from the tip of the crack, and dIAnd dmaxSatisfy the relationship 0 < dI≤2dmax。
In one embodiment, before the step of calculating the crack propagation direction according to the calculation result of the stress-strain field and the crack propagation criterion, the method comprises: judging whether the crack length reaches the maximum value; judging whether the position of the crack tip reaches the boundary of the problem domain of the material; stopping the numerical calculation process if the crack length reaches the maximum value and/or the crack tip position reaches the problem domain boundary; and if the crack length does not reach the maximum value and the crack tip position does not reach the problem domain boundary, executing the step of calculating the crack propagation direction according to the calculation result of the stress-strain field and by combining the crack propagation criterion.
In one embodiment, the prediction method includes the following steps: acquiring an initial characteristic matrix, initial displacement and initial speed, and obtaining initial acceleration according to the initial characteristic matrix and the initial displacement; selecting a time step length and calculating an integral constant; acquiring a stiffness characteristic matrix at the time t, and acquiring an equivalent stiffness matrix at the time t + delta t according to the stiffness characteristic matrix at the time t and the integral constant; acquiring an equivalent force characteristic matrix at the time t, the speed at the time t and the acceleration at the time t, and acquiring an equivalent force vector at the time t + delta t according to the equivalent force characteristic matrix at the time t, the speed at the time t, the acceleration at the time t and the integral constant; obtaining a displacement increment at the t moment according to the equivalent stiffness matrix at the t + delta t moment and the equivalent force vector at the t + delta t moment; acquiring the displacement at the time t, and acquiring the displacement at the time t + delta t according to the displacement at the time t and the displacement increment at the time t; obtaining the speed at the t + delta t moment and the acceleration at the t + delta t moment according to the displacement increment at the t moment, the speed at the t moment and the acceleration at the t moment; and (5) circulating the time step, wherein t is t + delta t, and returning to execute the step three.
The application also provides a device for predicting the crack propagation direction, which comprises a memory, a processor and a computer program stored in the memory and capable of running on the processor, wherein the processor realizes the steps of the prediction method when executing the computer program.
The present application also provides a computer-readable storage medium, in which a computer program is stored, which computer program, when being executed by a processor, realizes the steps of the above-mentioned prediction method.
The method provides a strain-rotation and decomposition theorem on the basis of a gridless method, forms a new crack propagation direction prediction method, can effectively improve the accuracy of numerical analysis of the nonlinear deformation fracture problem, and provides new ideas and means for the numerical analysis of the fracture problem.
Drawings
FIG. 1 is a schematic flow chart of a prediction method according to an embodiment of the present disclosure;
FIG. 2 is a schematic diagram of a tow coordinate system provided in an embodiment of the present application;
FIG. 3 is a schematic view of a crack tip and its corresponding node according to an embodiment of the present disclosure;
fig. 4 is a schematic flowchart of a prediction method according to a second embodiment of the present application;
fig. 5 is a schematic flowchart of a prediction method provided in the third embodiment of the present application;
fig. 6 is a schematic structural diagram of a prediction apparatus according to a fourth embodiment of the present application.
Detailed Description
The following description of the embodiments of the present application is provided for illustrative purposes, and other advantages and capabilities of the present application will become apparent to those skilled in the art from the present disclosure.
In the following description, reference is made to the accompanying drawings that describe several embodiments of the application. It is to be understood that other embodiments may be utilized and that mechanical, structural, electrical, and operational changes may be made without departing from the spirit and scope of the present application. The following detailed description is not to be taken in a limiting sense, and the scope of embodiments of the present application is defined only by the claims of the issued patent. The terminology used herein is for the purpose of describing particular embodiments only and is not intended to be limiting of the application.
Although the terms first, second, etc. may be used herein to describe various elements in some instances, these elements should not be limited by these terms. These terms are only used to distinguish one element from another.
Also, as used herein, the singular forms "a", "an" and "the" are intended to include the plural forms as well, unless the context indicates otherwise. It will be further understood that the terms "comprises," "comprising," "includes" and/or "including," when used in this specification, specify the presence of stated features, steps, operations, elements, components, items, species, and/or groups, but do not preclude the presence, or addition of one or more other features, steps, operations, elements, components, species, and/or groups thereof. The terms "or" and/or "as used herein are to be construed as inclusive or meaning any one or any combination. Thus, "A, B or C" or "A, B and/or C" means "any of the following: a; b; c; a and B; a and C; b and C; A. b and C ". An exception to this definition will occur only when a combination of elements, functions, steps or operations are inherently mutually exclusive in some way.
Fig. 1 is a schematic flowchart of a prediction method according to an embodiment of the present application. As shown in fig. 1, the prediction method provided by the present application includes the following steps:
step S101: calculating the displacement field of the material according to the strain-rotation and decomposition theorem and a gridless method;
in one embodiment, the strain-rotation and resolution theorem includes:
F＝S+R (1)
Figure BDA0002849582520000061
Figure BDA0002849582520000062
wherein,
Figure BDA0002849582520000063
Figure BDA0002849582520000071
wherein F is a reversible linear differential transformation function; s, S,
Figure BDA0002849582520000072
Is the strain tensor; r, R,
Figure BDA0002849582520000073
Is the tensor of rotation; u. ofi|jIs a displacement component uiTo coordinate xjA derivative of (a);
Figure BDA0002849582520000074
is ui|jThe transposed matrix of (2);
Figure BDA00028495825200000713
is a kronecker function; theta is the average angle of full rotation.
Taking a two-dimensional problem as an example, in order to describe the motion of any deformable body in space, two reference frames can be selected, as shown in fig. 2: (1) global fixed reference system { Xi(i ═ 1, 2); (2) natural towing coordinate system { x embedded in deformable bodyiAnd (i is 1,2), the coordinate line of which deforms along with the deformation of the object. The motion deformation process of the object is completely carried out by an initial towing coordinate system
Figure BDA0002849582520000075
And instantaneous tow coordinate system (g)i) Are defined and derived. The strain-rotation and decomposition theorem states that a unique direct and decomposition exists for any one reversible linear differential transformation function F. Wherein r, tr、t+ΔtrIs relative to the O pointA position vector of (a); u is t0The displacement to the time t, Δ u, is the displacement increment from the time t to the time t + Δ t.
In one embodiment, the step of calculating the displacement field of the material according to the strain-rotation and decomposition theorem and the gridless method comprises:
1) obtaining an incremental variational equation according to the strain-rotation and decomposition theorem:
Figure BDA0002849582520000076
wherein the prime symbol "-" represents the physical quantity under the initial tow at time t in the tow coordinate system,
Figure BDA0002849582520000077
is a stress tensor,
Figure BDA0002849582520000078
Is an elastic matrix; Δ t is the time increment;
Figure BDA0002849582520000079
is the strain rate;
Figure BDA00028495825200000710
is the rate of rotation;t+ΔtWineandt+ΔtWextrespectively is inertia force power and external force virtual power; delta is a variation sign; Ω is the area of the material;
2) according to the gridless method, a moving least square interpolation shape function is adopted for approximation, and displacement increment and speed at the time t are obtained:
Figure BDA00028495825200000711
Figure BDA00028495825200000712
wherein,tψk(x) A moving least squares function for node k at time t;tΔukandtVkrespectively the displacement increment and the speed of the node k at the time t;tΔ u is the displacement increment at time t,tV is the velocity at time t.
In one embodiment, the step of calculating the displacement field of the material according to the strain-rotation and decomposition theorem and the gridless method comprises:
substituting equations (7) and (8) into equation (6) to obtain a system discrete equation:
t+ΔtKe tΔu＝t+ΔtFe (9)
wherein,
Figure BDA0002849582520000081
Figure BDA0002849582520000082
wherein,t+ΔtKe、t+ΔtFerespectively an equivalent stiffness matrix and an equivalent force vector at the moment of t + delta t;tKL、tKN、tM、tf is a linear stiffness matrix, a nonlinear stiffness matrix, a mass matrix and a stress integral vector at the moment t respectively; kαFor the penalty function stiffness matrix, FαIs a penalty function force vector;t+Δtr is an external force vector at the moment of t + delta t;ta is the acceleration at the time t; the parameter β is determined according to the accuracy and stability requirements of the integration.
It is worth mentioning that for equation (9), if the stiffness matrix term is not considered, the solution of the statics problem can be achieved.
In one embodiment, the step of calculating the displacement field of the material according to the strain-rotation and decomposition theorem and the gridless method further comprises:
according to the Newmark method, the relationship among the displacement increment, the speed and the acceleration at different moments is obtained:
Figure BDA0002849582520000083
wherein,t+Δta is the acceleration at time t + Deltat,tA is the acceleration at time t,t+ΔtV is the speed at time t + Deltat,tV is the speed at time t,tAnd delta u is the displacement increment at the moment t, delta t is the time increment, and the parameters beta and gamma are determined according to the requirements of the integral precision and stability.
In one embodiment, the step of calculating the displacement field of the material according to the strain-rotation and decomposition theorem and the gridless method comprises:
obtaining an initial equivalent stiffness matrix, an initial mass matrix, an initial equivalent force vector and an initial displacement, and calculating an initial acceleration according to the following formula:
A＝M-1(Fe-Keu0) (13)
wherein A is initial acceleration, KeIs an initial equivalent stiffness matrix, M is an initial mass matrix, FeIs an initial equivalent force vector, u0Is the initial displacement.
Since the least squares approximation does not have the property of the kronecker function, the application of the intrinsic boundary condition can adopt a penalty function method or a lagrange multiplier method, and if the penalty function method is adopted, the formula (6) needs to be modified appropriately.
In one embodiment, the step of calculating the displacement field of the material according to the strain-rotation and decomposition theorem and the gridless method further comprises:
obtaining a correction increment variational equation according to the strain-rotation and decomposition theorem:
Figure BDA0002849582520000091
wherein the prime symbol "-" represents the physical quantity under the initial tow at time t in the tow coordinate system,
Figure BDA0002849582520000092
is a stress tensor,
Figure BDA0002849582520000093
Is an elastic matrix; Δ t is the time increment;
Figure BDA0002849582520000094
is the strain rate;
Figure BDA0002849582520000095
a rate of rotation;t+ΔtWineandt+ΔtWextrespectively is inertia force power and external force virtual power; delta is a variation sign; Ω is the area of the material;
Figure BDA00028495825200000914
is a penalty factor; s is a boundary; vi、VjIs the speed;
Figure BDA0002849582520000096
is the velocity determined by the boundary conditions.
Step S102: calculating a stress-strain field according to the calculation result of the displacement field;
in one embodiment, the step of calculating a stress-strain field from the calculation of the displacement field includes:
and acquiring a stress tensor according to the strain tensor by adopting the following formula:
Figure BDA0002849582520000097
wherein,
Figure BDA0002849582520000098
is the stress tensor;
Figure BDA0002849582520000099
is an elastic matrix;
Figure BDA00028495825200000910
is the strain tensor;
wherein the strain tensor
Figure BDA00028495825200000911
Calculated according to the formula (2).
Step S103: and calculating the crack propagation direction according to the calculation result of the stress-strain field and by combining a crack propagation criterion.
The prediction method provided by the application can adapt to different crack propagation criteria, and the maximum shear stress criteria are taken as an example for illustration. To obtain a smoother crack propagation path, the crack propagation direction is calculated using the weighted stresses of all integration points near the crack tip. As shown in fig. 3, the integration point corresponding to the black solid node is the integration point considered in the crack propagation direction calculation.
In one embodiment, the step of calculating the crack propagation direction according to the calculation result of the stress-strain field and the crack propagation criterion includes:
and obtaining the main stress and the main direction thereof by adopting the following formula according to the weighted stress of the crack tip:
Figure BDA00028495825200000912
wherein σ1＞σ2，
Figure BDA00028495825200000913
Wherein σtipWeighted stress for crack tip; sigma1、σ2Is principal stress, n1，n2Is the main direction.
In one embodiment, the weighted stress of the crack tip is calculated according to the following formula:
Figure BDA0002849582520000101
wherein,
Figure BDA0002849582520000102
wherein σISelecting the stress of a node near the tip of the crack; n isqThe number of integration points near the crack tip; w is aIIs the weight of the integral point; dmaxMaximum value of the distance between the crack tip and the selected node in the vicinity thereof, dIFor the corresponding integration point xiA distance from the tip of the crack, and dIAnd dmaxSatisfy the relationship 0 < dI≤2dmax。
Fig. 4 is a schematic flowchart of a prediction method according to a second embodiment of the present application. As shown in fig. 4, the prediction method provided by the present application includes the following steps:
step S201: performing node dispersion on the problem domain of the material;
step S202: creating a background integral grid and integral points according to the node discrete domain;
step S203: performing visual processing and display on the crack model;
step S204: calculating a displacement field;
step S205: calculating a stress-strain field by using the calculation result of the displacement field;
step S206: judging whether the crack length reaches the maximum value;
if the crack length reaches the maximum value, the process proceeds to step S207: stopping the numerical calculation process; if the crack length does not reach the maximum value, the step S208 is executed; optionally, whether the crack length reaches the maximum value is judged by using a crack propagation criterion, and if the crack length reaches the maximum value, the crack is judgedtip≤σ0The crack length reaches a maximum value, if σtip＞σ0If so, the crack length does not reach the maximum value; wherein σtipThe weighted stress of the crack tip is calculated according to the formula (17); sigma0The fracture critical shear stress of the material is measured by experiments;
step S208: judging whether the position of the crack tip reaches the boundary of the problem domain;
if the crack tip position reaches the problem domain boundary, the process proceeds to step S207: if the crack tip position does not reach the problem domain boundary, the step S209 is carried out; optionally, whether the crack tip position reaches the problem domain boundary is judged by using a coordinate system, and if the coordinate of the crack tip position coincides with the coordinate of any point on the problem domain boundary, the crack tip position reaches the problem domain boundary.
Step S209: calculating the crack propagation direction by using a crack propagation criterion based on the calculation result of the stress-strain field;
step S210: the position of the new crack tip is determined and the process returns to step S202.
With the loop from step S202 to step S210, the displacement field in step S204 is also calculated with loop iteration, which is specifically described with reference to the third embodiment; step S205 and step S209 can refer to the first embodiment, and are not described herein again.
Fig. 5 is a schematic specific flowchart of a prediction method according to a third embodiment of the present application. As shown in fig. 5, the prediction method provided by the present application includes the following steps:
step S301: acquiring an initial characteristic matrix, initial displacement and initial speed, and obtaining initial acceleration according to the initial characteristic matrix and the initial displacement;
wherein the initial characteristic matrix comprises an initial equivalent stiffness matrix KeInitial quality matrix M, initial equivalent force vector Fe(ii) a In particular, according to an initial equivalent stiffness matrix KeInitial quality matrix M, initial equivalent force vector FeInitial displacement u0Calculating the initial acceleration A by adopting a formula (13) in the first embodiment; the initial velocity belongs to the initial condition, and normally the structure is stationary at the initial moment, so the initial velocity V is typically zero.
Step S302: selecting a time step length and calculating an integral constant;
specifically, a time step delta t is selected, a parameter beta is combined, and an integral constant A is calculated1、A2、A3；
Wherein,
Figure BDA0002849582520000111
the parameter β is determined according to the accuracy and stability requirements of the integration.
Step S303: acquiring a stiffness characteristic matrix at the time t, and acquiring an equivalent stiffness matrix at the time t + delta t according to the stiffness characteristic matrix and an integral constant at the time t; acquiring an equivalent force characteristic matrix at the time t, the speed at the time t and the acceleration at the time t, and acquiring an equivalent force vector at the time t + delta t according to the equivalent force characteristic matrix at the time t, the speed at the time t, the acceleration at the time t and an integral constant;
specifically, according to the stiffness characteristic matrix and the integral constant at the time t, the formula (10) in the first embodiment is adopted to obtain an equivalent stiffness matrix at the time t + Δ t; and obtaining an equivalent force vector at the t + delta t moment by adopting the formula (11) in the first embodiment according to the equivalent force characteristic matrix at the t moment, the speed at the t moment, the acceleration at the t moment and the integral constant. Wherein the stiffness characteristic matrix at the time t comprises a linear stiffness matrix at the time ttKLTime t non-linear stiffness matrixtKNPenalty function stiffness matrix KαAnd the quality matrix at time ttM; the isodynamic characteristic matrix at the time t comprises an external force vector at the time t + delta tt+ΔtR, t stress integration vectortF. Penalty function force vector FαAnd the quality matrix at time ttM。
Step S304: obtaining a displacement increment at the t moment according to the equivalent stiffness matrix at the t + delta t moment and the equivalent force vector at the t + delta t moment; acquiring the displacement at the time t, and acquiring the displacement at the time t + delta t according to the displacement at the time t and the displacement increment at the time t; obtaining the speed at the t + delta t moment and the acceleration at the t + delta t moment according to the displacement increment at the t moment, the speed at the t moment and the acceleration at the t moment;
specifically, according to the equivalent stiffness matrix and the equivalent force vector at the time t + Δ t, the displacement increment at the time t is obtained by adopting the formula (9) in the first embodimenttΔ u; displacement and displacement increment according to t timetΔ u, resulting in time t + Δ tWherein the displacement at time t + Δ tt+ΔtDisplacement at time u-ttDisplacement increment at time u + ttΔ u; increment of displacement according to t timetΔ u, speedtV and accelerationtA, adopting the formula (12) in the first embodiment, the acceleration at the time of t + delta t is obtainedt+ΔtA. Speed of rotationt+ΔtV。
Step S305: and (4) a time step loop, wherein t is t + delta t, and the step S303 is executed in a returning mode.
Specifically, the displacement at time t + Δ t is updatedt+Δtu, accelerationt+ΔtA. Speed of rotationt+ΔtThe stiffness characteristic matrix and the equivalent effect characteristic matrix at the time V, t + delta t are obtained by executing the step S303 and the step S304 to obtain the node displacement at the time t +2 delta tt+2Δtu, accelerationt+2ΔtA. Speed of rotationt+2ΔtAnd V, repeating the steps to realize the cyclic iterative computation of the displacement field.
The prediction method provided by the application is based on a gridless method, combines the strain-rotation and decomposition theorem to form a new prediction method of the crack propagation direction, effectively improves the accuracy of numerical analysis of the nonlinear deformation fracture problem, and provides a new idea and means for the numerical analysis of the fracture problem.
Fig. 6 is a schematic structural diagram of a prediction apparatus according to a fourth embodiment of the present application. As shown in fig. 6, the prediction apparatus of this embodiment includes: a processor 110, a memory 111 and a computer program 112 stored in said memory 111 and executable on said processor 110. The processor 110, when executing the computer program 112, implements the steps in the various prediction method embodiments described above, such as the steps S101 to S103 shown in fig. 1.
The prediction device may include, but is not limited to, a processor 110, a memory 111. Those skilled in the art will appreciate that fig. 6 is only an example of a terminal and is not intended to be limiting and may include more or fewer components than those shown, or some components may be combined, or different components, e.g., the terminal may also include input-output devices, network access devices, buses, etc.
The Processor 110 may be a Central Processing Unit (CPU), other general purpose Processor, a Digital Signal Processor (DSP), an Application Specific Integrated Circuit (ASIC), an off-the-shelf Programmable Gate Array (FPGA) or other Programmable logic device, discrete Gate or transistor logic device, discrete hardware component, etc. A general purpose processor may be a microprocessor or the processor may be any conventional processor or the like.
The storage 111 may be an internal storage unit of the terminal, such as a hard disk or a memory of the terminal. The memory 111 may also be an external storage device of the terminal, such as a plug-in hard disk, a Smart Media Card (SMC), a Secure Digital (SD) Card, a Flash memory Card (Flash Card), and the like provided on the terminal. Further, the memory 111 may also include both an internal storage unit and an external storage device of the terminal. The memory 111 is used for storing the computer program and other programs and data required by the terminal. The memory 111 may also be used to temporarily store data that has been output or is to be output.
The present application also provides a computer-readable storage medium having stored thereon a computer program which, when executed by a processor, implements the steps of the prediction method as described above.
It should be understood that, the sequence numbers of the steps in the foregoing embodiments do not imply an execution sequence, and the execution sequence of each process should be determined by its function and inherent logic, and should not constitute any limitation to the implementation process of the embodiments of the present application.
The above embodiments are merely illustrative of the principles and utilities of the present application and are not intended to limit the application. Any person skilled in the art can modify or change the above-described embodiments without departing from the spirit and scope of the present application. Accordingly, it is intended that all equivalent modifications or changes which can be made by those skilled in the art without departing from the spirit and technical concepts disclosed in the present application shall be covered by the claims of the present application.

Claims (14)

1. A method for predicting a crack propagation direction, comprising:
calculating the displacement field of the material according to the strain-rotation and decomposition theorem and a gridless method;
calculating a stress-strain field according to the calculation result of the displacement field;
and calculating the crack propagation direction according to the calculation result of the stress-strain field and by combining a crack propagation criterion.
2. The prediction method of claim 1, wherein the strain-rotation and decomposition theorem comprises:
F＝S+R
Figure FDA0002849582510000011
Figure FDA0002849582510000012
wherein,
Figure FDA0002849582510000013
Figure FDA0002849582510000014
wherein F is a reversible linear differential transformation function; s, S,
Figure FDA0002849582510000015
Is the strain tensor; r, R,
Figure FDA0002849582510000016
Is the tensor of rotation; u. ofi|jIs a displacement component uiTo coordinate xjA derivative of (a);
Figure FDA0002849582510000017
is ui|jThe transposed matrix of (2);
Figure FDA0002849582510000018
is a kronecker function;
Figure FDA0002849582510000019
is a unit vector of the rotating shaft; theta is the average angle of full rotation.
3. The prediction method of claim 2, wherein the step of calculating the displacement field of the material according to the strain-rotation and decomposition theorem and a gridless method comprises:
obtaining an incremental variational equation according to the strain-rotation and decomposition theorem:
Figure FDA00028495825100000110
wherein the prime symbol "-" represents the physical quantity under the initial tow at time t in the tow coordinate system,
Figure FDA00028495825100000111
is a stress tensor,
Figure FDA00028495825100000112
Is an elastic matrix; Δ t is the time increment;
Figure FDA00028495825100000113
is the strain rate;
Figure FDA00028495825100000114
is the rate of rotation;t+Δ tWineandt+ΔtWextrespectively is inertia force power and external force virtual power; delta is a variation sign; Ω is the area of the material;
according to the gridless method, a moving least square interpolation function is adopted to obtain the displacement increment and the speed at the time t:
Figure FDA00028495825100000115
Figure FDA0002849582510000021
wherein,tψk(x) A moving least squares function for node k at time t;tΔukandtVkrespectively the displacement increment and the speed of the node k at the time t;tΔ u is the displacement increment at the time t,tAnd V is the speed at the time t.
4. The prediction method according to claim 3, wherein the step of calculating the displacement field of the material according to the strain-rotation and decomposition theorem and a gridless method comprises:
substituting the displacement increment of the time t and the speed of the time t into the increment variational equation to obtain a system discrete equation:
t+ΔtKe tΔu＝t+ΔtFe
wherein,
Figure FDA0002849582510000022
Figure FDA0002849582510000023
wherein,t+ΔtKe、t+ΔtFeat times t + Δ t, respectivelyAn effective stiffness matrix and an equivalent force vector;tKL、tKN、tM、tf is a linear stiffness matrix, a nonlinear stiffness matrix, a mass matrix and a stress integral vector at the moment t respectively; kαFor the penalty function stiffness matrix, FαIs a penalty function force vector;t+Δtr is an external force vector at the moment of t + delta t;ta is the acceleration at the time t; the parameter β is determined according to the accuracy and stability requirements of the integration.
5. The prediction method of claim 1, wherein the step of calculating the displacement field of the material according to strain-rotation and decomposition theorem and a gridless method further comprises:
according to the Newmark method, the relationship among the displacement increment, the speed and the acceleration at different moments is obtained:
Figure FDA0002849582510000024
wherein,t+Δta is the acceleration at time t + Deltat,tA is the acceleration at time t,t+ΔtV is the speed at time t + Deltat,tV is the speed at time t,tAnd delta u is the displacement increment at the moment t, delta t is the time increment, and the parameters beta and gamma are determined according to the requirements of the integral precision and stability.
6. The prediction method according to any one of claims 4 to 5, wherein the step of calculating the displacement field of the material according to the strain-rotation and decomposition theorem and a gridless method comprises:
obtaining an initial equivalent stiffness matrix, an initial mass matrix, an initial equivalent force vector and an initial displacement, and calculating an initial acceleration according to the following formula:
A＝M-1(Fe-Keu0)
wherein A is the initial acceleration, KeIs the initial equivalent stiffness matrix, M is the initial mass matrix, FeIs the initial equivalent force vector, u0Is the initial displacement.
7. The prediction method of claim 1, wherein the step of calculating the displacement field of the material according to strain-rotation and decomposition theorem and a gridless method further comprises:
obtaining a correction increment variational equation according to the strain-rotation and decomposition theorem:
Figure FDA0002849582510000031
wherein the prime symbol "-" represents the physical quantity under the initial tow at time t in the tow coordinate system,
Figure FDA0002849582510000032
is a stress tensor,
Figure FDA0002849582510000033
Is an elastic matrix; Δ t is the time increment;
Figure FDA0002849582510000034
is the strain rate;
Figure FDA0002849582510000035
a rate of rotation;t+Δ tWineandt+ΔtWextrespectively is inertia force power and external force virtual power; delta is a variation sign; Ω is the area of the material;
Figure FDA0002849582510000036
is a penalty factor; s is a boundary; vi、VjIs the speed;
Figure FDA0002849582510000037
is the velocity determined by the boundary conditions.
8. The prediction method according to claim 2, wherein the step of calculating a stress-strain field from the calculation result of the displacement field comprises:
and acquiring a stress tensor according to the strain tensor by adopting the following formula:
Figure FDA0002849582510000038
wherein,
Figure FDA0002849582510000039
is the stress tensor;
Figure FDA00028495825100000310
is an elastic matrix;
Figure FDA00028495825100000311
is the strain tensor.
9. The prediction method according to claim 1, wherein the step of calculating the crack propagation direction based on the calculation result of the stress-strain field in combination with the crack propagation criterion comprises:
and obtaining the main stress and the main direction thereof by adopting the following formula according to the weighted stress of the crack tip:
Figure FDA00028495825100000312
wherein σ1＞σ2，
Figure FDA00028495825100000313
Wherein σtipA weighted stress of the crack tip; sigma1、σ2Is the principal stress, n1，n2Is the primary direction.
10. The prediction method according to claim 9, wherein the step of calculating the crack propagation direction based on the calculation result of the stress-strain field in combination with the crack propagation criterion comprises:
the weighted stress of the crack tip is calculated according to the following formula:
Figure FDA0002849582510000041
wherein,
Figure FDA0002849582510000042
wherein σ1Selecting the stress of a node near the tip of the crack; n isqThe number of integration points near the crack tip; w is aIIs the weight of the integral point; dmaxMaximum value of the distance between the crack tip and the selected node in the vicinity thereof, dIFor the corresponding integration point xiA distance from the tip of the crack, and dIAnd dmaxSatisfy the relationship 0 < dI≤2dmax。
11. The prediction method according to claim 1, wherein before the step of calculating a crack propagation direction in combination with a crack propagation criterion based on the calculation of the stress-strain field, comprises:
judging whether the crack length reaches the maximum value;
judging whether the position of the crack tip reaches the boundary of the problem domain of the material;
stopping the numerical calculation process if the crack length reaches the maximum value and/or the crack tip position reaches the problem domain boundary;
and if the crack length does not reach the maximum value and the crack tip position does not reach the problem domain boundary, executing the step of calculating the crack propagation direction according to the calculation result of the stress-strain field and by combining the crack propagation criterion.
12. The prediction method according to claim 1, characterized in that the prediction method comprises the steps of:
acquiring an initial characteristic matrix, initial displacement and initial speed, and obtaining initial acceleration according to the initial characteristic matrix and the initial displacement;
selecting a time step length and calculating an integral constant;
acquiring a stiffness characteristic matrix at the time t, and acquiring an equivalent stiffness matrix at the time t + delta t according to the stiffness characteristic matrix at the time t and the integral constant; acquiring an equivalent force characteristic matrix at the time t, the speed at the time t and the acceleration at the time t, and acquiring an equivalent force vector at the time t + delta t according to the equivalent force characteristic matrix at the time t, the speed at the time t, the acceleration at the time t and the integral constant;
obtaining a displacement increment at the t moment according to the equivalent stiffness matrix at the t + delta t moment and the equivalent force vector at the t + delta t moment; acquiring the displacement at the time t, and acquiring the displacement at the time t + delta t according to the displacement at the time t and the displacement increment at the time t; obtaining the speed at the t + delta t moment and the acceleration at the t + delta t moment according to the displacement increment at the t moment, the speed at the t moment and the acceleration at the t moment;
and (5) circulating the time step, wherein t is t + delta t, and returning to execute the step three.
13. A prediction apparatus of crack propagation direction comprising a memory, a processor and a computer program stored in the memory and executable on the processor, characterized in that the processor implements the steps of the prediction method according to any of claims 1 to 12 when executing the computer program.
14. A computer-readable storage medium, in which a computer program is stored which, when being executed by a processor, carries out the steps of the prediction method according to any one of claims 1 to 12.
CN202011523316.0A 2020-12-21 2020-12-21 Crack propagation direction prediction method and device and storage medium Pending CN112861305A (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011523316.0A CN112861305A (en) 2020-12-21 2020-12-21 Crack propagation direction prediction method and device and storage medium

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011523316.0A CN112861305A (en) 2020-12-21 2020-12-21 Crack propagation direction prediction method and device and storage medium

Publications (1)

Publication Number Publication Date
CN112861305A true CN112861305A (en) 2021-05-28

Family

ID=75997827

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011523316.0A Pending CN112861305A (en) 2020-12-21 2020-12-21 Crack propagation direction prediction method and device and storage medium

Country Status (1)

Country Link
CN (1) CN112861305A (en)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113984515A (en) * 2021-10-18 2022-01-28 同济大学 Method for observing initiation and expansion of hydraulic fracture by using confocal microscope
CN115906583A (en) * 2022-12-16 2023-04-04 中国人民解放军陆军工程大学 Method and system for simulating and analyzing structural integrity of explosive column based on virtual unit method

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113984515A (en) * 2021-10-18 2022-01-28 同济大学 Method for observing initiation and expansion of hydraulic fracture by using confocal microscope
CN115906583A (en) * 2022-12-16 2023-04-04 中国人民解放军陆军工程大学 Method and system for simulating and analyzing structural integrity of explosive column based on virtual unit method
CN115906583B (en) * 2022-12-16 2023-08-01 中国人民解放军陆军工程大学 Grain structure integrity simulation analysis method and system based on virtual unit method

Similar Documents

Publication Publication Date Title
Farhat et al. Dimensional reduction of nonlinear finite element dynamic models with finite rotations and energy‐based mesh sampling and weighting for computational efficiency
CN111859766B (en) Lagrange integral point finite element numerical simulation system and method of variable calculation domain
Witteveen Explicit and robust inverse distance weighting mesh deformation for CFD
Biancolini et al. Static aeroelastic analysis of an aircraft wind-tunnel model by means of modal RBF mesh updating
US20230102815A1 (en) Turbulence field update method and apparatus, and related device thereof
CN112861305A (en) Crack propagation direction prediction method and device and storage medium
Liu et al. Impact force identification via sparse regularization with generalized minimax-concave penalty
Noventa et al. A high-order discontinuous Galerkin solver for unsteady incompressible turbulent flows
Moxey et al. Optimising the performance of the spectral/hp element method with collective linear algebra operations
CN116861822B (en) Cartesian grid-based object plane boundary processing method and device
CN116702571B (en) Numerical simulation method and device based on multiple smoothness measurement factors
Tamamidis A new upwind scheme on triangular meshes using the finite volume method
Bauer et al. Simulation of tropical-cyclone-like vortices in shallow-water ICON-hex using goal-oriented r-adaptivity
CN112733415B (en) Non-grid processing method and device for thin-wall elastomer boundary, terminal equipment and computing medium
CN114692529B (en) CFD high-dimensional response uncertainty quantification method and device, and computer equipment
CN113901552B (en) Method and device for establishing weir body model
Sonar Classical finite volume methods
JP3587827B2 (en) Airfoil performance estimation method and airfoil performance estimation program
Gan et al. Super-convergent second derivative recovery for lower-order strain gradient plasticity
CN115995278B (en) Method, device, equipment and readable storage medium for evaluating thermodynamic characteristics of material
Hashemabadi et al. Implicit second-order CUSP gridless method for unsteady moving boundary simulations
CN115995279B (en) Material mechanical property evaluation method, device, equipment and readable storage medium
CN116611369B (en) Interpolation method and device based on smoothness magnitude and candidate template point number
Van Der Helm et al. Comparison of artificial dissipation and limited flux schemes in Arbitrary Lagrangian–Eulerian finite element formulations
Leifsson et al. Fast multi-objective aerodynamic optimization using space-mapping-corrected multi-fidelity models and kriging interpolation

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