Continuous fiber reinforced ceramic composite cross-scale damage heterogeneous parallel method

Abstract

The invention discloses a continuous fiber reinforced ceramic composite cross-scale damage heterogeneous parallel method, which comprises the following steps of: establishing a macro-scale finite element model, and establishing a micro-fiber bundle scale finite element model through any weaving mode; calculating Gaussian point strain of all units under a macroscale; calculating displacement boundaries of all units under the macro Gaussian point strain under the micro scale; calculating phase field damage of all the units under the microscopic scale, and calculating displacement of all the microscopic units after damage through the phase field damage; calculating the average stress and the average modulus after damage of all the units under the micro-scale; internal forces of all units under the macro scale are calculated through the average stress of the micro units, and the process is iterated by using the internal forces and external forces, so that damage evolution is predicted.

Description

Technical Field

[0001] The present invention belongs to the technical field of high-performance finite element simulation, and specifically is a parallel method for cross-scale damage heterogeneity of continuous fiber reinforced ceramic composites. Technical Background

[0002] Composite materials have received extensive research and attention due to their superior performance. Compared to traditional materials, composite materials have higher strength, higher stiffness, and lower density. Therefore, composite materials are widely used in high-end equipment such as aerospace and medical devices.

[0003] Composite components are exposed to harsh service environments for long periods of time. These components often lose their mechanical properties due to failures caused by various reasons. Compared to these explicit failures that are easy to observe and repair, the damage and destruction generated within the composite material is not only difficult to observe but also has more serious consequences. The cause of internal damage is that composite materials have a complex microstructure. For example, woven composite materials are usually composed of fibers, which are then combined with a matrix to form a preform, and finally form a macroscopic equipment component. Due to the imperfect composite material preparation process, defects are inevitably introduced during forming. These defects include pores and cracks, and are difficult to observe with the naked eye. When large equipment components are in service, they are affected by long-term loads. After reaching the stress limit, the cracks and pores inside the components will begin to expand and evolve.

[0004] Therefore, a cross-scale damage heterogeneous parallel method for continuous fiber reinforced ceramic composites is established, and the mature finite element numerical calculation theory is applied. The advantages of different parallel methods are fully utilized to accelerate the calculation, which can more effectively utilize computing resources to achieve rapid solutions to large-scale models. At the same time, based on the cross-scale approach, the microstructure of the composite material can be better reflected, and the damage calculation of the composite material can be realized by using the phase field method, and the crack propagation of the composite material can be predicted. This method can not only be applied to fields such as aerospace, but can also provide more accurate solutions in other aspects, reduce the time cost required to predict the performance of composite materials, and efficiently feedback to the design, accelerating product iteration. Summary of the Invention

[0005] The present invention discloses a cross-scale damage heterogeneous parallel method for continuous fiber reinforced ceramic composites, belonging to the field of finite element high-performance computing simulation technology. The method performs damage calculation based on the phase field damage calculation method and further analyzes the microstructure of the composite material across scales. It also combines the parallel advantages of clusters, CPUs, and GPUs to provide an efficient and high-precision solution method, thereby realizing the rapid solution of large-scale composite material damage.

[0006] The application provides a continuous fiber reinforced ceramic composite trans-scale damage heterogeneous parallel method, which comprises the following steps:

[0007] S1, a macro-scale component finite element model is established, and a continuous fiber reinforced ceramic matrix composite material model in an arbitrary weaving mode under a micro-scale is established by using a B-spline curve mode;

[0008] S2, CPU parallel calculation is introduced, and the displacement and strain of all units of the macro model are calculated by using a finite element method;

[0009] S3, cluster parallel calculation is introduced, and the displacement boundary of the micro unit on the Gaussian point of the macro unit is calculated through the strain of the Gaussian point;

[0010] S4, according to the displacement boundary, the phase field value, the displacement field value and the average stress value after damage are calculated through a phase field damage mode under the micro-scale;

[0011] S5, the average modulus of each unit after damage under the micro-scale is calculated through a homogenization method;

[0012] S6, the internal force of the macro unit is calculated by using the average stress of the micro unit under the macro-scale;

[0013] S7, the stiffness matrix of the macro unit after damage is calculated by using the average modulus of each macro unit after damage;

[0014] S8, the external force obtained by using the total stiffness matrix after damage, the internal force and the applied load is used for iteration, so that the damage evolution is predicted;

[0015] S9, GPU parallel calculation is introduced, and the stress and strain of the macro model at this time are calculated according to the displacement value to realize post-processing.

[0016] The specific implementation method of the above step S1 is as follows:

[0017] S11, a finite element model under a macro-scale is established;

[0018] The macro model is regarded as a homogeneous material, and the material parameters are obtained by using a homogenization method to predict the modulus of the undamaged micro model. According to the results, the Poisson's ratio and the elastic modulus are set, and the boundary conditions and the constraints need to be set, and the positions of the displacement boundary and the fixed constraint are determined.

[0019] S12, a micro model of the continuous fiber reinforced ceramic matrix composite material in an arbitrary weaving mode is established;

[0020] The parameterized modeling is established by using a level set method and a non-uniform rational B-spline curve to establish the composite material micro element in the hot weaving mode.

[0021] The B-spline curve describes the fiber position of any weaving method, and the shortest distance between the node and the B-spline curve is calculated using the grid method and gradient descent method. It generates a curve by generating multiple rational polynomials through control points, and its C(t) is analytically defined as:

[0022]

[0023] Through the control points, a B-spline curve with arbitrary spatial distribution can be generated, and the curve can be used as the axis of the fiber. i :i=0,1,...,n} as control points, N i,j (t) is the basis function of the j-th order i B-spline curve described by m non-decreasing knot sequences t, which are t i The j-order piecewise polynomial of control is expressed as:

[0024]

[0025] The B-spline curve with any form of spatial distribution can be generated by the control points. The curve is used as the axis of the fiber. In order to determine whether the unit passes through the interface layer between the fiber and the matrix, it is necessary to obtain the P at any point in space. a The shortest distance d from (a, b, c) to the fiber curve, d is about point P on the fiber curve b The function of (x,y,z) can be expressed as:

[0026]

[0027] The above maximum value problem can be equivalently solved to find the minimum value of its square, which can be written as an equivalent distance function, which is:

[0028]

[0029] The grid method combined with the gradient descent method is used to solve the point P on the fiber axis. b (x, y, z) is a function of parameter t, so the equivalent distance function obtains the gradient of the fiber curve with respect to parameter t

[0030]

[0031] Find the approximate parameter t by grid point method rough , then use t rough Initialize the parameter t in the gradient descent method, use the gradient descent method to continuously approximate the value of t, and continuously modify the parameter t in the iteration until the point P on the fiber curve corresponding to t is bmin (x,y,z) and P a Shortest distance:

[0032]

[0033] The distribution of fibers and matrix is determined by using a level set method, so as to divide the fibers and the matrix. For a description of one or more fibers in a cubic model discretized by tetrahedral elements, the position of the fiber surface can be described by a level set function of the fiber. In the three-dimensional case, the level set method represents the fiber surface curve ζ as the zero level set of a three-dimensional auxiliary function φ:

[0034] ζ = {(x, y, z) | φ(x, y, z) = 0}

[0035] In order to represent the shape region Ω of the fiber, a level set function φ is defined in three-dimensional space, and the zero level set of φ represents the boundary Γ between the fiber and the matrix. When φ is less than zero, it represents the inside of the fiber region, and when φ is greater than zero, it represents the outside of the fiber region. The structure of the fiber also needs to define the cross-sectional shape of the fiber. First, the position of the cross section needs to be determined, and the first derivative of the expression of the fiber axis with respect to t is taken to obtain the tangent vector expression of the fiber axis

[0036]

[0037]

