CN108763841B - Elastic fracture simulation method based on dual boundary element and strain energy optimization analysis - Google Patents

Elastic fracture simulation method based on dual boundary element and strain energy optimization analysis Download PDF

Info

Publication number
CN108763841B
CN108763841B CN201810815571.9A CN201810815571A CN108763841B CN 108763841 B CN108763841 B CN 108763841B CN 201810815571 A CN201810815571 A CN 201810815571A CN 108763841 B CN108763841 B CN 108763841B
Authority
CN
China
Prior art keywords
boundary
fracture
crack
equation
stress
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201810815571.9A
Other languages
Chinese (zh)
Other versions
CN108763841A (en
Inventor
李帅
崔柳
杨晨
于洋
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Qingdao Research Institute Of Beihang University
Beihang University
Original Assignee
Qingdao Research Institute Of Beihang University
Beihang University
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 Qingdao Research Institute Of Beihang University, Beihang University filed Critical Qingdao Research Institute Of Beihang University
Priority to CN201810815571.9A priority Critical patent/CN108763841B/en
Publication of CN108763841A publication Critical patent/CN108763841A/en
Application granted granted Critical
Publication of CN108763841B publication Critical patent/CN108763841B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

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
    • G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

The invention provides an elastic fracture simulation method based on dual boundary elements and strain energy optimization analysis, which is used for processing two basic problems related to three-dimensional linear elastic fracture. Models containing overlapping fracture surfaces are represented by using a dual boundary element method, and a semi-analytic method is adopted to solve a super-singular boundary integral equation. The continuous properties from contact load and crack front strain energy release are represented by introducing a contact force model. The extension of the fracture end is accurately controlled by calculating a fracture factor for the extension distance of the fracture over the fracture duration.

Description

Elastic fracture simulation method based on dual boundary element and strain energy optimization analysis
Technical Field
The invention belongs to the technical field of physical simulation, and particularly relates to an elastic fracture simulation method based on dual boundary elements and strain energy optimization analysis.
Background
The boundary element method has been widely used in many fields of engineering application as a powerful numerical calculation tool. Researchers have made the Boundary element method have various applications (especially on the problems of infinite field and open surface) by substituting the basic solution for different applications into the Boundary Integral Equations (BIEs).
In terms of fracture simulation, on the basis of applying a traditional Displacement Boundary integral equation (Displacement BIEs) to one side of a fracture by a Dual Boundary Element Method (Dual-BEM), an independent stress Boundary integral equation (transformation BIEs) is applied to the other side of the fracture to solve the problem of overlapping fracture surfaces. The Dual-BEM only discretizes the surface of the model and sets different boundary integral equations for two sides of the overlapped surface at the fracture of the model to solve the problem of system degradation, which has obvious advantages in the application of linear elastic fracture simulation, but also has to face the situation of solving various singular boundary integrals.
Modeling of the model fracture process includes geometric and physical methods. The simulation of rigid body fracture can calculate how to fracture by predefined methods or dynamically calculating stress. The predefined method on the model can design a fracture pattern in advance, and the fracture pattern can be started from a certain fracture point, or the cracks are not concentrated to a certain fracture point, and the predefined fracture pattern is used for replacing the fracture pattern when the fracture occurs, so that the speed is high, but the authenticity is lacked. The method for dynamically calculating the fracture from the internal stress of the object can generate a very real fracture mode from a certain point, but needs a large amount of calculation, and is not suitable for real-time application. Brittle fracture generally refers to little elastic deformation between hard bodies during crack opening and propagation.
The main challenge of three-dimensional fracture simulation is the updating of the geometry during crack propagation. In general, the crack growth amplitude during fracture is typically two orders of magnitude smaller than the grid resolution, and moreover, given the freedom of fracture propagation, updating the geometric grid presents significant difficulties. More importantly, the complexity caused by grid updating can cause discontinuity and singularity frequently in the simulation process, which puts new requirements on the numerical solution method of fracture simulation. If the volumetric mesh subdivision method is adopted to represent the gap, the actual gap cannot be accurately met due to the limitation of the minimum granularity and the numerical continuity of the mesh, and the pathological mesh is reflected in a physical system to cause the degradation of the system. The finite element method, as a more common fracture modeling method, is very weak in supporting such problems. For a large-scale grid caused by the subdivided grid operation at the fracture, the efficiency of the finite element method is greatly reduced, and the convergence rate of the finite element method is even influenced by the grid with low quality. Furthermore, the shape function of the traditional finite element method cannot represent the singular result, and how to represent the singular stress is a big challenge faced by the finite element method.
Disclosure of Invention
In order to solve the problems in the prior art, the invention provides an elastic fracture simulation method based on dual boundary elements and strain energy optimization analysis. In the physical simulation process of the fracture, the extension of the fracture end is accurately controlled by introducing a contact force model to represent continuous attributes from contact load and crack front strain energy release, and calculating a fracture factor for the extension distance of the crack over the fracture duration.
The invention is realized by adopting the following technical scheme: an elastic fracture simulation method based on dual boundary elements and strain energy optimization analysis comprises the following steps:
step (1), carrying out dual boundary integral equation modeling on a target model and cracks thereof;
step (2), dispersing a boundary integral equation based on a segmented smooth surface;
step (3), singular integrand processing is carried out on the singular integrand in the boundary integral equation by adopting a singularity reduction technology;
step (4), solving a free term in a boundary integral equation;
step (5), in the simulation process, accurately modeling the fracture generation process by using a continuous contact model;
establishing a continuous contact time period by adopting a Hertz model contacted by an elastic ball, representing the change of stress in the contact process by defining an external force changing along with time in the contact time period, and representing the dynamic process of the external force acting on the surface of the model by placing a fracture in a continuous contact time period through a simulation time step appointed by a user and the obtained continuous contact time period;
step (6), in the simulation process, describing a stress field near the crack by using a stress intensity factor;
carrying out deformation simulation on the model under the action of continuous external force to obtain the displacement of boundary elements on the surface of the crack so as to obtain the opening displacement of the crack, establishing a local coordinate system for an extension point at the front edge of the crack, projecting the opening displacement onto a local coordinate basis vector, calculating stress intensity factors, and averaging the stress intensity factors in a plurality of cracks to obtain a final stress intensity factor when the extension point is shared by the plurality of cracks;
step (7), in the process of generating fracture, initializing fracture according to surface stress analysis, and extending on the crack based on a stress intensity factor;
aiming at the surface stress analysis of boundary elements, firstly obtaining the deformation gradient of the current element, obtaining the First Piola-Kirchhoff stress tensor of a triangular patch, and generating a new crack on the current triangular patch when the maximum stress exceeds the strength of a material; and calculating the stress intensity factor of the point on the crack, extending the crack when the effective force intensity is greater than the fracture toughness of the material, and calculating the extending speed and the extending angle in the extending process of the crack.
Further, the step (1) comprises: establishing a boundary integral equation for the region omega of the target model by adopting a Betfi-Somigliana equation; the boundary integral equation at the source point is expressed by a Kelvin basis solution formed by the distance between the source point and the field point of the domain boundary; and modeling the cracks in the domain into an open surface, modeling the three groups of surfaces into independent displacement boundary integral equations by using dual boundary elements, and deriving the source points by using a displacement boundary integral method on the overlapped surfaces to obtain a stress boundary integral equation.
Further, the step (2) comprises: performing linear dispersion on a segmented smooth surface model widely existing in the field of computer graphics; solving the boundary integral equation on the triangular patch set by adopting a normal numerical integration method; for solving the boundary integral equation of the set of the surface patches with numerical singularity, a culling domain needs to be established at a source point.
Further, the step (3) comprises: unified processing is carried out on singular integrand functions in the boundary integral equation by adopting a singularity reduction technology; taylor expansion is carried out on singular integrands under a polar coordinate system, so that singular terms of the integrands are explicitly expressed, and a normal integrand without the singular terms is obtained by subtracting the singular terms from the integrands; their analytical solutions are solved for the removed singular terms and the results are brought back to the boundary integral equation.
Further, the step (4) comprises: adopting a spherical elimination domain, and performing intersection operation on a tangent plane set where boundary elements around a source point are positioned and the spherical elimination domain to obtain a vertex set and a profile arc set which are positioned on a spherical surface; the following set of variables is defined: a tangent plane set, an outer normal set of the tangent plane set and an included angle radian set of two intersected tangent planes; and obtaining a coefficient matrix related to the set structure of the local surface where the source point is located.
Compared with the prior art, the invention has the advantages and positive effects that:
1. the method is characterized in that displacement boundary integral equations and stress boundary integral equations are respectively used for expressing two sides of overlapped fracture surfaces, and sampling points of the two overlapped surfaces are independently modeled. And for various singular integral equations generated therewith, the influence of the singular integral equations on numerical solution is eliminated by adopting a semi-analytic method and a robust integral transformation method.
2. The continuous contact force model is introduced, the process of applying the dominant load in the fracture simulation can be explicitly expressed into a continuously-changed time interval, the continuously-changed load can be obtained, the time sequence continuous control is carried out on the fracture generation process, and the fracture application occasion is expanded.
3. A method for analyzing the significance characteristics of the elongation and strain energy of the fracture front end in the fracture extension process is provided. On one hand, the calculation of parameters in the crack generation process is simplified, and the calculation efficiency is improved, and on the other hand, the influence of crack generation on the prior data can be reflected.
Drawings
FIG. 1 is a fracture simulation flow diagram;
FIG. 2 is a schematic diagram of a boundary element integral equation domain;
FIG. 3 is a schematic view of a piecewise smooth surface mesh model;
FIG. 4 is a schematic view of a local coordinate projection of a non-smooth surface;
FIG. 5 is a local coordinate system of crack propagation;
FIG. 6 is a split load mode;
FIG. 7 is a schematic diagram of global coordinate system, local coordinate system, and conformal coordinate system transformation;
FIG. 8 is a schematic of the geometry for calculating a free term consisting of five tangent planes;
fig. 9 is a graph of human femoral (thigh) fracture simulation results;
fig. 10 is a graph of the simulation result of fracture of the fifth lumbar vertebra of a human body with different fracture toughness.
Detailed Description
The invention provides an elastic fracture simulation method based on dual boundary elements and strain energy optimization analysis, and the principle of the invention is introduced below by combining specific method steps.
1) Dual boundary element method modeling of segmented smooth surface
And (1) establishing a boundary integral equation for the region omega of the target model by adopting a Betfi-Somigliana equation. The boundary integral equation at the source point is a Kelvin basis solution U formed by the distance between the source point and the field point of the domain boundaryijAnd TijThe Kelvin basis solution shows the displacement of the reaction to the site after being influenced by the material after a Unit Load (Unit Load) is applied to the source site. When the site approaches the source, the base solution UijWith Weak Singularity (Weak Singularity), such Weak singular terms can transform U by Polar Coordinate Transformation (PCT)ijInto a Integrable Integrand. T isijHas Strong Singularity (Strong singular), so Cauchy Principal Value integral (CPV) needs to be introduced for processing. When the domain Ω of the target model contains cracks, the cracks in the domain are modeled as an open surface, and the domain Ω boundary Γ becomes a triple (Γ ═ Γr+Γc++Γc-) Wherein gamma isc+,Γc-Are overlapping surfaces that are in the same position and in opposite directions, where we assume that at Γc+,Γc-No stress (Traction-Free) was applied. Using dual boundary elements(Dual-BEM) modeling the three sets of surfaces as independent Displacement Boundary Integral Equations (DBIEs), and deriving the source points by the displacement boundary integral method on the overlapped surfaces to obtain force boundary integral equations (TBIEs).
The step (2) and the boundary element method used by the invention mainly aim at the field of computer graphics, and Linear Discretization (Linear Discretization) is needed to be carried out on segmented smooth surface models widely existing in the field of computer graphics. Since a piecewise smooth surface is used as the boundary of the model, the source point s is a discontinuous (Non-smooth) point, and is shared by multiple triangular patches (see fig. 3 (a)). Since the integrand F in the boundary integral equation has an odd anisotropy with respect to the distance of the field point x (integration point) to the source point s, the piecewise continuous surface Γ is divided into two parts: singular facet set ΓsAnd a set of triangular patches ΓrAmong which a singular patch set ΓsIs composed of patches sharing a source point s (e.g. green part in fig. 3(a)), and the remaining gray part constitutes a triangular patch set Γr. The solution of the boundary integral equation on the triangular patch set adopts a normal numerical integration method, and good calculation accuracy can be obtained. For solving the boundary integral equation of the singular patch set, the method needs to establish a removed Domain (Excluded Domain) at the source point s, the region is a sphere (as shown in fig. 3(b)) with the radius of epsilon and taking the source point s as the center of the circle, and the spherical removed Domain is adopted in the method so as to utilize the removed Domain to the free term cijAnd (4) the advantages in calculation. The invention relates to a method for modeling dual boundary elements of a segmented smooth surface, which is usually used together with a spline function (NURBS), wherein the spline function has higher requirement on the representation of a model grid, and a segmented smooth grid (triangular grid) is generally used in virtual operation model grid data (or a grid model widely used in the field of graphics).
And (3) obtaining a boundary integral equation defined on the segmented smooth surface through the processing of the step (2). The equations contain various singular integrand functions. In order to be able to uniformly process these singular integrand, Singular Subtraction Techniques (SST) are used. Singular integrants in a polar coordinate system are Taylor expanded so that singular terms of the integrants can be Explicitly (Explicitly) expressed, and a Regular (Regular) integrant without the singular terms is obtained by subtracting the singular terms from the integrants. Their analytical solutions (analytical solutions) are solved for the removed singular terms and the results are brought back to the boundary integral equation. Extreme operations in the formula boundary integral equation can be avoided by using the singularity reduction method.
2) Numerical solution of boundary element equation
In the boundary element modeling process, the Somigliana identity equation generally comprises two parts, wherein one part is used for solving a kelvin basic solution with singular attributes through a boundary integral equation, and the other part is used for solving a coefficient matrix C related to the geometric structure of the local surface where the source point is located. For modeling the three-dimensional elastic fracture problem on a piecewise smooth surface to which the present invention relates, an explicit computation technique is employed that supports explicitly computing a coefficient matrix C at boundary vertices shared by any plurality of tangent planes on the piecewise smooth surface. Compared with the traditional application scene, such as CAD design, the set model in the field is basically fixed and is spliced by various basic geometric shapes, so the method for calculating the free items comprises the steps of (exhaustively) calculating in advance (using a method with large operation amount such as integration) and using the method to quote (look up table), and the method can continuously and quickly calculate any polygonal geometric structure and is well suitable for the complex geometric conditions in the fracture process.
3) Fracture simulation method based on continuous contact model
And (5) introducing a continuous contact model in order to accurately model the generation process of the fracture. A Hertz model of elastic ball contact was used to establish the continuous contact time period. The change of the stress during the Contact process is represented by defining the external force which changes along with the time, and since the Hertz Contact Model is a Non-adhesive Contact Model, the external force is regarded as a sine Function (Sinussoid Function). And (3) placing the fracture in a continuous contact time period through a simulation time step specified by a user and the obtained continuous contact time period to represent the dynamic process of the external force acting on the surface of the model.
And (6) for the crack propagation, because the Stress has singularity at the crack tip, Stress Intensity Factors (SIFs) are used for describing the Stress field near the crack. And performing deformation simulation on the model under the action of continuous external force to obtain the Displacement of the boundary element on the surface of the Crack, thereby obtaining the Crack Opening Displacement (COD). Calculating a stress intensity factor, namely establishing a local coordinate system for an extension point at the front edge of the crack, establishing mutually vertical local basis vectors for the extension point c according to the crack where the extension point c is located, the normal of a boundary element to which the crack belongs and a tangent parallel to the crack, projecting the opening displacement onto the local coordinate basis vectors, and calculating the stress intensity factor; when the extension point is shared by a plurality of cracks, the final stress intensity factor is obtained by averaging the stress intensity factors in the plurality of cracks.
Step (7), modeling of the object fracture process requires handling initial Crack Initiation and Propagation. The generation of fracture is based on the analysis of surface stress, the principal stress and the principal stress direction of all boundary elements owned by the model are calculated one by one, firstly, the Deformation Gradient (Deformation Gradient) of the current element is obtained, the First Piola-Kirchhoff stress tensor of the triangular patch is obtained, and when the maximum stress exceeds the material strength, a new fracture is inserted into the centroid of the current triangular patch by taking the direction of the maximum stress as the normal line of the fracture surface. And calculating the stress intensity factor of the point on the crack, extending the crack when the effective stress intensity is larger than the fracture toughness of the material, and calculating the extending speed and the extending angle in the extending process of the crack.
The present invention will be described in further detail with reference to the accompanying drawings and specific embodiments.
Fig. 1 shows an overall processing flow of an elastic fracture simulation process based on dual boundary elements and strain energy optimization analysis, and an embodiment of the present invention provides an elastic fracture simulation method based on dual boundary elements and strain energy optimization analysis, which includes the steps of:
for the problem of three-dimensional linear elasticity, the method adopts a Betfi-Somigliana equation to establish a boundary integral equation (see formula 1) for an area omega (shown in figure 2) of a target model.
In equation 1, the boundary integral equation at Source Point (Source Point) s can be represented by the s and field Ω boundaries
Figure BDA0001740229700000081
The distance r (r ═ s-x |) between Field points (Field points) x of (c) a Kelvin-based solution UijAnd
Figure BDA0001740229700000082
to indicate. Wherein, the Kelvin basic solution UijThe displacement of the reaction to the field point x after the influence of the material is shown after a Unit Load (Unit Load) is applied to the source point s. In the present method, when the field point x approaches the source point s (r → 0), U in equation 1ijHas Weak Singularity (Weak Singularity), Uij＝O(r-1) Such weak singular terms can be converted into Integrable Integrand (Integrable Integrand) by Polar Coordinate Transformation (PCT). Unfortunately, when r → 0, TijHaving a Strong Singularity (Strong Singularity), Tij＝O(r-2) Therefore, it is necessary to introduce Cauchy Principal Value integral (CPV, symbolized by
Figure BDA0001740229700000083
Represented) to perform the process. U and t in equation 1 represent displacement and force at the field point, respectively, and the subscript i, j represents the direction in three-dimensional space (i, j ═ 1,2, 3).
Figure BDA0001740229700000084
When the target model contains a crack in the domain Ω (as shown in fig. 2 (b)), the method models the crack as an Open Surface (Open Surface). For the case in fig. 2(b), the domain Ω boundary Γ becomes a three-part (Γ ═ Γr+Γc++Γc-) Wherein gamma isc+,Γc-Are overlapping surfaces, which are in the same position and in opposite directions, where we assume that at Γc+,Γc-No stress (Traction-Free) was applied. Dual-boundary element (Dual-BEM) is obtained by modeling three sets of surfaces as independent Displacement Boundary Integral Equations (DBIEs) and superimposing surfaces Γc-The above displacement boundary integration method derives the source point s to obtain a stress boundary integration equation (TBIEs) to solve the problem (as shown in equation 2 and equation 3). It can be seen that the integrand in the force boundary integral equation
Figure BDA0001740229700000091
And
Figure BDA0001740229700000092
since the derivation operation makes them have strong singularity at the source point s, Vij＝O(r-2) And super Singularity (Hyper Singularity), Wij＝O(r-3) Using Cauchy principal value integral sign respectively
Figure BDA0001740229700000093
And Hadamard Finite Part (HFP) integral sign
Figure BDA0001740229700000094
To indicate.
Figure BDA0001740229700000095
Figure BDA0001740229700000096
Γ in equation 2c+Subscript "+" of (a) and Γ in equation 3c-The subscript "-" is to distinguish the overlapping surfaces from both sides of the fracture. What appears to the left in equations 2 and 3 is called the Free term (Free Terms) that results from the derivation of the Somigliana equation, cijReflecting the geometry at the source point s.
Step (2) since the used boundary element method mainly aims at the field of computer graphics, it is necessary to perform Linear Discretization (Linear Discretization) on a segmented smooth surface model (as shown in fig. 3(a)) widely existing in the field of computer graphics. Assume that the model used is represented by C1A piecewise smooth surface of continuous triangular Patches (Patches) and subjected to linear fracture analysis. Since the segmented smooth surface is used as the boundary of the model, the source point s is a discontinuous (Non-smooth) point, and is shared by multiple triangular patches (as shown in fig. 3 (a)). Furthermore, since the integrand F (where F can be any of U, T, W, V) in the boundary integral equation has singularity on the distance of the field point x (the integration point) to the source point s, it is natural to divide the piecewise continuous surface Γ into two parts: singular facet set ΓsAnd a set of triangular patches ΓrAmong which a singular patch set ΓsIs composed of patches sharing a source point s (e.g. green part in fig. 3(a)), and the remaining gray part constitutes a triangular patch set Γr. When the boundary integral equation is solved on the triangular surface slice set, a normal numerical integration method (such as Gauss Quadrature Rule) can be adopted, and good calculation accuracy can be obtained. For the solution of the boundary integral equation of the singular patch set, the method needs to establish a culled Domain (Excluded Domain) at the source point, where the Excluded Domain is a sphere with a radius of epsilon and a center at the source point s (as shown in fig. 3(b)), and ρ ═ epsilon is the form under polar coordinates. The spherical culling domain is adopted to utilize the spherical culling domain in the free item cijAnd (4) the advantages in calculation. After the above processing, the singular surface patch set gamma issThe above equation 1 becomes the form of equation 4.
Figure BDA0001740229700000101
Equation 4 is the Theoretical Form of the Somigliana equation (Theoretical Form), and two modifications are required to convert equation 4 into a numerically usable Form. The first is to satisfy each of C1Consecutive arbitrary triangular patches are mapped from global Space (Mapping) to canonical domain in Intrinsic geometry Space (the invention uses right triangles) to be numerically integrated. To handle the singular integrand, we use Polar Coordinates (Polar Coordinates) to represent the parametric space. In FIG. 4, culling domains T for source points s and s whose vertices at non-smooth boundaries are shared by multiple patches are shownεState in coordinate transformation. In fig. 4 b, it can be seen that the origin is a map (Image) of the source point s in the geometric space in the polar coordinate system. It should be noted that the Spherical culling Domain (Spherical exposed Domain) at the source point s does not necessarily maintain a Spherical shape in the implicit geometric space after coordinate transformation (as shown in fig. 4 (b)). The integrand U of formula 4 is transformed by polar coordinatesij(s,x)uj(x) And Tij(s,x)tj(x) Become into
Figure BDA0001740229700000102
And
Figure BDA0001740229700000103
where N represents a shape function of the field point x in the implicit geometric space with respect to the source point s, J represents a Jacobian Transformation (Jacobian Transformation) from the global space to the implicit geometric space, ρ represents a Polar Radius (Polar Radius), and θ represents a Polar Angle (Polar Angle). The second is to change equation 4 into the boundary integral equation (see equation 5) in the case that multiple patches share the source point s, and to maintain the continuity of the mapping form between adjacent patches during the transition. In equation 5, N represents the number of shared patches.
Figure BDA0001740229700000104
And (3) obtaining a boundary integral equation defined on the segmented smooth surface through the processing of the step (2), and obtaining a formula 5, wherein the equation contains various singular integrand. To enable uniform processing of these singular integrand, we use Singular Subtraction Techniques (SST). The singularity reduction method firstly performs Taylor expansion on singular integrand under a polar coordinate system, so that singular terms of the integrand can be expressed Explicitly (explicit), and a Regular (Regular) integrand without the singular terms is obtained by subtracting the singular terms from the integrand. And finally solving the Analytic solutions (Analytic solutions) of the removed singular terms, and bringing the result back to the boundary integral equation. The limiting operation in equation 5 can be avoided by using the singularity reduction method. For convenience of description, the method is only applied to the Hyper-singular (Hyper-singular) integrand
Figure BDA0001740229700000111
The singularity reduction operation is performed, see equation 6. To use the singularity reduction method, the method assumes that the Displacement Vector (Displacement Vector) u satisfies when the field point x approaches the source point s
Figure BDA0001740229700000112
The continuous condition is recorded as u ∈ C0,α(s). Satisfied at the displacement vector u
Figure BDA0001740229700000113
On the basis of continuous conditions, the method can carry out singular integral terms in the formula 6
Figure BDA0001740229700000114
A series expansion operation is performed, thereby obtaining equation 7.
Figure BDA0001740229700000115
Figure BDA0001740229700000116
In the formula (5-7)
Figure BDA0001740229700000117
And
Figure BDA0001740229700000118
are two functions related only to Polar Angle (Polar Angle) θ. The process of performing the singularity reduction operation can be expressed as formula 8, and the Double Integral (Double Integral) equation at the right side in formula 8 has reduced (singular) singular terms
Figure BDA0001740229700000119
The limit operation and the integral equation at the outer side are decoupled, and the decoupling can be used as a normal integral equation to carry out numerical integration processing.
Figure BDA00017402297000001110
The Single layer Integral (Single Integral) on the right side of equation 8 is the portion of the singular term that needs to be solved analytically. The final equation 5 from step (2) can be processed into a form containing only conventional integral equations, see equation 9. For Strong Singular terms (Strong Singular Term) and Weak Singular terms (Weak Singular Term), they are processed in the same way as the super Singular terms, and are simpler to process.
Figure BDA0001740229700000122
And (4) providing a stable and efficient numerical solving method for solving the three-dimensional elastic fracture boundary element modeling problem on the segmented smooth surface faced by the method. From the numerical valueBoth the integration and the free term solution provide further numerical optimization to the existing methods. After the processing of step (3), the initial boundary integral equation containing the over-singularity product function, formula 1, has been converted into a regular numerical integral equation by a singularity subtraction method (SST), see formula 9. In equation 9, it can be seen that in the singularity reduction processing of the boundary integral equation containing the singularity integrand, polar coordinates are used to represent the reference space (see step (2)) to take advantage of the polar coordinates in processing the singularity integrand. It should be noted that, although the regular numerical integration method expressed by polar coordinates can also use a Standard Gaussian numerical integration method (Standard Gaussian Quadrature components) to solve the method, the numerical calculation accuracy needs to be satisfied by increasing the sampling density of the integration points. Furthermore, for the case that the boundary Element (Triangle Element used in the present invention) has a large Aspect Ratio (Aspect Ratio), the integrand corresponding to the Element still has a singular attribute. This occurs because in the Laurent series expansion form of the integrand obtained by equation 7, the function
Figure BDA0001740229700000123
Although decoupled from the polar axis variable ρ, singularities from the polar angle θ still appear during the integration operation for the polar angle θ, since it is still a function of the polar angle θ. In view of the above problem, the numerical solution stability of the existing boundary integral equation is improved by transforming the function projection of the polar angle θ obtained by equation 7 into a stable conformal space.
For boundary elements with larger aspect ratios, the resulting function in equation 7 has the following expression:
Figure BDA0001740229700000131
wherein V at the moleculei(θ) is a function of θ, and p is a constant determined by the index i. The function A (θ) is a reflectionWith the function of the boundary element geometry information expanded as shown in equation 11, it can be seen that in equation 11
Figure BDA0001740229700000132
Global to local coordinates (ξ) representing boundary elements1,ξ2) The basis vectors of the upper map. We measure the projection space by defining the ratio λ of the lengths of the two space basis vectors and the cosine cos (γ) of the included angle γ of the two space basis vectors, as detailed in equation 13. From equation 13, two auxiliary coefficients can be obtained
Figure BDA0001740229700000133
And μ, see equation 14. Using the double angle formula for formula 11 and substituting formula 13, a new expression for function a (θ), formula 12, can be obtained. It can be seen that the coefficient θ can be determined by making μ approach 1
Figure BDA0001740229700000134
Results in a function A (θ) → 0, and thus in a function
Figure BDA0001740229700000135
Possesses singular properties. From a geometric point of view, there are two factors that contribute to the above-mentioned situation (μ → 1): (1) due to projection of the spatial basis vector u1And u2They are caused by either a close overlap (angle γ → 0) or a parallel (angle γ → π); (2) due to | u1I and I u2The ratio of | tends to 0 or infinity, i.e., the aspect ratio is too large, resulting in (λ + λ)-1)2→∞。
Figure BDA0001740229700000136
Figure BDA0001740229700000137
Figure BDA0001740229700000138
Figure BDA0001740229700000139
The method for solving the problems is to perform Conformal Transformation (Conformal Transformation) again on the basis of the projection of the existing global coordinate system (fig. 7) to the local coordinate system (fig. 7b), as shown in fig. 7c, and make a new projection space basis vector through Conformal Transformation
Figure BDA0001740229700000141
The following conditions are satisfied:
Figure BDA0001740229700000142
when the new basis vectors are substituted into equation 12, the function A (θ) becomes a constant term from the previous function with respect to the polar angle θ
Figure BDA0001740229700000143
The advantages that result from this are: (1) since A (θ) is a constant term, in equation 10
Figure BDA0001740229700000144
The singularity for the polar angle θ has been decoupled; (2) a simple equation 10 reduces the complexity of numerical integration along the polar angle theta.
Figure BDA0001740229700000145
Figure BDA0001740229700000146
To achieve this goal, two auxiliary coefficients need to be calculated
Figure BDA0001740229700000147
And
Figure BDA0001740229700000148
(see equation 16) to represent the coordinates of the new point in projection space, as shown in FIG. 7(c), where λ and cos (γ) are defined by equations (5-23). Then, a new projection matrix T (see formula 17) can be obtained, and a new local coordinate (see fig. 7c) can be obtained by projecting the existing local coordinate system (see fig. 7b) using the projection matrix T, and it can be seen that the basis vector of the new conformal space
Figure BDA0001740229700000149
And
Figure BDA00017402297000001410
equation 15 is satisfied. Since a new conformal space projection is performed, equation 9 needs to be modified to obtain equation 18. | T in the formula-1| is defined by equation 17, the polar angle θ becomes the polar angle in the new projection space.
Figure BDA0001740229700000151
For the modeling of the three-dimensional fracture problem on the segmented smooth surface related by the method, a free term explicit computing technology is adopted to support explicit computation of a free term coefficient matrix at the boundary vertex shared by any multiple tangent planes on the segmented smooth surface.
For the sake of description, the calculation process of the coefficient matrix C is described below by way of an example, which can be applied to the general case of piecewise smooth surfaces, including the case of overlapping fracture surfaces. As shown in FIG. 8, assume in the example that the current source point s is bounded by five boundary elements (satisfy C)1Arbitrary triangular patches of contiguous conditions). Because the spherical culling domain is adopted (see step (2)), the tangent plane set pi where the 5 boundary elements are located and the spherical culling domain are intersected to obtain a vertex set v ═ { v } located on the spherical surface1,v2,v3,v4,v5The sum of the contour arcs is gamma-gamma1,2,γ2,3,γ3,4,γ4,5,γ5,1}. Here, the vertices and contour arcs are ordered clockwise, as shown in FIG. 8. For the purpose of facet description, the following set of variables are defined: (1) defining a tangent plane set pi ═ pi1,2,π2,3,π3,4,π4,5,π5,1}; (2) defining the set of external normals for the set of planes pi as n ═ n1,2,n2,3,n3,4,n4,5,n5,1}; (3) the radian set of included angles of two intersecting planes is defined as alpha ═ alpha1,α2,α3,α4,α5For α }iSee equation 19 for the calculation of:
αi＝π+sgn((ni-1,i×ni,i+1)·ri)arccos(ni-1,i·ni,i+1) (19)
in the formula 19, ri＝r(vi)＝viS, s denotes the source point. sgn (x) is a sign function.
To this end, we can obtain a general formula 20 of the coefficient matrix C corresponding to the local geometry formed by the intersection of n tangent planes. Using equation 20, we obtain a 3 × 3 matrix C ═ C for the source point sijWhere i, j ═ 1,2, 3.
Figure BDA0001740229700000152
And (5) in order to better accurately model the generation process of the fracture, the method introduces a continuous contact model, and places the fracture in a continuous contact time period to represent the dynamic process of the external force acting on the surface of the model. And during the contact time, the change of the force applied in the contact process is represented by defining the external force which changes along with the time. To establish the continuous contact period, the present invention employs a Hertz model of elastic ball contact, see equation 21.
Figure BDA0001740229700000161
Where m represents mass, E represents elastic modulus, R represents radius of the contact ball, v represents velocity of the model when in contact with the outside, and c is constant coefficient of the Hertz model (parameter used in the present invention, c ═ 2.87). The breaking process can be regarded as a simulated operation within a time period by the user-specified simulated time step Δ t and the continuous contact time period obtained by equation 21.
Since the Hertz Contact Model is a Non-adhesive Contact Model, the method regards the external Contact force as a sine Function (sinussoidal Function), see equation 22.
Figure BDA0001740229700000162
Wherein f ismaxRepresenting the maximum force to which the surface of the model is subjected during contact.
Step (6) is directed to crack propagation, and Stress Intensity Factors (SIFs) are used to describe the Stress field near the crack because the Stress has singularity at the crack tip. And performing deformation simulation on the model under the action of continuous external force to obtain the Displacement of boundary elements on the surface of the Crack, thereby obtaining the Crack Opening Displacement (COD), and calculating a Stress Intensity Factor (SIF) at the Crack according to the Crack Opening Displacement. It follows that crack propagation is a process of propagation based on the crack tip, according to the stress intensity factor and the material properties of the model.
Calculating the stress intensity factor requires establishing a local coordinate system for the extension point at the leading edge of the crack. For the local coordinate system of the extension point c, firstly, the crack ct of the extension point c and the boundary element elt of the crack are determinedctTo obtain the boundary element eltctNormal n to1And a tangent n parallel to the crack ct3Through n1And n3The Cross Product (Cross Product) yields a normal n perpendicular to the crack ct2As shown in FIG. 5. In FIG. 5, n1，n2，n3Are three mutually perpendicular local basis vectors, where n1Perpendicular to the boundary element eltctForce applied in the opening direction (see Modei in FIG. 6), n2Perpendicular to the crack ct, corresponding to the shear direction (see ModeII in FIG. 6), n3Parallel to the crack ct, corresponds to the force in the sliding direction (see ModeIII in fig. 6).
Figure BDA0001740229700000171
Figure BDA0001740229700000172
Figure BDA0001740229700000173
For Δ u may pass the boundary element eltctIs obtained from the vertex p not belonging to the crack ct. After the opening displacement delta u is obtained, the opening displacement delta u needs to be projected onto a local coordinate base vector to obtain delta uI，ΔuIIAnd Δ uIIIThen, the stress intensity factor K is calculated according to equation 13I，KIIAnd KIIIIn formula 23, μ represents Shear Modulus (Shear module), μ ═ E/(2+2v), E represents young's Modulus, v represents poisson's ratio, and r represents the distance from the vertex p to the crack ct. When the propagation point is shared by a plurality of cracks, the final stress intensity factor may be obtained by averaging the stress intensity factors in the plurality of cracks.
Step (7) modeling of the object fracture process requires handling of Crack Initiation and Propagation. For the boundary element method, the generation of the fracture is based on the analysis of the surface Stress, and the maximum Principal Stress (maximum Principal Stress) of the boundary element (triangular patch) is calculated, which causes the fracture phenomenon when the maximum Principal Stress exceeds the material strength. And calculating the principal stress and the principal stress direction of all boundary elements owned by the model one by one, then checking whether the current boundary elements need to generate cracks, initializing the crack by taking the current principal stress direction as the normal of the surface of the crack when the crack generation condition is met, and extending the crack on the basis of the cracks. The crack formation ends when no new crack surface is produced.
In the crack initialization phase, in order to perform surface stress analysis of the boundary element, firstly, a Deformation Gradient (Deformation Gradient) of the current element is obtained, and considering that the boundary element adopted by the invention is a space triangle, the states of three vertexes of the boundary element before and after Deformation cannot completely determine affine transformation of the space triangle (because of the lack of dimension perpendicular to the space triangle). To solve this problem, the present invention adds its vertical points to the space triangle, see equation 24. Let v be the three pre-deformation vertices of a spatial triangleiI is 1,2,3, and the vertexes after three deformations are
Figure BDA0001740229700000181
V in equation 244That is, adding a new point perpendicular to the spatial triangular patch, the transformed point can be obtained by the same method as in equation 24
Figure BDA0001740229700000182
After the points in the vertical direction are obtained, an edge vector matrix V ═ V of the space triangular patch can be constructed2-v1 v3-v1 v4-v1]And
Figure BDA0001740229700000183
using equation 24 to obtain the deformation gradient F3×3. Obtaining F by performing Singular Value Decomposition (SVD) on the deformation gradient F3×3＝U∑VTSubstituting the singular value matrix sigma into equation 26 calculates the First Piola-Kirchhoff stress tensor P for the triangular patch, where P represents the stress per unit area, and the direction of the stress corresponds to the column vector of V.
Figure BDA0001740229700000184
Figure BDA0001740229700000185
P＝2μ(∑-I)+λTr(∑-I)I (26)
When the maximum stress in the stress tensor P exceeds the material strength, a new split needs to be created on the current triangular face. The centroid of the current triangular patch is first located and a new break is inserted with the direction of maximum stress (column vector from V) as normal to the break surface.
When the crack is initiated, the crack may be extended during the simulation. First, the stress intensity factor K is calculated for the point on the crack using equation 23I，KIIAnd KIII. Then, K is addedI， KIIAnd KIIISubstituting equation 27 to obtain the effective stress intensity Keff(Effective Stress Intensity). By applying an effective stress intensity KeffFracture toughness K of materialcFor comparison, when KeffGreater than KcThe crack is propagated.
Figure BDA0001740229700000186
Figure BDA0001740229700000187
Figure BDA0001740229700000191
Figure BDA0001740229700000192
During the propagation of the crack, the propagation needs to be calculatedVelocity vcpAnd the angle theta of extensioncpThis can be obtained by equation 28 and equation 30. In equation 28 cRThe speed of the Rayleigh wave is represented,
Figure BDA0001740229700000193
ρ is the density of the material. To calculate the extension angle thetacpFirst, K is requiredIAnd KIIIAre combined by equation 29, since KIAnd KIIIPredominate the forward propagation of cracks, KIIResulting in lateral rotational propagation of the crack. Will KI,IIIAnd KIISubstituting equation 30 yields θcpFinally, the extension angle theta needs to be extendedcpAnd converted to a direction under the global coordinate.
The software platforms used for realizing the prototype system are Microsoft visual studio 2010 and OpenGL, and the hardware platforms are Intel Core i73.60GHz CPU, NVIDIA GeForce GTX 770GPU and 4G memory. The method effect is shown in fig. 9 and fig. 10.
The above description is only a preferred embodiment of the present invention, and not intended to limit the present invention in other forms, and any person skilled in the art may apply the above modifications or changes to the equivalent embodiments with equivalent changes, without departing from the technical spirit of the present invention, and any simple modification, equivalent change and change made to the above embodiments according to the technical spirit of the present invention still belong to the protection scope of the technical spirit of the present invention.

