CN115083541A - Method for calculating stress intensity factor of orthotropic material based on ICVEFG method - Google Patents

Method for calculating stress intensity factor of orthotropic material based on ICVEFG method Download PDF

Info

Publication number
CN115083541A
CN115083541A CN202210528523.8A CN202210528523A CN115083541A CN 115083541 A CN115083541 A CN 115083541A CN 202210528523 A CN202210528523 A CN 202210528523A CN 115083541 A CN115083541 A CN 115083541A
Authority
CN
China
Prior art keywords
nodes
matrix
calculating
crack
icvefg
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
CN202210528523.8A
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.)
Wuhan University of Technology WUT
Original Assignee
Wuhan University of Technology WUT
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 Wuhan University of Technology WUT filed Critical Wuhan University of Technology WUT
Priority to CN202210528523.8A priority Critical patent/CN115083541A/en
Publication of CN115083541A publication Critical patent/CN115083541A/en
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16C—COMPUTATIONAL CHEMISTRY; CHEMOINFORMATICS; COMPUTATIONAL MATERIALS SCIENCE
    • G16C60/00—Computational 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
    • 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
    • G06F2113/00—Details relating to the application field
    • G06F2113/26—Composites
    • 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
    • Y—GENERAL 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
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling 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 relates to the technical field of simulation calculation in computer aided engineering, and discloses a method for calculating stress intensity factors of orthotropic materials based on an ICVEFG method, which comprises the following steps: generating a calculation model, establishing a flexibility matrix under a material coordinate system, converting the flexibility matrix into a flexibility matrix under a whole coordinate system through a conversion matrix, calculating an elastic constant matrix of the orthotropic material, establishing a crack model, processing displacement discontinuity by using a visibility criterion, obtaining a system equation and solving the system equation, outputting stress and displacement vectors of all nodes, and analyzing a stress intensity factor of the orthotropic crack-containing structure based on an ICVEFG method. The method for calculating the stress intensity factor of the orthotropic material based on the ICVEFG method solves the problems of large number of used nodes, low calculation efficiency and poor calculation accuracy in the problem of analyzing the fracture of the orthotropic material by other numerical methods.

Description