[0038] A point o(t) is taken on the fiber axis, and the tangent vector at this point is taken as the x` axis at this point. The global z-axis direction is taken as the z` axis direction at this point, and the y` axis of the elliptical cross section is calculated using the cross product. Then the y-axis vector is normalized, and finally the z` axis normalized vector at this point is obtained by cross product of x-axis and y-axis. The level set function value of the point to the fiber interface is calculated by the equation, and the position relationship between any point on the plane and the ellipse, rectangle is obtained.

[0039] Thus, the finite element model of the macro-scale component is completed, and a B-spline curve method is used to establish a model of a continuous fiber reinforced ceramic matrix composite material with any weaving method at the micro-scale;

[0040] The specific implementation method of the above step S2 is: after dividing a continuous domain into multiple discrete domains, that is, dividing a complete model into elements and nodes in the finite element method. For any discrete domain Ω in space, the following three basic equations need to be satisfied, that is, the balance equation, the geometric equation and the physical equation.

[0041] σ ij,j +F i = 0 (i = 1, 2, 3; j = 1, 2, 3)

[0042] Where, σ ijFor the second order symmetric stress tensor, the equilibrium equation establishes the physical relationship between the principal stress and shear stress in three directions of the elastic body in space

[0043]

[0044] where u, v, w represent the displacement in x, y, z directions of the Cartesian coordinate system respectively. By taking the partial derivative of the displacement in the corresponding direction, the geometric equation establishes the physical relationship between the strain component and the displacement component.

[0045] σ ij = D ijkl ε kl (i, j, k, l = 1, 2, 3)

[0046] where σ ij and ε kl are two-dimensional stress and strain tensors respectively. D ijkl is the elastic four-dimensional tensor in the domain Ω. By the symmetric four-dimensional elastic tensor, the physical equation establishes the physical relationship between the strain component and the stress component.

[0047] So far, by introducing three basic equations, there are a total of 15 unknowns for the elastic body in any three-dimensional space, including six directions of stress, six directions of strain and three directions of displacement. In order to solve the above equations, it is necessary to introduce force boundary conditions and displacement boundary conditions to specify the values and normal derivatives on the boundary.

[0048]

[0049] where n is the direction cosine of the normal of the boundary of the domain Ω, is the area force tensor on the boundary. The Dirichlet boundary condition is to directly specify the function value on the boundary, and thus the basic equations of elastic mechanics are introduced. Based on the equilibrium equation, the strong form of the differential equation can be derived and written in the form of integral:

[0050]

[0051] For the finite element calculation of the real physical model domain, it is necessary to discretize the domain into a physical model domain with finite degrees of freedom. And for each discretized physical domain, the relationship between the local displacement and the global displacement needs to be established, and the physical equation and the geometric equation need to be satisfied:

[0052] u = NU e

[0053] ε = BU e

[0054] σ = Dε

[0055] e represents each discretized physical domain, and Ue is the displacement of each node on a single element. B is the element geometry matrix, N is the element shape function matrix, and the above integral is converted into a summation method:

[0056]

[0057] Using the variational principle to take the minimum value of the potential energy variation to simplify the above form, the unit stiffness matrix K is obtained e :

[0058]

[0059] K e =∫ Ω B T DBdΩ e

[0060]

[0061] Each time the stiffness matrix of the microscopic model is assembled, due to the occurrence of damage, the stiffness matrix at this time needs to be updated using the following formula:

[0062] K=∫ Ω (1-d) 2 B T DBdΩ

[0063] Among them, the damage parameter d is the average value of the damage value of each node on the unit.

[0064] The stiffness matrix of each unit is only related to its own node coordinates and the material properties of the unit, and there is no data association between the units. Since the serial program of the unit stiffness matrix is ​​parallelized in a loop, OpenMP can be used to parallelize the serial program to calculate the stiffness matrix of each unit.

[0065] At this point, the introduction of CPU parallel computing is completed, and the finite element method is used to calculate the displacement and strain of all units in the macro model;

[0066] The specific implementation method of the above step S3 is:

[0067] S31, distribute macro unit processes through MPI;

[0068] After the macro-model calculation is completed, the micro-model calculation is only related to the strain of the macro-model unit in which the micro-model is located. There is no data exchange between the micro-models of different macro-units. Therefore, MPI parallel technology is used to split the macro-units according to the number of processes. Each process calculates a part of the unit of the complete macro-model. In order to ensure that the macro-units allocated to each process are as equal as possible, it is prevented that some processes have to wait for other processes after completing their calculations.

[0069] S32, calculating micro-unit model displacement boundary using Gaussian point strain;

[0070] For any unit, the actual calculation is the stress and strain at the Gaussian point of the unit, and the complete integral value is calculated by weighting each Gaussian point, and then the stress and strain values at the remaining positions of the unit are obtained by interpolation. That is, for any Gaussian point on the unit, the strain tensor at the Gaussian point is:

[0071]

[0072] The micro-scale model is regarded as a point of the macro model, and the strain at the point is used as the complete strain of the micro-scale model for calculation. The micro model strain is the product of the coordinates of all boundary nodes of the micro model and the strain. For any node on the boundary, there is:

[0073]

[0074] By this, the CPU parallel computing is completed to distribute each macro model unit, and the displacement boundary of the micro unit at the Gaussian point of the macro unit is calculated by the strain at the Gaussian point;

[0075] The specific implementation method of the above step S4 is:

[0076] S41, phase field calculation;

[0077] The phase field damage in the composite micro model is realized by using the phase field method, and then fed back to the material properties of the macro. The core idea is to introduce a continuous phase field variable d(x), the value range of which is [0, 1], which represents the state of the material from no damage to complete damage, where d=1 represents complete damage of the material, and d=0 represents no damage of the material.

[0078] The damage field in one dimension can be described by the following formula:

[0079]

[0080] Where l∈R is a scale parameter controlling the crack, not the actual length of the crack. When using finite element method to calculate the phase field, the scale parameter l in the crack density function is usually related to the grid size h, and after derivation and simplification, the expression is:

[0081]

[0082] d(x)-l 2 d″(x)=0

[0083] Since the differential equation is subject to Dirichlet boundary conditions and the differential equation can be written as the Euler control equation of the following formula based on the variational principle,

[0084]

[0085] Where W Γ is the Dirichlet boundary condition, that is, W Γ ={d|d(0)=1; d(x≠0)=0}; Arg represents the value of the independent variable when the formula in the brackets reaches the minimum value.

[0086] For a two-dimensional finite domain containing a crack, the expression for the phase field value d(x) is:

[0087]

[0088] in, represents the gradient operator, and vector n is the boundary normal. Since the two-dimensional case still needs to satisfy the value that minimizes the crack surface density function, the phase field value at this time also satisfies the Euler control equation. At this time, the Dirichlet boundary can be expressed as W Γ ={d|d(x)=1on x∈Γ}. In addition, the unit volume crack density function can be expressed as,

[0089]

[0090] The functional of the internal energy of the damaged body can be expressed as:

[0091] Π(u,Γ)=Π d (u,Γ)+Π s (Γ)

[0092] Among them, Π(u,Γ) is the total energy inside the object, Π d (u,Γ) is the elastic potential energy of the cracked elastic body, Π s (Γ) is the fracture energy required to generate a new crack. And the fracture elastic energy is the integral of the unit crack density function in the domain, that is,

[0093]

[0094] Among them, G c is the critical energy release rate. d (u,Γ) is usually expressed as the integral of the strain energy density function inside the elastic body within the domain, and its expression is:

[0095] Π d (u,Γ)=∫ Ω ψ(ε(u),d(x))dΩ

[0096] Where ψ(ε(u),d(x)) is the elastic energy density function, ε(u) is the strain tensor, and u is the displacement vector. When the elastic body is intact, that is, without damage, the elastic strain energy density is expressed as:

[0097]

[0098] Where C represents the fourth-order elasticity tensor.

[0099] Considering that crack generation and expansion mainly come from tensile loads rather than compressive loads, energy dissipation only occurs under tension. Based on this, the free energy density function is decomposed into tension and compression, and only the energy dissipation under the tensile model is considered. Its expression is:

[0100] ψ(ε)=[g(d)+k]ψ + (ε)+ψ - (ε)

[0101] Among them, ψ + / - (ε) represents the strain energy density function in tension and compression mode, respectively. g(d) is the fracture toughness function related to the phase field value.

[0102] S42, displacement field calculation and mean stress calculation;

[0103] After the displacement field calculation is obtained based on the previous displacement boundary calculation of the micro-model, the new displacement of the micro-model at this time will be obtained after the damage calculation is completed. This displacement is used to calculate the average stress of the micro-model at this time:

[0104]

[0105] After the mandatory boundary conditions are applied, the finite element method and homogenization theory are used to calculate the mean stress under each mandatory periodic boundary condition:

[0106]

[0107] S43, process function calculation;

[0108] The variable of the history function is introduced. The history function always stores the maximum strain energy density in the deformation history of the material, ensuring that the strain energy density of any point x in space will not decrease with time t. It can be expressed as:

[0109]

[0110] At this point, the calculation of phase field values, displacement field values, and average stress values ​​based on the displacement boundary and the introduction of phase field method and history function method at the microscopic scale is completed;

[0111] The specific implementation method of step S5 is as follows: the overall modulus of the microscopic model will decrease, and the modulus of the microscopic model at this time needs to be calculated. At this time, the microscopic composite unit cell model is regarded as an elastic body. After the damage is completed, mandatory boundary conditions are applied to the unit cell model to calculate the new modulus of the microscopic model after the damage. The imposed boundary conditions are:

[0112]

[0113] Among them, u + / - (x + / - ) represents the displacement vector applied to the outer surface node perpendicular to the X-axis of the unit cell model, and the node needs to be symmetrical about the YOZ plane, that is, each node to which the displacement is applied needs to find a corresponding node on the symmetry plane. Indicates the mandatory boundaries in 6 different directions, namely:

[0114]

[0115] in To apply a load in a single direction, To apply loads in two orthogonal directions. For the composite unit cell model, there is the following constitutive relationship:

[0116]

[0117] After applying the above six boundary conditions, the constants in the elastic matrix can be obtained one by one. For example,

[0118]

[0119] At this time, only the displacement boundary in the positive direction of the X axis is applied, then,

[0120]

[0121] The specific implementation method of the above step S6 is: after the average stress calculation and modulus of the micro model are updated, the internal force of the unit is calculated based on the average stress calculated at each Gauss point on each unit:

[0122] f int =∫ Ω B T SdΩ

[0123] Among them, S is the average stress inside the unit, and the internal force of each unit is calculated by numerical integration. After completing the calculation of the internal force of each unit, the overall internal force F of the macro model can be obtained according to the index of the node on the unit. in .

[0124] The overall external force on the macro model is:

[0125]

[0126] The specific implementation method of the above step S7 is: according to the incremental form of the equilibrium equation:

[0127] K tan ΔU=F ext -F in

[0128] where K tan is the tangent stiffness matrix for each iteration, which is the element stiffness matrix obtained by combining the modulus calculated from the micromodel at the macro Gaussian point with the element geometry matrix of the macromodel:

[0129]

[0130] Among them, D lower Represents the average modulus of the microscopic model calculated at each Gaussian point on the macroscopic scale.

[0131] The specific implementation method of the above step S9 is: based on the large deformation theory of hyperelastic materials, derive the weak form equilibrium equation under the complete Lagrangian format:

[0132]

[0133] Among them, the internal force virtual work δW int and external force virtual work δW ext are expressed by the second Piola-Kirchhoff stress S and external load, respectively.

[0134] Linearizing the weak form equations, we get the incremental form of the governing equations:

[0135]

[0136] The tangent stiffness matrix K is composed of the linear stiffness matrix K L (material response) and the geometric stiffness matrix K NL (Stress stiffness) composition:

[0137] K=K L +K NL

[0138] The incremental iteration starts from the initial iteration, that is, from the known equilibrium state, and calculates the initial residual:

[0139]

[0140] During the iteration process, the displacement variables are gradually corrected using the Newton-Raphson method:

[0141]

[0142] Update the strain, stress, and internal force vectors after each iteration:

[0143]

[0144] After iteration, the convergence is judged, and the residual norm or displacement increment norm is used as the convergence standard:

[0145]

[0146] As outlined in the text, the residuals are determined by comparing the displacements of the new iteration step with those of the previous iteration step. At the initial iteration step, the residuals are compared with the displacement results of the macro model. If convergence is achieved, the calculation is completed at this displacement load step; if not, the next iteration is performed. The new displacement results replace the displacement results of the second step, and the macro Gaussian point strain is recalculated. The history function values ​​and displacement field values ​​of each micro model are retained and used in the phase field calculation in the next iteration step to achieve phase field evolution.

[0147] The specific implementation method of the above step S9 is as follows: During post-processing, since the variables only involve the elasticity matrix, stress matrix, and strain matrix, and the calculation is only the multiplication of the design matrices, using the GPU's CUDA can more effectively utilize computing resources. The process is as follows:

[0148] Copying CPU data to the CPU: During post-processing, copying the matrix of each unit to the GPU multiple times based on the number of units will reduce computational efficiency. Therefore, first assembling the stress and strain vectors of each unit into a complete matrix in the CPU memory according to the unit order and then copying it to the GPU memory at once can greatly reduce the efficiency loss caused by copying data back and forth.

[0149] GPU kernel configuration: After completing the CPU-to-GPU data transfer, the GPU must also configure the number of thread blocks (Blocks) and the number of threads per block (Threads). To prevent some threads from performing multiple tasks and impacting overall parallelism, the number of threads is generally set to be greater than or equal to the number of elements to be executed. After the GPU kernel is launched, the runtime will launch the appropriate number of threads in a thread grid based on the launch parameters, all executing in parallel.

[0150] The GPU data is copied to the CPU. After the aforementioned matrix calculations are performed on the GPU, the data is copied to the CPU. Since the data is still stored in memory as a one-dimensional vector, the first step is reversed to convert this one-dimensional data into a two-dimensional matrix. After all calculations are completed, the displacement results of the macro model are obtained. The unit deformation matrix and material matrix stored when calculating the unit stiffness matrix are used to parallelize the GPU to calculate the stress and strain contours for this displacement result.

[0151] The beneficial effects of the present invention are as follows: the present invention realizes an efficient algorithm for composite materials based on the advantages of different parallel methods. First, a cross-scale calculation method is proposed based on the finite element theory, connecting the macroscale and the microscale through internal and external forces. In addition, combined with phase field damage and composite modulus prediction, the cross-scale calculation of the mean macro model and anisotropic micro composite materials is realized with the help of homogenization theory, and a heterogeneous parallel acceleration strategy is provided through the parallel method of CPU+GPU+cluster. The CPU's OpenMP is used in parallel for the calculation of the unit stiffness matrix, the cluster's MPI is used in parallel for the allocation of macro units in cross-scale calculations, and the GPU is used in parallel for the stress and strain post-processing after the finite element calculation. Calculations are performed in two-dimensional conditions for different composite micro models. The relationship between different initial pores, cracks and crack extension is discussed, the displacement-load curves under each model are compared, and the fracture work is calculated through the displacement-load curve, which also provides a solution for three-dimensional calculations. Using modern large-scale computing, a large-scale heterogeneous parallel simulation calculation method based on regional decomposition and combined with multi-process, multi-threading, and GPU parallelism is proposed. The "divide and conquer" idea is used to solve the problems of insufficient computing resources and long computing time, thereby realizing the calculation of large-scale composite material damage prediction. BRIEF DESCRIPTION OF THE DRAWINGS

[0152] Figure 1 This is a flow chart of the cross-scale damage heterogeneity parallel method for continuous fiber reinforced ceramic composites of the present invention;