Claims (5)

1. An elastic fracture simulation method based on dual boundary elements and strain energy optimization analysis is characterized in that: the method comprises the following steps:
step (1), carrying out dual boundary integral equation modeling on a target model and cracks thereof;
step (2), dispersing a boundary integral equation based on a segmented smooth surface;
step (3), singular integrand functions in the boundary integral equation are processed uniformly by a singularity reduction technology;
step (4), solving a free term in a boundary integral equation;
step (5), in the simulation process, accurately modeling the fracture generation process by using a continuous contact model;
establishing a continuous contact time period by adopting a Hertz model contacted by an elastic ball, representing the change of stress in the contact process by defining an external force changing along with time in the contact time period, and representing the dynamic process of the external force acting on the surface of the model by placing a fracture in a continuous contact time period through a simulation time step specified by a user and the obtained continuous contact time period;
step (6), in the simulation process, describing a stress field near the crack by using a stress intensity factor;
carrying out deformation simulation on the model under the action of continuous external force to obtain the displacement of boundary elements on the surface of the crack so as to obtain the opening displacement of the crack, establishing a local coordinate system for an extension point at the front edge of the crack, projecting the opening displacement onto a local coordinate basis vector, calculating stress intensity factors, and averaging the stress intensity factors in a plurality of cracks to obtain a final stress intensity factor when the extension point is shared by the plurality of cracks;
step (7), in the process of generating fracture, initializing fracture according to surface stress analysis, and extending on the crack based on a stress intensity factor;
aiming at the surface stress analysis of boundary elements, firstly obtaining the deformation gradient of the current element, obtaining the First Piola-Kirchhoff stress tensor of the triangular patch, and generating a new crack on the current triangular patch when the maximum stress exceeds the strength of a material; and calculating the stress intensity factor of the point on the crack, extending the crack when the effective stress intensity is larger than the fracture toughness of the material, and calculating the extending speed and the extending angle in the extending process of the crack.
2. The elastic fracture simulation method based on dual boundary elements and strain energy optimization analysis according to claim 1, wherein the step (1) comprises: establishing a boundary integral equation for the region omega of the target model by adopting a Betfi-Somigliana equation; the boundary integral equation at the source point is expressed by a Kelvin base solution formed by the distance between the source point and the field point of the domain boundary; and modeling the cracks in the domain into an open surface, modeling the three groups of surfaces into independent displacement boundary integral equations by using dual boundary elements, and deriving the source points by using a displacement boundary integral method on the overlapped surfaces to obtain a stress boundary integral equation.
3. The elastic fracture simulation method based on dual boundary elements and strain energy optimization analysis according to claim 1, wherein the step (2) comprises: performing linear dispersion on a segmented smooth surface model widely existing in the field of computer graphics; solving the boundary integral equation on the triangular patch set by adopting a normal numerical integration method; for solving the boundary integral equation of the set of the surface patches with numerical singularity, a culling domain needs to be established at a source point.
4. The elastic fracture simulation method based on dual boundary elements and strain energy optimization analysis according to claim 1, wherein the step (3) comprises: unified processing is carried out on singular integrand functions in the boundary integral equation by adopting a singularity reduction technology; taylor expansion is carried out on singular integrands under a polar coordinate system, so that singular terms of the integrands are explicitly expressed, and a normal integrand without the singular terms is obtained by subtracting the singular terms from the integrands; their analytical solutions are solved for the removed singular terms and the results are brought back to the boundary integral equation.
5. The elastic fracture simulation method based on dual boundary elements and strain energy optimization analysis according to claim 1, wherein the step (4) comprises: adopting a spherical elimination domain, and performing intersection operation on a tangent plane set where boundary elements around a source point are located and the spherical elimination domain to obtain a vertex set and a contour arc set which are located on a spherical surface; the following set of variables is defined: a tangent plane set, an outer normal set of the tangent plane set and an included angle radian set of two intersected tangent planes; and obtaining a coefficient matrix related to the set structure of the local surface where the source point is located.
CN201810815571.9A 2018-07-24 2018-07-24 Elastic fracture simulation method based on dual boundary element and strain energy optimization analysis Active CN108763841B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810815571.9A CN108763841B (en) 2018-07-24 2018-07-24 Elastic fracture simulation method based on dual boundary element and strain energy optimization analysis

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810815571.9A CN108763841B (en) 2018-07-24 2018-07-24 Elastic fracture simulation method based on dual boundary element and strain energy optimization analysis