Method for calculating stress intensity factor of orthotropic material based on ICVEFG method
Technical Field
The invention relates to the technical field of simulation calculation in computer aided engineering, in particular to a method for calculating stress intensity factors of orthotropic materials based on an ICVEFG method.
Background
The composite material is formed by combining two or more materials with different physicochemical property characteristics according to a certain mode. The combined new material not only has the advantages of each component material, but also can complement each other to obtain certain characteristics which are not possessed by a single material. Therefore, the material components can be reasonably selected, and the target composite material can be designed according to the requirement. The material has the characteristics of light weight, high strength, high temperature resistance, fatigue resistance, wear resistance and the like, and plays an irreplaceable role in the fields of aerospace, ocean engineering, mechanical manufacturing, civil engineering and the like. Because the composite material is formed by combining different materials, the composite material belongs to an anisotropic material, and certain fiber reinforced composite materials even belong to an orthotropic material. However, the process for manufacturing such materials is not mature, and the component generates micro-cracks during the forming process, and the micro-cracks are continuously developed under the influence of loads and boundary conditions during the use process. In general, high specific strength materials have little plasticity and toughness and will break rapidly upon cracking. The stress intensity factor is an important parameter in linear elastic fracture mechanics. When the stress intensity factor exceeds the fracture toughness of the material itself, the crack begins to propagate. In the fields of aerospace and the like, the occurrence and the propagation of cracks directly jeopardize the service life of the component. Therefore, in order to avoid brittle fracture, it is important to study the fracture behavior of the crack-containing member and calculate the stress intensity factor of the crack-containing member.
For some simple fragmentation problems, analytical methods are the most effective methods. Irwin and Sih et al obtained stress and displacement fields near the crack tips of isotropic and anisotropic materials by a complex function method. But theoretical studies can only deal with a simple special case. Numerical methods are a better choice for solving complex fracture problems.
Common numerical methods for analyzing the fracture problem include extended finite element method (XFEM) and extended geometric analysis (XIGA). XFEM and XIGA can be seen as an improvement over the finite element method. XFEM takes advantage of the concept of unit segmentation to account for the discontinuity of the crack, making the crack independent of the mesh. In contrast to the finite element method, XFEM can characterize any form of crack and does not require re-meshing as the crack propagates. The XIGA weakens the influence of elements in the FEM on the geometry by improving the basis function of the FEM. Commonly used basis functions are non-uniform rational B-splines (NURBS) functions, T-splines, etc. At present, these two types of methods have been used to analyze Stress Intensity Factor (SIF) calculations and crack propagation simulations, etc. for isotropic and anisotropic solids. Although these methods have had great success in studying the fracture problem, these methods still have many drawbacks, either the quality of the grid has a large impact on the calculation results, the continuity between cells is low, or the calculation efficiency is too low.
The analysis of the fracture problems by the gridless method is one of the effective means for solving the problems. The problem domain is represented by nodes in a meshless method, and the discrete process is free from the constraint of units. Based on this feature, the gridless method has been widely used to analyze the fracture problem. Ventura et al proposed a new vector level set method that was used to analyze the crack propagation problem. This method is easier to simulate branch cracks and cross cracks than other methods. Nineyman et al correct the weight function in the EFG method (cell-free Galerkin method) and also make branch cracks easier to handle. Liew et al expressed cracks in unconnected cohesive segments and studied crack propagation paths in homogeneous isotropic bodies. The Yuantong and the like develop an enrichment radial base point interpolation method to research the stress field and the stress intensity factor of the crack tip. However, these studies are directed to the problem of breaking isotropic materials, and the non-grid method for breaking anisotropic materials is rarely studied. Ghorashi and Mohammadi et al correct the basis functions in EFG, apply orthotropic enrichment functions as the propagation basis functions, and analyze the stress intensity factors at the crack tips of orthotropic materials, and in their studies, too many nodes were used to discretize the problem domain. Again, this problem was studied by Fallah and Nikraftar using meshless finite volume methods. In addition, the orthotropic enrichment EFG method is also used for analyzing the fracture problem of the functional gradient material. However, these methods can be regarded as derivation methods of the EFG method, and the disadvantage that when the EFG method calculates the shape function, more nodes are needed to avoid matrix singularity is retained.
A complex variable gridless method that approximates a two-dimensional trial function using a one-dimensional basis function is an effective method for solving this problem. Compared with the traditional non-grid method, the complex variable non-grid method has the characteristics of high precision, good stability and the like due to less number of basis function terms. Currently, many Complex Variable meshless methods have been proposed, such as the modified Complex Variable Element-free Galerkin (ICVEFG), the Complex Variable quasi-convex reconstruction Kernel Particle Method (CVRKPM), and the Complex Variable meshless Boundary Element Method (CVBEM). The ICVEFG error functional has definite mathematical and physical significance, is a very effective complex variable grid-free method, and has very wide application prospect. The current ICVEFG method is mainly used to analyze the problems of thermal conduction, elastostatics and large deformations. Usually, the research objects are all isotropic materials, and the research on anisotropic materials is very rare, not to mention the method for researching the fracture problem of the anisotropic materials.
Disclosure of Invention
The invention aims to provide a method for calculating the stress intensity factor of an orthotropic material based on an ICVEFG method aiming at the defects of the technology, and solves the problems of more nodes, low calculation efficiency and poor calculation precision in the analysis of the fracture problem of the orthotropic material by other numerical methods.
In order to achieve the purpose, the method for calculating the stress intensity factor of the orthotropic material based on the ICVEFG method comprises the following steps:
A) preprocessing the model based on an ICVEFG method: inputting material parameters of an orthotropic material, generating a calculation model, arranging nodes in the calculation model, determining coordinates of boundary nodes, introducing load borne by the calculation model, generating a background integral grid in a problem domain, arranging Gaussian points in the background integral grid, and establishing coordinates, weight and Jacobian of the Gaussian points;
B) establishing a flexibility matrix under a material coordinate system, converting the flexibility matrix into a flexibility matrix under a whole coordinate system through a conversion matrix, and calculating an elastic constant matrix of the orthotropic material;
C) establishing a crack model, and processing displacement discontinuity by using a visibility criterion;
D) obtaining a system equation and solving, and outputting stress and displacement vectors of all nodes;
E) the interaction integral is calculated and the stress intensity factor is extracted.
Preferably, in the step a), the material parameters include a primary poisson's ratio, a primary elastic modulus, a secondary elastic modulus, a shear modulus, a material principal axis angle β and a size.
Preferably, the step B) includes the following steps:
B1) establishing a compliance matrix C under a material coordinate system by using the input material parameters 1 ：
Figure BDA0003645297670000041
In the formula, E 1 Principal modulus of elasticity, E 2 Is the sub-modulus of elasticity, v 12 Is a dominant poisson's ratio, G 12 Is the shear modulus;
B2) the included angle between the material coordinate system and the whole coordinate system is the material main shaft angle beta, and the flexibility matrix C under the material coordinate system is used 1 Converting the flexibility matrix C to T under the integral coordinate system T C 1 T, wherein T is a conversion matrix:
Figure BDA0003645297670000042
B3) calculating the elastic constant matrix D ═ C of orthotropic material -1 。
Preferably, the step C) includes the following steps:
C1) acquiring coordinates of starting and stopping positions of crack segments, and inputting the coordinates into a calculation model;
C2) traversing all nodes in the problem domain, judging whether the nodes are on the crack section, and if the nodes are on the crack section, deleting the nodes on the crack section;
C3) traversing all Gaussian points in the problem domain, judging whether the Gaussian points are on the crack section, and if the Gaussian points are on the crack section, deleting the Gaussian points on the crack section to obtain a crack model;
C4) if the connecting line of the node and the Gaussian point is intersected with the crack section, correcting the weight function in the ICVEFG method by adopting a visibility criterion so as to process the discontinuity of displacement;
preferably, the step D) includes the steps of:
D1) traversing all Gaussian points and nodes, calculating a corresponding shape function phi and the load of each node, and assembling a load vector and a rigidity matrix K;
D2) and (3) processing the essential boundary conditions by adopting a penalty function method: traversing all boundary nodes, calculating shape function phi of ICVEFG method of the boundary nodes, and calculating punishment rigidity matrix K of each boundary node by combining penalty factor alpha α And local penalty force vector F α And the system is superposed with a rigidity matrix K and a force vector F to form a system equation:
[K+K α ]U＝F+F α
wherein U is a node displacement parameter vector,
K＝∫ Ω B T DBdΩ
Figure BDA0003645297670000051
K α ＝α∫ Γu Φ T SΦdΓ u
Figure BDA0003645297670000052
wherein B is a strain matrix, D is an elastic constant matrix, B is a physical vector,
Figure BDA0003645297670000053
for surface force vector, Ω is problem domain, Φ is shape function, Φ (z) ═ p T (z)A- 1 (z) b (z), z ═ x + iy, (x, y) are node coordinates, p T (z) is a vector of basis functions,
Figure BDA0003645297670000054
p is a vector P of basis functions supported by nodes in the domain T (z) a matrix of (z) elements,
Figure BDA0003645297670000055
is the conjugate form of P, W (z) is a matrix of weight functions, Γ t For load boundaries, α is a penalty factor, Γ u In order to be a boundary of the displacement,
Figure BDA0003645297670000056
when a displacement boundary condition needs to be applied, s i Is 1, otherwise is 0;
D3) solving a system equation, and calculating a node displacement parameter u of each node * ；
D4) Traversing all nodes, calculating the stress and displacement of the nodes, and fitting the displacement of the nodes by the node displacement parameters of the nodes in the support domain:
Figure BDA0003645297670000057
where Re represents the real part and Im represents the imaginary part, wherein,
Figure BDA0003645297670000058
the stress sigma of any node is obtained by displacement derivation and multiplication of an elastic constant matrix D, wherein sigma (z) is DBu, and u is a node displacement vector;
D5) and outputting stress and displacement vectors of all the nodes.
Preferably, the step E) includes the following steps:
E1) inputting the size of an interaction integration region, generating a background integration grid in the interaction integration region, arranging Gaussian points in the background integration grid, and determining the coordinates, weight and Jacobian of the Gaussian points;
E2) judging whether the Gaussian point is on the crack or not, and if so, deleting the Gaussian point;
E3) traversing all Gaussian points, and calculating the shape function and the derivative of the Gaussian points;
E4) calculating the strain epsilon (z) of a Gaussian point to Bu and the stress sigma (z) to D epsilon (z), wherein the interaction integral is carried out under a crack tip coordinate system, the obtained strain and stress of the Gaussian point are values under a global coordinate system, and the stress and strain under the global coordinate system are converted into the stress and strain under a crack tip local coordinate system;
E5) calculate the value of the interaction integral:
Figure BDA0003645297670000061
a denotes the area of the M integration region, q denotes a weight function, q is 0 at the integration boundary, q is 1 at the center point, aux denotes the state of the auxiliary field, the analytical solution of stress and displacement at the tip of the auxiliary field crack, W (1,2) In order to interact with the strain energy,
Figure BDA0003645297670000062
E6) solving a characteristic equation: the characteristic equation form of the fourth-order governing partial differential equation derived from the compatibility equation and the balance equation is as follows:
Figure BDA0003645297670000063
wherein, c ij For elements in the compliance matrix C, the two imaginary parts obtained are denoted as r 1 And r 2 ；
E7) The value of the stress intensity factor is calculated from the interaction integrals,
Figure BDA0003645297670000064
in the formula (I), the compound is shown in the specification,
Figure BDA0003645297670000065
M (1) and M (2) Representing the interaction integral values obtained from different auxiliary fields, and separating the I-type stress intensity factor K from the interaction integral by using the relation between the mixed mode stress intensity factor and the interaction integral I And type II stress intensity factor K II Namely, the I-type stress intensity factor K of the orthogonal anisotropic material crack-containing member based on the ICVEFG method can be output I And type II stress intensity factor K II 。
Compared with the prior art, the invention has the following advantages:
1. the problems that other methods are low in calculation efficiency, poor in calculation stability due to the fact that grid discrete problem domains are relied on and the like are solved;
2. compared with other non-grid methods, the convergence rate of related parameters of the non-grid ICVEFG method is higher;
3. compared with other numerical methods, the grid-free ICVEFG method has the advantages that the calculated result is better matched with the reference result, the calculation precision is higher, the number of used nodes is less, the calculation can be effectively completed at high speed, and the method has good application prospect;
4. by modifying the material parameters and the material direction main axis of the orthotropic material, the fracture problems of different orthotropic materials can be flexibly and simply analyzed, so that the fracture problems of the orthotropic material in the actual engineering can be well analyzed, and the method has good practical engineering significance.
Drawings
FIG. 1 is a schematic flow chart of a method for calculating stress intensity factors of orthotropic materials based on an ICVEFG method according to the present invention;
FIG. 2 is a schematic diagram of a model used in the present embodiment;
FIG. 3 is a node layout diagram of the discrete problem domain in the present embodiment;
FIG. 4 is a schematic diagram of a material coordinate system and a global coordinate system in the present embodiment;
FIG. 5 is a graph illustrating the I-shaped regular stress intensity factors calculated in this embodiment when the elastic principal axes of the material are different;
fig. 6 is a type II canonical stress intensity factor of the example calculated in this embodiment when the material elastic principal axis direction is different.
Detailed Description
The invention is described in further detail below with reference to the figures and the specific embodiments.
As shown in fig. 1, a method for calculating the stress intensity factor of orthotropic material based on ICVEFG method includes the following steps:
A) pretreating the model based on an ICVEFG method: inputting material parameters of the orthotropic material, generating a calculation model, wherein the material parameters comprise a main Poisson ratio, a main elastic modulus, a secondary elastic modulus, a shear modulus, a material main shaft angle beta and a size, the object in the example is an infinite large unidirectional tension flat plate containing a crack, the structural geometric model and boundary conditions are shown in figure 2, the length of the model is l, the width of the model is w, the length of the crack is a, an interaction integration area is a square area with the side length of b, the material main shaft direction angle beta is 0 degree, 30 degrees, 45 degrees, 60 degrees and 90 degrees, as shown in figure 3, in the calculation model, nodes of a discrete problem domain are arranged, coordinates of boundary nodes are determined, load borne by the calculation model is led in, a background integral grid in the problem domain is generated, Gaussian points in the background integral grid are arranged, and coordinates, weight and Jacobian of the Gaussian points are established;
B) establishing a compliance matrix under a material coordinate system, converting the compliance matrix into a compliance matrix under an integral coordinate system through a conversion matrix, and calculating an elastic constant matrix of the orthotropic material:
B1) establishing a compliance matrix C under a material coordinate system by using the input material parameters 1 ：
Figure BDA0003645297670000081
In the formula, E 1 Principal modulus of elasticity, E 2 Is the sub-modulus of elasticity, v 12 Is a dominant poisson's ratio, G 12 Is the shear modulus;
B2) as shown in fig. 4, the included angle between the material coordinate system and the global coordinate system is the material principal axis angle β, and the compliance matrix C under the material coordinate system is represented by 1 Converting the flexibility matrix C to T under the global coordinate system T C 1 T, wherein T is a conversion matrix:
Figure BDA0003645297670000082
B3) calculating the elastic constant matrix D ═ C of orthotropic material -1 ；
C) A crack model is established and displacement discontinuities are treated using visibility criteria:
C1) acquiring coordinates of starting and stopping positions of crack segments, and inputting the coordinates into a calculation model;
C2) traversing all nodes in the problem domain, judging whether the nodes are on the crack section, and if the nodes are on the crack section, deleting the nodes on the crack section;
C3) traversing all Gaussian points in the problem domain, judging whether the Gaussian points are on the crack section, and if the Gaussian points are on the crack section, deleting the Gaussian points on the crack section to obtain a crack model;
C4) if the connecting line of the node and the Gaussian point is intersected with the crack section, correcting the weight function in the ICVEFG method by adopting a visibility criterion so as to process the discontinuity of displacement;
D) obtaining a system equation and solving, and outputting stress and displacement vectors of all nodes;
D1) traversing all Gaussian points and nodes, calculating a corresponding shape function phi and the load of each node, and assembling a load vector and a rigidity matrix K;
D2) and (3) processing the essential boundary conditions by adopting a penalty function method: traversing all boundary nodes, calculating shape function phi of ICVEFG method of the boundary nodes, and calculating punishment rigidity matrix K of each boundary node by combining penalty factor alpha α And local penalty force vector F α And the system is superposed with a rigidity matrix K and a force vector F to form a system equation:
[K+K α ]U＝F+F α
wherein U is a node displacement parameter vector,
K＝∫ Ω B T DBdΩ
Figure BDA0003645297670000091
K α ＝α∫ Γu Φ T SΦdΓ u
Figure BDA0003645297670000092
wherein B is a strain matrix, D is an elastic constant matrix, B is a physical vector,
Figure BDA0003645297670000093
for surface force vector, Ω is problem domain, Φ is shape function, Φ (z) ═ p T (z)A- 1 (z) b (z), z ═ x + iy, (x, y) are node coordinates, p T (z) is a vector of basis functions,
Figure BDA0003645297670000094
p is a vector P of basis functions supported by nodes in the domain T (z) a matrix of a plurality of,
Figure BDA0003645297670000096
is PW (z) is a weight function matrix, Γ t For load boundaries, α is a penalty factor, Γ u In order to be a boundary of the displacement,
Figure BDA0003645297670000095
when a displacement boundary condition needs to be applied, s i Is 1, otherwise is 0;
D3) solving a system equation, and calculating a node displacement parameter u of each node * ；
D4) Traversing all nodes, calculating the stress and displacement of the nodes, and fitting the displacement of the nodes by the node displacement parameters of the nodes in the support domain:
Figure BDA0003645297670000101
where Re represents the real part and Im represents the imaginary part, wherein,
Figure BDA0003645297670000102
the stress sigma of any node is obtained by displacement derivation and multiplication of an elastic constant matrix D, wherein sigma (z) is DBu, and u is a node displacement vector;
D5) outputting stress and displacement vectors of all nodes;
E) calculating an interaction integral and extracting a stress intensity factor:
E1) inputting the size of an interaction integration region, generating a background integration grid in the interaction integration region, arranging Gaussian points in the background integration grid, and determining the coordinates, weight and Jacobian of the Gaussian points;
E2) judging whether the Gaussian point is on the crack or not, and if so, deleting the Gaussian point;
E3) traversing all Gaussian points, and calculating the shape function and the derivative of the Gaussian points;
E4) calculating the strain epsilon (z) of a Gaussian point to Bu and the stress sigma (z) to D epsilon (z), wherein the interaction integral is carried out under a crack tip coordinate system, the obtained strain and stress of the Gaussian point are values under a global coordinate system, and the stress and strain under the global coordinate system are converted into the stress and strain under a crack tip local coordinate system;
E5) calculate the value of the interaction integral:
Figure BDA0003645297670000103
a denotes the area of the M integration region, q denotes a weight function, q is 0 at the integration boundary, q is 1 at the center point, aux denotes the state of the auxiliary field, the analytical solution of stress and displacement at the tip of the auxiliary field crack, W (1,2) In order to interact with the strain energy,
Figure BDA0003645297670000104
E6) solving a characteristic equation: the characteristic equation form of the fourth-order governing partial differential equation derived from the compatibility equation and the balance equation is as follows:
c 11 r 4 -2c 13 r 3 +(2c 12 +c 33 )r 2 -2c 23 r+c 22 ＝0，
wherein, c ij For elements in the compliance matrix C, the two imaginary parts obtained are denoted as r 1 And r 2 ；
E7) The value of the stress intensity factor is calculated from the interaction integral,
Figure BDA0003645297670000111
in the formula (I), the compound is shown in the specification,
Figure BDA0003645297670000112
M (1) and M (2) Representing the interaction integral values obtained from different auxiliary fields, and separating the I-type stress intensity factor K from the interaction integral by using the relation between the mixed mode stress intensity factor and the interaction integral I And type II stress intensity factor K II That is to sayOutputting I-type stress intensity factor K of orthogonal anisotropic material crack-containing member based on ICVEFG method I And type II stress intensity factor K II 。
Fig. 3 shows a layout of nodes for analyzing the fracture problem of the orthotropic material based on the ICVEFG method and the EFG method, and table 1 shows the numbers of nodes and integration grids used for analyzing the fracture problem of the orthotropic material based on the ICVEFG method and the EFG method. As can be seen from fig. 3 and table 1, the number of nodes and the number of integration grids used in the method are smaller than those used in other non-grid methods, so that the method can have higher efficiency in analyzing the fracture problem of the orthotropic material.
TABLE 1 non-grid method for calculating parameter values
Figure BDA0003645297670000113
FIG. 5 and FIG. 6 show the type I stress intensity factor K in different material principal axis directions calculated based on ICVEFG method I And type II stress intensity factor K II And comparison with other numerical methods. It can be seen that the values of the stress intensity factor calculated based on the ICVEFG method are well within the range of the results calculated by other numerical methods, and thus, the ICVEFG method can obtain higher accuracy in analyzing the orthotropic fracture problem.
The method for calculating the stress intensity factor of the orthotropic material based on the ICVEFG method solves the problems that other methods are low in calculation efficiency, poor in calculation stability due to the fact that grid discrete problem domains are relied on, and the like; compared with other non-grid methods, the convergence rate of related parameters of the non-grid ICVEFG method is higher; compared with other numerical methods, the grid-free ICVEFG method has the advantages that the calculated result is better matched with the reference result, the calculation precision is higher, the number of used nodes is less, the calculation can be effectively completed at high speed, and the method has good application prospect; in addition, the fracture problem of different orthotropic materials can be flexibly and simply analyzed by modifying the material parameters and the material direction main axis of the orthotropic material, so that the fracture problem of the orthotropic material in the actual engineering can be well analyzed, and the method has good practical engineering significance.