[0153] Figure 2 Schematic diagram of macro and micro structures, where the left side is the macro model and the right side is the micro model;

[0154] Figure 3 The left side is the two-dimensional displacement result cloud map, and the right side is the three-dimensional displacement result cloud map;

[0155] Figure 4 The figure is a parallel diagram of the element stiffness matrix using OpenMP;

[0156] Figure 5 It is a schematic diagram of the parallelization of MPI to macro units;

[0157] Figure 6 Schematic diagram of modulus prediction;

[0158] Figure 7 This is a parallel diagram of GPU post-processing;

[0159] Figure 8 It is a parallel overall framework diagram;

[0160] Figure 9Schematic diagram of damage evolution of two-dimensional and three-dimensional models, where the left side shows the two-dimensional damage evolution process and the right side shows the three-dimensional damage evolution process in the YOZ plane;

[0161] Figure 10 These are the force and displacement result diagrams of the two-dimensional and three-dimensional models. The left side is the damage force and displacement change curve of the two-dimensional model, and the right side is the damage force and displacement change curve of the three-dimensional model. DETAILED DESCRIPTION

[0162] The technical solution of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the embodiments described are only part of the embodiments of the invention, not all of them. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making any creative efforts are within the scope of protection of the present invention.

[0163] like Figure 1 As shown in the figure, a parallel method for cross-scale damage heterogeneity of continuous fiber reinforced ceramic composites is provided. The specific implementation method of step S1 is as follows:

[0164] S11. Establish a finite element model at a macro scale;

[0165] A macroscopic model is constructed, treating the material as a homogeneous material. The material parameters are derived by using a homogenization method to predict the modulus of the undamaged microscopic model. Based on the results, the Poisson's ratio and elastic modulus are set. Boundary conditions and constraints are also set, determining the locations of displacement boundaries and fixed constraints.

[0166] S12. Establish a microscopic model of continuous fiber reinforced ceramic matrix composites with arbitrary weaving patterns;

[0167] Parametric modeling is established using the level set method and non-uniform rational B-spline curve to establish the composite material element for hot-discussing braided driving.

[0168] The B-spline curve describes the fiber position of any weaving method, and the shortest distance between the node and the B-spline curve is calculated using the grid method and gradient descent method. It generates a curve by generating multiple rational polynomials through control points, and its C(t) is analytically defined as:

[0169]

[0170] Through the control points, a B-spline curve with arbitrary spatial distribution can be generated, and the curve can be used as the axis of the fiber. i :i=0,1,...,n} as control points, N i,j (t) is the basis function of the j-th order i B-spline curve described by m non-decreasing knot sequences t, which are t i The j-order piecewise polynomial of control is expressed as:

[0171]

[0172] The B-spline curve with any form of spatial distribution can be generated by the control points. The curve is used as the axis of the fiber. In order to determine whether the unit passes through the interface layer between the fiber and the matrix, it is necessary to obtain the P at any point in space. a The shortest distance d from (a, b, c) to the fiber curve, d is about point P on the fiber curve b The function of (x,y,z) can be expressed as:

[0173]

[0174] The above maximum value problem can be equivalently solved to find the minimum value of its square, which can be written as an equivalent distance function, which is:

[0175]

[0176] The grid method combined with the gradient descent method is used to solve the point P on the fiber axis. b (x, y, z) is a function of parameter t, so the equivalent distance function obtains the gradient of the fiber curve with respect to parameter t

[0177]

[0178] Find the approximate parameter t by grid point method rough , then use t rough Initialize the parameter t in the gradient descent method, use the gradient descent method to continuously approximate the value of t, and continuously modify the parameter t in the iteration until the point P on the fiber curve corresponding to t is bmin (x,y,z) and P a Shortest distance:

[0179]

[0180] The level set method is used to determine the distribution of fibers and the matrix, thereby dividing the fibers and the matrix. For a cubic model discretized using tetrahedral elements, there are one or more fibers described. The position of the fiber surface can be described by the fiber level set function. In the three-dimensional case, the level set method represents the fiber surface curve ζ as the zero level set of the three-dimensional auxiliary function φ:

[0181] ζ={(x,y,z)|φ(x,y,z)=0}

[0182] In order to represent the shape area Ω of the fiber, a level set function φ is defined in three-dimensional space. The zero level set of φ represents the boundary Γ where the fiber and the matrix meet. When φ is less than zero, it represents the inside of the fiber area, and when φ is greater than zero, it represents the outside of the fiber area. To describe the structure of the fiber, it is also necessary to define the cross-sectional shape of the fiber. To define the cross-sectional shape of the fiber, it is necessary to first determine the position of the cross-sectional shape. The first derivative of t is taken from the expression of the fiber axis to obtain the tangent vector expression of the fiber axis.

[0183]

[0184] Take a point o(t) on the fiber axis and use the tangent vector as the x-axis at that point. Use the global z-axis as the z-axis at that point. Use the cross product to calculate the y-axis of the ellipse's cross section. Then, normalize the y-axis vector. Finally, multiply the x-axis by the y-axis to obtain the normalized z-axis vector at that point. Use this equation to calculate the level set function from that point to the fiber interface, and then determine the positional relationship between any point on the plane and the ellipse or rectangle.

[0185] At this point, the finite element model of the macro-scale component is completed, and the B-spline curve method is used to establish the model of the continuous fiber reinforced ceramic matrix composite material with arbitrary weaving method at the micro scale, such as Figure 2 As shown;

[0186] S2, introduce CPU parallel computing and use finite element method to calculate the displacement and strain of all units in the macro model;

[0187] After dividing a continuous domain into multiple discrete domains, that is, dividing a complete model into elements and nodes in the finite element method, for any discrete domain Ω in space, the following three basic equations must be satisfied: equilibrium equations, geometric equations, and physical equations.

[0188] σ ij,j +F i =0(i=1,2,3;j=1,2,3)

[0189] Among them, σ ij is a second-order symmetric stress tensor. The equilibrium equation establishes the physical relationship between the principal stress and shear stress in the three directions of the elastic body in space:

[0190]

[0191] Where u, v, and w represent the displacements in the x, y, and z directions, respectively, in the Cartesian coordinate system. By taking the partial derivatives of the displacements in the corresponding directions, the geometric equations establish the physical relationship between the strain components and the displacement components.

[0192] σ ij =D ijkl εkl (i,j,k,l = 1,2,3)

[0193] where σ ij and ε kl are two-dimensional stress and strain tensors respectively. D ijkl is the elastic four-dimensional tensor in domain Ω. By symmetric four-dimensional elastic tensor, the physical equation establishes the physical relationship between strain components and stress components.

[0194] So far, three basic equations are introduced, which have six directions of stress, six directions of strain and three directions of displacement for elastic body in any three-dimensional space. The calculated displacement results are shown in Figure 3 , which totals 15 unknowns. In order to solve the above equations, force boundary conditions and displacement boundary conditions are also introduced to specify the values and normal derivatives on the boundary.

[0195]

[0196] where n is the direction cosine of the normal of the boundary of domain Ω, is the area force tensor on the boundary. Dirichlet boundary condition is to directly specify the function value on the boundary. So far, the basic equations of elasticity are introduced. Based on the equilibrium equation, the strong form of the differential equation can be derived and written in the form of integral:

[0197]

[0198] For finite element calculation of real physical model domain, the domain needs to be discretized into physical model domain with finite degrees of freedom. And for each discretized physical domain, the relationship between local displacement and global displacement needs to be established, and the physical equation and geometric equation need to be satisfied:

[0199] u = N U e

[0200] ε = B U e

[0201] σ = D ε

[0202] e represents each discretized physical domain, U e is the displacement of each node on the single element. B is the element geometry matrix, and N is the element shape function matrix, which converts the above integral into summation:

[0203]

[0204] Using the variational principle to take the minimum value of potential energy variation simplifies the above form to get the element stiffness matrix K e :

[0205]