Publications (2)

Publication Number Publication Date
CN108763841A CN108763841A (en) 2018-11-06
CN108763841B true CN108763841B (en) 2022-04-05

Family

ID=63971288

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810815571.9A Active CN108763841B (en) 2018-07-24 2018-07-24 Elastic fracture simulation method based on dual boundary element and strain energy optimization analysis

Country Status (1)

Country Link
CN (1) CN108763841B (en)

Families Citing this family (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110705057B (en) * 2019-09-19 2021-05-18 武汉大学 Method and device for solving static thermoelasticity problem of isotropic solid material
CN112784495B (en) * 2021-01-28 2021-09-24 郑州轻工业大学 Mechanical structure real-time fatigue life prediction method based on data driving
CN113158531B (en) * 2021-02-07 2022-06-21 南开大学 Single-component and multi-component incompressible fluid simulation method utilizing deformation gradient
CN112926213B (en) * 2021-03-11 2024-01-16 郑州轻工业大学 Method, system, medium, equipment and terminal for measuring thermal damage boundary element
CN114511534B (en) * 2022-01-28 2023-05-05 江苏泰和木业有限公司 PC board crack judging method and system based on image processing
CN114818197B (en) * 2022-05-10 2024-04-12 西安交通大学 High-speed motorized spindle thermoelastic deformation simulation method and system based on boundary element model

Family Cites Families (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7509245B2 (en) * 1999-04-29 2009-03-24 Schlumberger Technology Corporation Method system and program storage device for simulating a multilayer reservoir and partially active elements in a hydraulic fracturing simulator
US7505885B2 (en) * 2003-01-24 2009-03-17 The Boeing Company Method and interface elements for finite-element fracture analysis
US8494827B2 (en) * 2009-09-25 2013-07-23 Exxonmobil Upstream Research Company Method of predicting natural fractures and damage in a subsurface region
CN101788425A (en) * 2010-02-09 2010-07-28 浙江工业大学 Method for determining separation and distribution of structural-member composite crack front stress intensity factors
CN102332046B (en) * 2011-09-30 2012-12-19 北京工业大学 Gear crack propagation simulated wavelet extension finite element simulation analysis method
CN103020426B (en) * 2012-11-23 2015-11-18 北京航空航天大学 A kind of short-cut method of rectangular slab Plate with Inclined Center Crack propagation life of fatigue prediction
CN105302974B (en) * 2015-11-06 2018-06-08 北京航空航天大学 A kind of real-time cutting simulation method of flexible article analyzed based on finite element and time-varying modal

Also Published As

Publication number Publication date
CN108763841A (en) 2018-11-06

Similar Documents

Publication Publication Date Title
CN108763841B (en) Elastic fracture simulation method based on dual boundary element and strain energy optimization analysis
Kineri et al. B-spline surface fitting by iterative geometric interpolation/approximation algorithms
Dietrich et al. Edge transformations for improving mesh quality of marching cubes
WO2021203711A1 (en) Isogeometric analysis method employing geometric reconstruction model
Farin et al. Geometric Modeling
Yoshihara et al. Topologically robust B-spline surface reconstruction from point clouds using level set methods and iterative geometric fitting algorithms
Wang et al. Hole filling of triangular mesh segments using systematic grey prediction
Aubry et al. A three-dimensional parametric mesher with surface boundary-layer capability
Aubry et al. A robust conforming NURBS tessellation for industrial applications based on a mesh generation approach
Leng et al. Parameterized modeling and optimization of reusable launch vehicles based on reverse design approach
Abdul-Rahman et al. Freeform texture representation and characterisation based on triangular mesh projection techniques
Yang et al. An open source, geometry kernel based high-order element mesh generation tool
Douglass et al. Current views on grid generation: summaries of a panel discussion
Scheirer et al. Fiber layup generation on curved composite structures
Shojaee et al. Combination of isogeometric analysis and extended finite element in linear crack analysis
Ren et al. A Review of T-spline Surfaces
Lin et al. A new interpolation subdivision scheme for triangle/quad mesh
Bock et al. Generation of high-order polynomial patches from scattered data
Ba et al. Research on 3D medial axis transform via the saddle point programming method
Yu et al. Smooth geometry generation in additive manufacturing file format: problem study and new formulation
Zhao et al. Error analysis of reparametrization based approaches for curve offsetting
Lin et al. Fusion of disconnected mesh components with branching shapes
Xie et al. A triangulation-based hole patching method using differential evolution
Jiang et al. G 2-continuity extension algorithm of ball B-Spline curves
Kang et al. Application of morphing technique with mesh-merging in rapid hull form generation

Legal Events

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