Claims (6)

1. A method for calculating stress intensity factors of orthotropic materials based on an ICVEFG method is characterized in that: the method comprises the following steps:
A) pretreating the model based on an ICVEFG method: inputting material parameters of an orthotropic material, generating a calculation model, arranging nodes in the calculation model, determining coordinates of boundary nodes, importing load borne by the calculation model, generating a background integral grid in a problem domain, arranging Gaussian points in the background integral grid, and establishing coordinates, weight and Jacobian of the Gaussian points;
B) establishing a flexibility matrix under a material coordinate system, converting the flexibility matrix into a flexibility matrix under an integral coordinate system through a conversion matrix, and calculating an elastic constant matrix of the orthotropic material;
C) establishing a crack model, and processing displacement discontinuity by using a visibility criterion;
D) obtaining a system equation and solving, and outputting stress and displacement vectors of all nodes;
E) the interaction integral is calculated and the stress intensity factor is extracted.
2. The method of claim 1 for calculating the stress intensity factor of orthotropic materials based on the ICVEFG method, wherein: in the step A), the material parameters comprise a main Poisson's ratio, a main elastic modulus, a secondary elastic modulus, a shear modulus, a material main axis angle beta and a size.
3. The method of claim 2 for calculating the stress intensity factor of orthotropic material based on ICVEFG method, wherein: the step B) comprises the following steps:
B1) establishing a compliance matrix C under a material coordinate system by using the input material parameters 1 ：
Figure FDA0003645297660000011
In the formula, E 1 Principal modulus of elasticity, E 2 Is the sub-modulus of elasticity, v 12 Is a dominant poisson's ratio, G 12 Is the shear modulus;
B2) the included angle between the material coordinate system and the whole coordinate system is the material main shaft angle beta, and the flexibility matrix C under the material coordinate system is used 1 Converting the flexibility matrix C to T under the integral coordinate system T C 1 T, wherein T is a conversion matrix:
Figure FDA0003645297660000021
B3) calculating the elastic constant matrix D ═ C of orthotropic material -1 。
4. The method of claim 3 for calculating the stress intensity factor of orthotropic material based on ICVEFG method, wherein: the step C) comprises the following steps:
C1) acquiring coordinates of starting and stopping positions of crack segments, and inputting the coordinates into a calculation model;
C2) traversing all nodes in the problem domain, judging whether the nodes are on the crack section, and if the nodes are on the crack section, deleting the nodes on the crack section;
C3) traversing all Gaussian points in the problem domain, judging whether the Gaussian points are on the crack section, and if the Gaussian points are on the crack section, deleting the Gaussian points on the crack section to obtain a crack model;
C4) and if the connecting line of the node and the Gaussian point is intersected with the crack section, correcting the weight function in the ICVEFG method by adopting a visibility criterion so as to process the displacement discontinuity.
5. The method of claim 4 for calculating the stress intensity factor of orthotropic materials based on the ICVEFG method, wherein: the step D) comprises the following steps:
D1) traversing all Gaussian points and nodes, calculating a corresponding shape function phi and the load of each node, and assembling a load vector and a rigidity matrix K;
D2) and (3) processing the essential boundary conditions by adopting a penalty function method: traversing all boundary nodes, calculating the shape function phi of the ICVEFG method of the boundary nodes, and calculating the punishment rigidity matrix K of each boundary node by combining the penalty factor alpha α And local penalty force vector F α And the system is superposed with a rigidity matrix K and a force vector F to form a system equation:
[K+K α ]U＝F+F α
wherein U is a node displacement parameter vector,
K＝∫ Ω B T DBdΩ
Figure FDA0003645297660000022
Figure FDA0003645297660000023
Figure FDA0003645297660000024
wherein B is a strain matrix, D is an elastic constant matrix, B is a physical vector,
Figure FDA0003645297660000025
for the surface force vector, Ω is the problem domain, Φ is the shape function, Φ (z) is p T (z)A -1 (z) b (z), z ═ x + iy, (x, y) are node coordinates, p T (z) is a vector of basis functions,
Figure FDA0003645297660000031
p is a vector P of basis functions supported by nodes in the domain T (z) a matrix of (z) elements,
Figure FDA0003645297660000032
is the conjugate form of P, W (z) is a matrix of weight functions, Γ t For load boundaries, α is a penalty factor, Γ u In order to be a boundary of the displacement,
Figure FDA0003645297660000033
when a displacement boundary condition needs to be applied, s i Is 1, otherwise is 0;
D3) solving a system equation, and calculating a node displacement parameter u of each node * ；
D4) Traversing all nodes, calculating the stress and displacement of the nodes, and fitting the displacement of the nodes by the node displacement parameters of the nodes in the support domain:
Figure FDA0003645297660000034
where Re represents the real part and Im represents the imaginary part, wherein,
Figure FDA0003645297660000035
the stress sigma of any node is obtained by displacement derivation and multiplication of an elastic constant matrix D, wherein sigma (z) is DBu, and u is a node displacement vector;
D5) and outputting stress and displacement vectors of all the nodes.
6. The method of claim 5 for calculating the stress intensity factor of orthotropic material based on ICVEFG method, wherein: the step E) comprises the following steps:
E1) inputting the size of an interaction integration region, generating a background integration grid in the interaction integration region, arranging Gaussian points in the background integration grid, and determining the coordinates, weight and Jacobian of the Gaussian points;
E2) judging whether the Gaussian point is on the crack or not, and if so, deleting the Gaussian point;
E3) traversing all Gaussian points, and calculating the shape function and the derivative of the Gaussian points;
E4) calculating the strain epsilon (z) of a Gaussian point to Bu and the stress sigma (z) to D epsilon (z), wherein the interaction integral is carried out under a crack tip coordinate system, the obtained strain and stress of the Gaussian point are values under a global coordinate system, and the stress and strain under the global coordinate system are converted into the stress and strain under a crack tip local coordinate system;
E5) calculate the value of the interaction integral:
Figure FDA0003645297660000041
a denotes the area of the M integration region, q denotes a weight function, q is 0 at the integration boundary, q is 1 at the center point, aux denotes the state of the auxiliary field, the analytical solution of stress and displacement at the tip of the auxiliary field crack, W (1,2) In order to interact with the strain energy of the strain,
Figure FDA0003645297660000042
E6) solving a characteristic equation: the characteristic equation form of the fourth-order governing partial differential equation derived from the compatibility equation and the balance equation is as follows:
c 11 r 4 -2c 13 r 3 +(2c 12 +c 33 )r 2 -2c 23 r+c 22 ＝0，
wherein, c ij For elements in the compliance matrix C, the two imaginary parts obtained are denoted as r 1 And r 2 ；
E7) The value of the stress intensity factor is calculated from the interaction integral,
Figure FDA0003645297660000043
in the formula (I), the compound is shown in the specification,
Figure FDA0003645297660000044
M (1) and M (2) Representing the interaction integral values obtained from different auxiliary fields, and separating the I-type stress intensity factor K from the interaction integral by using the relation between the mixed mode stress intensity factor and the interaction integral I And type II stress intensity factor K II Namely, the I-type stress intensity factor K of the orthogonal anisotropic material crack-containing member based on the ICVEFG method can be output I And type II stress intensity factor K II 。
CN202210528523.8A 2022-05-16 2022-05-16 Method for calculating stress intensity factor of orthotropic material based on ICVEFG method Pending CN115083541A (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210528523.8A CN115083541A (en) 2022-05-16 2022-05-16 Method for calculating stress intensity factor of orthotropic material based on ICVEFG method

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210528523.8A CN115083541A (en) 2022-05-16 2022-05-16 Method for calculating stress intensity factor of orthotropic material based on ICVEFG method

Publications (1)

Publication Number Publication Date
CN115083541A true CN115083541A (en) 2022-09-20

Family

ID=83246671

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210528523.8A Pending CN115083541A (en) 2022-05-16 2022-05-16 Method for calculating stress intensity factor of orthotropic material based on ICVEFG method

Country Status (1)

Country Link
CN (1) CN115083541A (en)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116090113A (en) * 2022-10-08 2023-05-09 盐城工业职业技术学院 Variable node suspension unit grid finite element model of block structure
CN117371272A (en) * 2023-09-22 2024-01-09 天津大学 Method for calculating crack length and fracture performance of clamping type unilateral notch tensile test sample applicable to different anisotropic materials and sizes
CN117371271A (en) * 2023-09-22 2024-01-09 天津大学 Anisotropic material crack length and fracture performance testing method based on pin shaft type unilateral notch tensile test sample
CN117494486A (en) * 2024-01-03 2024-02-02 南通泰胜蓝岛海洋工程有限公司 Combined beam refinement stress displacement analysis method under local concentrated load effect

Cited By (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116090113A (en) * 2022-10-08 2023-05-09 盐城工业职业技术学院 Variable node suspension unit grid finite element model of block structure
CN116090113B (en) * 2022-10-08 2023-09-08 盐城工业职业技术学院 Variable node suspension unit grid finite element model of block structure
CN117371272A (en) * 2023-09-22 2024-01-09 天津大学 Method for calculating crack length and fracture performance of clamping type unilateral notch tensile test sample applicable to different anisotropic materials and sizes
CN117371271A (en) * 2023-09-22 2024-01-09 天津大学 Anisotropic material crack length and fracture performance testing method based on pin shaft type unilateral notch tensile test sample
CN117371272B (en) * 2023-09-22 2024-04-19 天津大学 Method for calculating crack length and fracture performance of clamping type unilateral notch tensile test sample applicable to different anisotropic materials and sizes
CN117494486A (en) * 2024-01-03 2024-02-02 南通泰胜蓝岛海洋工程有限公司 Combined beam refinement stress displacement analysis method under local concentrated load effect
CN117494486B (en) * 2024-01-03 2024-04-02 南通泰胜蓝岛海洋工程有限公司 Combined beam refinement stress displacement analysis method under local concentrated load effect

Similar Documents

Publication Publication Date Title
CN115083541A (en) Method for calculating stress intensity factor of orthotropic material based on ICVEFG method
Moorthy et al. A Voronoi cell finite element model for particle cracking in elastic-plastic composite materials
Drugan et al. Asymptotic analysis of growing plane strain tensile cracks in elastic-ideally plastic solids
Chen et al. The enhanced extended finite element method for the propagation of complex branched cracks
CN107169191A (en) A kind of fan blade modeling method
CN112380650B (en) Design method of structural part of working device
CN104123400B (en) Global Local details finite element methods based on force method
Yu-qiu et al. Generalized conforming element for bending and buckling analysis of plates
Feng et al. A gradient weighted extended finite element method (GW-XFEM) for fracture mechanics
Gupta et al. Geometrically nonlinear bending analysis of variable stiffness composite laminated shell panels with a higher-order theory
Ma et al. A smoothed enriched meshfree Galerkin method with two-level nesting triangular sub-domains for stress intensity factors at crack tips
Wu et al. A twice-interpolation finite element method (TFEM) for crack propagation problems
Liao et al. Crack propagation modelling using the weak form quadrature element method with minimal remeshing
Xiao et al. Incremental‐secant modulus iteration scheme and stress recovery for simulating cracking process in quasi‐brittle materials using XFEM
Houari et al. Using finite element analysis to predict the damage in FGM-3D notched plate under tensile load; Different geometric concept
Bai et al. A probabilistic combined high and low cycle fatigue life prediction framework for the turbine shaft with random geometric parameters
Sun et al. Research of large scale mechanical structure crack growth method based on finite element parametric submodel
CN115114816A (en) Numerical simulation method for crack propagation of multi-interface non-uniform material under strong transient thermal load
Muravin et al. Multiple crack weight for solution of multiple interacting cracks by meshless numerical methods
CN109388833A (en) A kind of elastic element structure optimum design method based on fatigue life
CN114861494B (en) C/C composite material elastic property prediction method considering pyrolytic carbon anisotropy characteristics
CN108595724A (en) Composite material revolving meber design method
Jin et al. A numerical study on damage characteristics in composite tapered laminates under cyclic loading with different stress ratios
Tian et al. Calculation of dynamic stress intensity factors and T-stress using an improved SBFEM
Ma et al. A conforming A-FEM for modeling arbitrary crack propagation and branching in solids

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