[0206] K e =∫ Ω B T DBdΩ e

[0207]

[0208] The stiffness matrix of each unit is only related to its own node coordinates and the material properties of the unit, and there is no data association between the units. Since the serial program of the unit stiffness matrix is ​​parallelized in a loop, OpenMP can be used to parallelize the serial program to calculate the stiffness matrix of each unit, such as Figure 4 shown.

[0209] Each time the stiffness matrix of the microscopic model is assembled, due to the occurrence of damage, the stiffness matrix at this time needs to be updated using the following formula:

[0210] K=∫ Ω (1-d) 2 B T DBdΩ

[0211] Among them, the damage parameter d is the average value of the damage value of each node on the unit.

[0212] S3. Introducing cluster parallel computing, the displacement boundary of the micro unit at the Gauss point is calculated through the strain of the macro unit Gauss point;

[0213] S31. Distribute macro unit processes through MPI:

[0214] After completing the calculation of the macro model, since the calculation of the micro model is only related to the strain of the macro model unit where the micro model is located, the micro models between different macro units will not generate data exchange. Therefore, the MPI parallel technology is used to split the macro unit according to the number of processes, such as Figure 5 As shown, each process calculates a portion of the units of the macro-complete model. This is done to make the macro-units allocated to each process as equal as possible, to prevent some processes from waiting for other processes after completing their calculations.

[0215] S32. Use Gaussian point strain to calculate the displacement boundary of the micro-element model:

[0216] For any element, the actual calculation is the stress and strain at the Gauss point of the element. The value of the complete integral is calculated by weighting each Gauss point, and then the stress and strain values ​​of the rest of the element are obtained by interpolation. That is, for any Gauss point on the element, the strain tensor at the Gauss point is:

[0217]

[0218] The micro-scale model takes a point of the macro-scale model, and uses the strain at the point to calculate the complete strain of the micro-scale model. The micro-scale model strain is the product of the coordinates of all boundary nodes of the micro-scale model and the strain. For any node on the boundary, the following equation is used:

[0219]

[0220] The reference CPU parallelly calculates and distributes each macro-scale model unit, and calculates the micro-scale unit displacement boundary at the Gaussian point of the macro-scale unit through the strain at the Gaussian point;

[0221] S4. According to the displacement boundary, the phase field value, displacement field value and average stress value after damage are calculated through the phase field damage method at the micro-scale;

[0222] S41. Phase field calculation

[0223] The phase field damage in the micro-scale model of the composite material is realized through the phase field method, and then fed back to the material properties of the macro-scale. The core idea is to introduce a continuous phase field variable d(x), which ranges from 0 to 1, representing the state of the material from no damage to complete damage, where d=1 represents complete damage of the material, and d=0 represents no damage of the material.

[0224] The damage field in one dimension can be described by the following formula:

[0225]

[0226] Where l∈R is a scale parameter that controls the crack, not the actual length of the crack. When using the finite element method to calculate the phase field, the scale parameter l in the crack density function is usually related to the grid size h. After derivation and simplification, the expression is:

[0227]

[0228] d(x)-l 2 d″(x)=0

[0229] Since the differential equation is subject to Dirichlet boundary conditions and can be written as the Euler control equation based on the variational principle as follows:

[0230]

[0231] In the formula, W Γ is the Dirichlet boundary condition, i.e. W Γ ={d|d(0)=1;d(x≠0)=0}; Arg represents the value of the independent variable in the formula when the formula in the bracket reaches the minimum value.

[0232] For a two-dimensional finite domain containing a crack, the expression for the phase field value d(x) is:

[0233]

[0234] in, represents the gradient operator, and vector n is the boundary normal. Since the two-dimensional case still needs to satisfy the value that minimizes the crack surface density function, the phase field value at this time also satisfies the Euler control equation. At this time, the Dirichlet boundary can be expressed as W Γ ={d|d(x)=1on x∈Γ}. In addition, the crack density function per unit volume can be expressed as:

[0235]

[0236] The functional of the internal energy of the damaged body can be expressed as:

[0237] Π(u,Γ)=Π d (u,Γ)+Π s (Γ)

[0238] Among them, Π(u,Γ) is the total energy inside the object, Π d (u,Γ) is the elastic potential energy of the cracked elastic body, Π s (Γ) is the fracture energy required to generate a new crack. And the fracture elastic energy is the integral of the unit crack density function in the domain, that is,

[0239]

[0240] Among them, G c is the critical energy release rate. d (u,Γ) is usually expressed as the integral of the strain energy density function inside the elastic body within the domain, and its expression is:

[0241] Π d (u,Γ)=∫ Ω ψ(ε(u),d(x))dΩ

[0242] Where ψ(ε(u),d(x)) is the elastic energy density function, ε(u) is the strain tensor, and u is the displacement vector. When the elastic body is intact, that is, without damage, the elastic strain energy density is expressed as:

[0243]

[0244] Where C represents the fourth-order elasticity tensor.

[0245] Considering that crack generation and expansion mainly come from tensile loads rather than compressive loads, energy dissipation only occurs under tension. Based on this, the free energy density function is decomposed into tension and compression, and only the energy dissipation under the tensile model is considered. Its expression is:

[0246] ψ(ε)=[g(d)+k]ψ + (ε)+ψ - (ε)

[0247] Among them, ψ + / - (ε) represents the strain energy density function in tension and compression mode, respectively. g(d) is the fracture toughness function related to the phase field value, and the phase field comparison is performed.

[0248] S42. Displacement field calculation and mean stress calculation

[0249] After the displacement field calculation is obtained based on the previous displacement boundary calculation of the micro-model, the new displacement of the micro-model at this time will be obtained after the damage calculation is completed. This displacement is used to calculate the average stress of the micro-model at this time:

[0250]

[0251] After the mandatory boundary conditions are applied, the finite element method and homogenization theory are used to calculate the mean stress under each mandatory periodic boundary condition:

[0252]

[0253] S43, process function calculation

[0254] The variable of the history function is introduced. The history function always stores the maximum strain energy density in the deformation history of the material, ensuring that the strain energy density of any point x in space will not decrease with time t. It can be expressed as:

[0255]

[0256] At this point, the calculation of phase field values, displacement field values, and average stress values ​​based on the displacement boundary and the introduction of phase field method and history function method at the microscopic scale is completed;

[0257] S5. Calculate the average modulus of each unit after damage at the micro scale by homogenization method;

[0258] The specific implementation steps are as follows: the overall modulus of the microscopic model will decrease, and the modulus of the microscopic model at this time needs to be calculated. At this time, the microscopic composite unit cell model is regarded as an elastic body. After the damage is completed, mandatory boundary conditions are applied to the unit cell model to calculate the new modulus of the microscopic model after damage. The imposed boundary conditions are:

[0259]

[0260] where u + / - (x + / - ) represents the displacement vector applied on the outer surface node of the unit cell model perpendicular to the X-axis, and the node needs to be symmetric about the YOZ plane, that is, each node with an applied displacement needs to find a corresponding node on the symmetry plane. represents the forced boundary in 6 different directions, respectively:

[0261]

[0262] where is the load applied in a single direction, is the load applied in two orthogonal directions. For the composite unit cell model, the following constitutive relationship exists:

[0263]

[0264] After applying the above six boundary conditions, each constant in the elastic matrix can be obtained one by one. Taking the boundary condition as an example, we have:

[0265]

[0266] And at this time, only the displacement boundary in the positive direction of the X-axis is applied, then we have:

[0267]

[0268] The predicted modulus results can be obtained, and the results are compared as shown in Figure 3 .

[0269] S6, return to the macro scale, use the micro scale unit average stress to calculate the internal force in the macro unit at this time;

[0270] The specific implementation steps are: after the average stress calculation and modulus update of the micro model, return to the macro model to calculate the internal force of the unit according to the average stress calculated by each Gauss point on each unit:

[0271] f int =∫ Ω B T SdΩ

[0272] where S is the average stress inside the unit, and the internal force of each unit is calculated in the same way using numerical integration. After completing the internal force calculation of each unit, the overall internal force F in of the macro model can be obtained according to the index of the node on the unit.

[0273] And the overall external force on the macro model is:

[0274]

[0275] S7, calculating the post-damage macro-element stiffness matrix using the average post-damage modulus of each macro-element;

[0276] The specific implementation steps are: According to the incremental form of the equilibrium equation:

[0277] K tan ΔU=F ext -F in

[0278] where K tan is the tangent stiffness matrix for each iteration, which is the element stiffness matrix obtained by combining the modulus calculated from the micromodel at the macro Gaussian point with the element geometry matrix of the macromodel:

[0279]

[0280] Among them, D lower Represents the average modulus of the microscopic model calculated at each Gaussian point on the macroscopic scale.

[0281] S8, iteratively predict damage evolution using the post-damage global stiffness matrix, internal forces, and external forces calculated from applied loads;

[0282] The specific implementation steps are as follows: Based on the large deformation theory of hyperelastic materials, the weak form equilibrium equation under the complete Lagrangian format is derived:

[0283]

[0284] Among them, the internal force virtual work δW int and external force virtual work δW ext are expressed by the second Piola-Kirchhoff stress S and external load, respectively.

[0285] Linearizing the weak form equations, we get the incremental form of the governing equations:

[0286]

[0287] The tangent stiffness matrix K is composed of the linear stiffness matrix K L (material response) and the geometric stiffness matrix K NL (Stress stiffness) composition:

[0288] K=K L +K NL

[0289] The incremental iteration starts from the initial iteration, that is, from the known equilibrium state, and calculates the initial residual:

[0290]

[0291] During the iteration process, the displacement variables are gradually corrected using the Newton-Raphson method:

[0292]

[0293] Update the strain, stress, and internal force vectors after each iteration:

[0294]

[0295] After iteration, the convergence is judged, and the residual norm or displacement increment norm is used as the convergence standard:

[0296]

[0297] In summary, the residuals are determined by comparing the displacements of the new iteration step with those of the previous iteration step. At the initial iteration step, the residuals are compared with the displacement results of the macro model. If convergence is achieved, the calculation is completed at this displacement load step; if not, the next iteration is performed. The new displacement results replace the displacement results of the second step, and the macro Gaussian point strain is recalculated. The history function values ​​and displacement field values ​​of each micro model are retained and used in the phase field calculation in the next iteration step to achieve phase field evolution.

[0298] S9. GPU parallel computing is introduced to calculate the stress and strain of the macro model according to the displacement value to achieve post-processing.

[0299] The specific implementation method is: in post-processing, since the variables only involve the elasticity matrix, stress matrix and strain matrix, and the calculation is only the multiplication of the design matrix. Using GPU CUDA can more effectively utilize computing resources. The specific process is as follows Figure 6 As shown, the text flow is as follows:

[0300] Copying CPU data to the CPU: During post-processing, copying the matrix of each unit to the GPU multiple times based on the number of units will reduce computational efficiency. Therefore, first assembling the stress and strain vectors of each unit into a complete matrix in the CPU memory according to the unit order and then copying it to the GPU memory at once can greatly reduce the efficiency loss caused by copying data back and forth.

[0301] GPU kernel configuration: After completing the CPU-to-GPU data transfer, the GPU must also configure the number of thread blocks (Blocks) and the number of threads per block (Threads). To prevent some threads from performing multiple tasks and impacting overall parallelism, the number of threads is generally set to be greater than or equal to the number of elements to be executed. After the GPU kernel is launched, the runtime will launch the appropriate number of threads in a thread grid based on the launch parameters, all executing in parallel.

[0302] The GPU data is copied to the CPU. After the calculation of the above matrix is ​​realized in the GPU, the data is copied to the CPU. Since the data is still stored in the memory in the form of a one-dimensional vector at this time, the one-dimensional data is converted into a two-dimensional matrix using the opposite method of the first step. After completing all the calculations, the displacement result of the macro model at this time is obtained. The unit deformation matrix and material matrix stored when calculating the unit stiffness matrix are used to use the GPU to parallelly calculate the stress and strain cloud map under this displacement result, and all the results are exported. The complete parallel calculation block diagram is as follows Figure 8 shown.

[0303] The present invention applies the phase field method and heterogeneous parallel computing strategy to the finite element calculation of composite materials. Based on the advantages of different parallel methods, the numerical computing capabilities of the computer are fully utilized to realize the calculation of large-scale models, and the rapid calculation of damage prediction of large-scale composite materials can be realized. For example, the damage evolution results of two-dimensional and three-dimensional model cases are shown in Figure 2. Figure 9 The force and displacement curves are shown in Figure 10 shown.

[0304] The above description does not limit the present invention in any form. Although the present invention has been disclosed through the above embodiments, it is not intended to limit the present invention. Any person skilled in the art can use the above disclosed technical content to make changes or modifications to equivalent embodiments without departing from the scope of the technical solution of the present invention. However, any simple modifications, equivalent changes, and modifications made to the above embodiments based on the technical essence of the present invention without departing from the content of the technical solution of the present invention are still within the scope of the technical solution of the present invention.

Claims

1. A parallel method for cross-scale damage heterogeneity of continuous fiber reinforced ceramic composites, characterized in that: The following steps are involved: S1. Establish a finite element model of a macro-scale component and use the B-spline curve method to establish a micro-scale model of a continuous fiber reinforced ceramic matrix composite material with arbitrary weaving pattern; S2, introduce CPU parallel computing and use finite element method to calculate the displacement and strain of all units in the macro model; S3. Introducing cluster parallel computing, the displacement boundary of the micro unit at the Gauss point is calculated through the strain of the macro unit Gauss point; S4. Calculate the phase field value, displacement field value and average stress value after damage at the microscopic scale based on the displacement boundary by using the phase field damage method; S5. Calculate the average modulus of each unit after damage at the micro scale by homogenization method; S6. Return to the macro scale and use the average stress of the micro scale unit to calculate the internal force of the macro unit at this time; S7, calculating the post-damage macro-element stiffness matrix using the average post-damage modulus of each macro-element; S8, iteratively predict damage evolution using the post-damage global stiffness matrix, internal forces, and external forces calculated from applied loads; S9. GPU parallel computing is introduced to calculate the stress and strain of the macro model according to the displacement value to achieve post-processing.

2. The method for cross-scale damage heterogeneity parallel treatment of continuous fiber reinforced ceramic composites according to claim 1, characterized in that: The specific implementation method of step S1 is: S11. Establish a finite element model at a macro scale; A macroscopic model that is treated as a homogeneous material is established. The material parameters are obtained by using the homogenization method to predict the modulus of the undamaged microscopic model. The Poisson's ratio and elastic modulus are set according to the results. Boundary conditions and constraints need to be set to determine the locations of the displacement boundaries and fixed constraints. S12. Establish a microscopic model of continuous fiber reinforced ceramic matrix composites with arbitrary weaving patterns; The parametric modeling is established by using the level set method and non-uniform rational B-spline curve to establish the composite material microelement of the hot discussion weaving travel. The B-spline curve describes the fiber position of any weaving method. The grid method and gradient descent method are used to calculate the shortest distance between the node and the B-spline curve. The curve is constructed by generating multiple rational polynomials through control points. Its C(t) is analytically defined as: The B-spline curve with any form of spatial distribution can be generated by controlling the points, and the curve is used as the axis of the fiber. i :i=0,1,...,n} as control points, N i,j (t) is the basis function of the j-th order i B-spline curve described by m non-decreasing knot sequences t, which are t i The j-order piecewise polynomial of control is expressed as: The B-spline curve with any form of spatial distribution can be generated by the control points. The curve is used as the axis of the fiber. In order to determine whether the unit passes through the interface layer between the fiber and the matrix, it is necessary to obtain the P at any point in space. a The shortest distance d from (a, b, c) to the fiber curve, d is about point P on the fiber curve b The function of (x,y,z) can be expressed as: The above maximum value problem can be equivalently solved to find the minimum value of its square, which can be written as an equivalent distance function, which is: The grid method combined with the gradient descent method is used to solve the point P on the fiber axis. b (x, y, z) is a function of parameter t, so the equivalent distance function obtains the gradient of the fiber curve with respect to parameter t Find the approximate parameter t by grid point method rough , then use t rough Initialize the parameter t in the gradient descent method, use the gradient descent method to continuously approximate the value of t, and continuously modify the parameter t in the iteration until the point P on the fiber curve corresponding to t is bmin (x,y,z) and P a Shortest distance: The level set method is used to determine the distribution of fibers and matrix, thereby dividing the fibers and matrix. For a cubic model discretized using tetrahedral elements, there are one or more fiber descriptions. The position of the fiber surface can be described by the fiber level set function. In the three-dimensional case, the level set method represents the fiber surface surface ζ as the zero level set of the three-dimensional auxiliary function φ: ζ={(x,y,z)|φ(x,y,z)=0} In order to represent the shape area Ω of the fiber, a level set function φ is defined in three-dimensional space. The zero level set of φ represents the boundary Γ between the fiber and the matrix. When φ is less than zero, it means it is inside the fiber area, and when φ is greater than zero, it means it is outside the fiber area. To describe the structure of the fiber, it is also necessary to define the cross-sectional shape of the fiber. To define the cross-sectional shape of the fiber, it is necessary to first determine the position of the cross-sectional shape. The first derivative of t is obtained from the expression of the fiber axis to obtain the tangent vector expression of the fiber axis. Take a point o(t) on the fiber axis and use the tangent vector as the x'axis at the point, the global z-axis direction as the z'axis direction at the point, use the cross product to calculate the y'axis of the elliptical section, and then normalize the y-axis vector. Finally, multiply the y-axis by the x-axis cross product to get the z'axis normalized vector at the point. Calculate the level set function value from the point to the fiber interface through the equation, and then obtain the positional relationship between any point on the plane and the ellipse or rectangle. At this point, the finite element model of the macro-scale component is completed, and the B-spline curve method is used to establish a model of continuous fiber reinforced ceramic matrix composite materials with arbitrary weaving methods at the micro scale.

3. The cross-scale damage heterogeneity parallel method for continuous fiber reinforced ceramic composites according to claim 1, characterized in that: The specific implementation method of step S2 is: after dividing a continuous domain into multiple discrete domains, that is, dividing a complete model into units and nodes in the finite element method, for any discrete domain Ω in space, the following three basic equations need to be satisfied, namely, the equilibrium equation, the geometric equation and the physical equation: s ij,j +F i =0(i=1,2,3;j=1,2,3) Among them, σ ij It is a second-order symmetric stress tensor. The equilibrium equation establishes the physical relationship between the principal stress and shear stress in the three directions of the elastic body in space. Among them, u, v, and w represent the displacement in the x, y, and z directions in the Cartesian coordinate system respectively. By taking the partial derivatives of the displacement in the corresponding directions, the geometric equation establishes the physical relationship between the strain component and the displacement component. s ij =D ijkl e kl (i,j,k,l=1,2,3) Among them, σ ij and ε kl are the two-dimensional stress and strain tensors, D ijkl It is the elastic four-dimensional tensor in the domain Ω. Through the symmetric four-dimensional elastic tensor, the physical equation establishes the physical relationship between the strain component and the stress component. So far, three basic equations have been introduced. In any three-dimensional space, there are 15 unknowns for the elastic body with stress in six directions, strain in six directions and displacement in three directions. In order to solve the above equations, it is necessary to introduce force boundary conditions and displacement boundary conditions to specify the values ​​and normal derivatives on the boundaries. where n is the direction cosine of the boundary normal of the domain Ω, is the area force tensor on the boundary, and the Dirichlet boundary condition directly specifies the function value on the boundary. So far, the basic equations of elasticity have been introduced. Based on the equilibrium equation, the strong form of the differential equation can be derived and written in the form of an integral: Finite element calculations for real physical model domains require discretizing the domain into physical model domains with finite degrees of freedom. For each discretized physical domain, the relationship between local displacement and global displacement must be established, and the physical and geometric equations must be satisfied: you=NOW e ε=THIS e σ=Dε e represents each physical domain after being discretized, U e is the displacement of each node on a single element, B is the element geometry matrix, N is the element shape function matrix, and the above integral is converted into a summation method: Using the variational principle to take the minimum value of the potential energy variation to simplify the above form, the unit stiffness matrix K is obtained e : K e =∫ Ω B T DBdΩ e Each time the stiffness matrix of the microscopic model is assembled, due to the occurrence of damage, the stiffness matrix at this time needs to be updated using the following formula: K=∫ Ω (1-d) 2 B T DBdΩ Among them, the damage parameter d is the average value of the damage value of each node on the unit. The stiffness matrix of each unit is only related to its own node coordinates and the material properties of the unit, and there is no data association between the units. Since the serial program of the unit stiffness matrix is ​​parallelized in a loop, OpenMP can be used to parallelize the serial program to calculate the stiffness matrix of each unit. At this point, the introduction of CPU parallel computing is completed, and the finite element method is used to calculate the displacement and strain of all units in the macro model.

4. The method for cross-scale damage heterogeneity parallel treatment of continuous fiber reinforced ceramic composites according to claim 1, characterized in that: The specific implementation method of step S3 is: S31, distribute macro unit processes through MPI; After completing the calculation of the macro model, since the calculation of the micro model is only related to the strain of the macro model unit where the micro model is located, the micro models between different macro units will not generate data exchange. Therefore, MPI parallel technology is used to split the macro unit according to the number of processes. Each process calculates a part of the unit of the macro complete model. In order to make the macro units allocated to each process as equal as possible, it is prevented that some processes have to wait for other processes after the calculation is completed. S32. Use Gaussian point strain to calculate the displacement boundary of the micro-element model; For any element, the actual calculation is the stress and strain at the Gaussian point of the element. The value of the complete integral is weighted by each Gaussian point, and then the stress and strain values ​​of the rest of the element are obtained by interpolation. That is, for any Gaussian point on the element, the strain tensor at the Gaussian point is: The micro-scale model is regarded as a point in the macro-model, and the strain at this point is used as the complete strain of the micro-scale model for calculation. The strain of the micro-model is the coordinates of all boundary nodes of the micro-model multiplied by the strain. For any node on the boundary, At this point, the CPU parallel calculation distribution of each macro model unit is completed, and the displacement boundary of the micro unit at the Gauss point is calculated through the Gauss point strain of the macro unit.

5. The method for cross-scale damage heterogeneity parallel treatment of continuous fiber reinforced ceramic composites according to claim 1, characterized in that: The specific implementation method of step S4 is: S41, phase field calculation; The phase field method is used to realize phase field damage in the microscopic model of composite materials, and then feedback to the macroscopic material properties. The core idea is to introduce a continuous phase field variable d(x). The value range of this variable is [0, 1], which represents the state of the material from no damage to complete damage, where d = 1 means the material is completely damaged, and d = 0 means the material is not damaged. The damage field in one dimension can be described by the following formula: Among them, l∈R is the scale parameter that controls the crack rather than the actual crack length. When using the finite element method to calculate the phase field, the scale parameter l in the crack density function is usually related to the grid size h. After derivation and simplification, the expression is: d(x)-l 2 d″(x)=0 Since the differential equation obeys the Dirichlet boundary conditions and can be written as the following Euler governing equation based on the variational principle, Where W Γ is the Dirichlet boundary condition, that is, W Γ ={d|d(0)=1; d(x≠0)=0}; Arg represents the value of the independent variable when the formula in the brackets reaches the minimum value. For a two-dimensional finite domain containing a crack, the expression for the phase field value d(x) is: in, represents the gradient operator, and vector n is the boundary normal. Since the two-dimensional case still needs to satisfy the value that minimizes the crack surface density function, the phase field value at this time also satisfies the Euler control equation. At this time, the Dirichlet boundary can be expressed as W Γ ={d|d(x)=1on x∈Γ}, In addition, the unit volume crack density function can be expressed as, The functional of the internal energy of the damaged body can be expressed as: Π(u,Γ)=Π d (u,C)+P s (C) Among them, Π(u,Γ) is the total energy inside the object, Π d (u,Γ) is the elastic potential energy of the cracked elastic body, Π s (Γ) is the fracture energy required to generate a new crack, and the fracture elastic energy is the integral of the unit crack density function in the domain, that is, Among them, G c is the critical energy release rate, Π d (u,Γ) is usually expressed as the integral of the strain energy density function inside the elastic body within the domain, and its expression is: P d (u,Γ)=∫ Ω ψ(ε(u),d(x))dΩ Among them, ψ(ε(u),d(x)) is the elastic energy density function, ε(u) is the strain tensor, and u is the displacement vector. When the elastic body is intact, that is, without damage, the elastic strain energy density is expressed as: Where C represents the fourth-order elastic tensor, Considering that crack generation and expansion mainly come from tensile loads rather than compressive loads, energy dissipation only occurs under tension. Based on this, the free energy density function is decomposed into tension and compression, and only the energy dissipation under the tensile model is considered. Its expression is: ψ(ε)=[g(d)+k]ψ + (e)+ψ - (e) Among them, ψ + / - (ε) represents the strain energy density function in tension and compression mode, g(d) is the fracture toughness function related to the phase field value, S42, displacement field calculation and mean stress calculation; After the displacement field calculation is obtained based on the previous displacement boundary calculation of the micro-model, the new displacement of the micro-model at this time will be obtained after the damage calculation is completed. This displacement is used to calculate the average stress of the micro-model at this time: After the mandatory boundary conditions are applied, the finite element method and homogenization theory are used to calculate the mean stress under each mandatory periodic boundary condition: S43, process function calculation; The variable of the history function is introduced. The history function always stores the maximum strain energy density in the deformation history of the material, ensuring that the strain energy density of any point x in space will not decrease with time t. It can be expressed as: At this point, the calculation of phase field values, displacement field values ​​and average stress values ​​based on the displacement boundary and the introduction of phase field method and history function method at the microscopic scale is completed.

6. The method for cross-scale damage heterogeneity parallel treatment of continuous fiber reinforced ceramic composites according to claim 1, characterized in that: The specific implementation method of step S5 is as follows: the overall modulus of the microscopic model will decrease, and the modulus of the microscopic model at this time needs to be calculated. At this time, the microscopic composite unit cell model is regarded as an elastic body. After the damage is completed, a mandatory boundary condition is applied to the unit cell model to calculate the new modulus of the microscopic model after the damage. The imposed boundary condition is: Among them, u + / - (x + / - ) represents the displacement vector applied to the outer surface node of the unit cell model perpendicular to the X axis, and the node needs to be symmetrical about the YOZ plane, that is, each node to which the displacement is applied needs to find a corresponding node on the symmetry plane. Indicates the mandatory boundaries in 6 different directions, namely: in To apply a load in a single direction, To apply loads in two orthogonal directions, the following constitutive relations are used for the composite unit cell model: After applying the above six boundary conditions, the constants in the elastic matrix can be obtained one by one to apply the boundary conditions For example, At this time, only the displacement boundary in the positive direction of the X axis is applied, then, 7. The method for cross-scale damage heterogeneity parallel treatment of continuous fiber reinforced ceramic composites according to claim 1, characterized in that: The specific implementation method of step S6 is: after the average stress calculation and modulus of the micro model are updated, the internal force of the unit is calculated based on the average stress calculated at each Gauss point on each unit. f int =∫ Ω B T SdΩ Among them, S is the average stress inside the unit. The internal force of each unit is calculated by numerical integration. After the internal force calculation of each unit is completed, the overall internal force F of the macro model can be obtained according to the index of the node on the unit. in , The overall external force on the macro model is:

8. The method for cross-scale damage heterogeneity parallel treatment of continuous fiber reinforced ceramic composites according to claim 1, characterized in that: The specific implementation method of step S7 is: according to the incremental form of the equilibrium equation: K tan ΔU=F ext -F in where K tan is the tangent stiffness matrix for each iteration, which is the element stiffness matrix obtained by combining the modulus calculated from the micromodel at the macro Gaussian point with the element geometry matrix of the macromodel: Among them, D lower Represents the average modulus of the microscopic model calculated at each Gaussian point on the macroscopic scale.

9. The method for cross-scale damage heterogeneity parallel treatment of continuous fiber reinforced ceramic composites according to claim 1, characterized in that: The specific implementation method of step S9 is: based on the large deformation theory of hyperelastic materials, deriving the weak form equilibrium equation under the complete Lagrangian format: Among them, the internal force virtual work δW int and external force virtual work δW ext Expressed by the second Piola-Kirchhoff stress S and external load respectively, Linearizing the weak form equations, we get the incremental form of the governing equations: The tangent stiffness matrix K is composed of the linear stiffness matrix K L (material response) and the geometric stiffness matrix K NL (Stress stiffness) composition: K=K L +K NL The incremental iteration starts from the initial iteration, that is, from the known equilibrium state, and calculates the initial residual: During the iteration process, the displacement variables are gradually corrected using the Newton-Raphson method: Update the strain, stress, and internal force vectors after each iteration: After iteration, the convergence is judged, and the residual norm or displacement increment norm is used as the convergence standard: According to the text summary, the new iterative step displacement and the previous iterative step displacement are used to judge the residual. At the initial iterative step, the residual is judged with the macro model displacement result. If converged, the calculation is completed in this displacement load step. If not converged, enter the next iterative step, and replace the displacement result of the second step with the new displacement result at this time to recalculate the macro Gaussian point strain. The history function value and displacement field value of each micro model at this time are retained, and the above values ​​are used to calculate the phase field value in the next iterative step to realize the evolution of the phase field value.

10. The method for cross-scale damage heterogeneity parallel treatment of continuous fiber reinforced ceramic composites according to claim 1, characterized in that: The specific implementation method of step S9 is as follows: in post-processing, since the variables only involve the elasticity matrix, the stress matrix and the strain matrix, and the calculation is only the multiplication of the design matrices, using the CUDA of the GPU can more effectively utilize the computing resources. The process is as follows: Copying CPU data to CPU: During post-processing calculations, if the matrix of each unit is copied to the GPU multiple times according to the number of units and then calculated, the calculation efficiency will be reduced. Therefore, first assemble the stress and strain vectors of each unit into a complete matrix in the CPU memory according to the unit order, and then copy it to the GPU memory at one time, which can greatly reduce the efficiency reduction caused by copying data back and forth. Setting up the GPU kernel function: After completing the data transfer from the CPU to the GPU, the number of thread blocks (Block) and the number of threads on each thread block (Thread) need to be set in the GPU. In order to prevent some threads from performing multiple tasks and affecting the overall parallel capability, the number of threads is generally greater than or equal to the number of elements to be executed. After the GPU kernel is started, the runtime will start the corresponding number of threads in a thread grid according to the startup parameters, and they are all executed in parallel. The GPU data is copied to the CPU. After the above matrix calculation is implemented in the GPU, the data is copied to the CPU. Since the data is still stored in the memory in the form of a one-dimensional vector at this time, the one-dimensional data is converted into a two-dimensional matrix using the opposite method of the first step. After completing all calculations, the displacement result of the macro model at this time is obtained. The unit deformation matrix and material matrix stored when calculating the unit stiffness matrix are used to use the GPU to parallelly calculate the stress and strain cloud diagrams under this displacement result.

Images