A multi-domain and multi-level collaborative simulation and life prediction method for complex customized products

Abstract

The present invention proposes a multi-domain and multi-level collaborative simulation and life prediction method for complex customized products, and the method comprises the following steps: S1, establishing a finite element model of the customized product; S2, importing the model for system dynamics simulation for subsequent model evaluation; S3, importing the model for modal analysis, and obtaining modal result-related files; S4, importing the modal results into the dynamics part, establishing and analyzing a rigid-flexible coupling model of the customized product, and combining the system dynamics simulation results to obtain a performance evaluation of the customized product model and a time-load history used for subsequent life prediction; S5, importing the model and establishing a fluid area for flow field analysis; S6, obtaining the fluid analysis results and combining the customized product finite element model, calculating the fluid pressure as the boundary condition of the surface of the customized product in fluid-solid coupling; S7, loading the fluid-solid coupling boundary conditions and importing it into LiToSim for fluid-solid coupling analysis; S8, importing the modal analysis file for acoustic simulation.

Description

Technical Field

[0001] The present invention belongs to the research field of multi-domain and multi-level simulation, and specifically is a method for realizing multi-domain and multi-level collaborative simulation and life prediction of customized products. Background Art

[0002] With the integration and application of new generation information technology and manufacturing, the environment, users, technical requirements, etc. are constantly changing, and product demand is becoming more diversified and personalized. Researching and building a set of innovative product R&D technology systems and technical processes to effectively support efficient and high-quality product development while creating a good user experience and supporting the comprehensive upgrade of R&D products is an urgent problem to be solved in the transformation and upgrading of the manufacturing industry.

[0003] Customized products have many uncertain factors, and traditional simulation methods are insufficient in data information interaction and user experience, and the verification cycle is long and the cost is high. A standardized and unified model is used throughout the entire design process, model information is extracted for different functions / performance needs, a physical simulation model is established, and mixed reality simulation technology is used to jointly build a real virtual environment using model-based multi-physics fields (mechanical, electrical, hydraulic, thermal, magnetic, etc.) in CPS to achieve modeling and simulation of customized products in the information space, and iterate with the verification information in the physical space to achieve comprehensive product optimization.

[0004] Establishing a multi-domain and multi-level collaborative simulation and life prediction method for customized products and forming a corresponding simulation process can not only simplify the simulation process of customized products, but also clarify their operating conditions in complex environments and form iterative optimization of customized products.

[0005] Therefore, there is an urgent need for a multi-domain and multi-level collaborative simulation and life prediction method for complex customized products. Summary of the invention

[0006] The purpose of the present invention is to provide a multi-domain and multi-level collaborative simulation and life prediction method for customized products. The method, combined with the customized product structure model, can realize its simulation analysis in multiple different physical domains and combine some of the physical domains to achieve collaborative simulation; through the multi-domain and multi-level collaborative simulation and life prediction results, it is judged whether the design of the customized product is reasonable and optimized.

[0007] The technical solution of the present invention is implemented as follows: a multi-domain and multi-level collaborative simulation and life prediction method for complex customized products, the method comprising the following steps

[0008] S1. Establish a finite element model of customized products;

[0009] S2, import the model for system dynamics simulation for subsequent model evaluation;

[0010] S3. Import the model for modal analysis and obtain the modal result related files;

[0011] S4. Import the modal results into the dynamics part, establish and analyze the rigid-flexible coupling model of the customized product, and combine the system dynamics simulation results to obtain the performance evaluation of the customized product model and the time-load history used for subsequent life prediction;

[0012] S5, import the model and establish the fluid region for flow field analysis;

[0013] S6. Obtain the fluid analysis results and combine them with the finite element model of the customized product to calculate the fluid pressure as the boundary condition of the surface of the customized product in the fluid-solid coupling;

[0014] S7, load the fluid-structure coupling boundary conditions and import them into LiToSim for fluid-structure coupling analysis;

[0015] S8. Import modal analysis files for acoustic simulation;

[0016] S9. Import the model for static analysis to obtain strength and displacement results;

[0017] S10, introduce the time-load history and static strength simulation results, establish the rain flow counting method and various load spectrum compilation methods;

[0018] S11. Based on the results of load spectrum compilation, various damage models are established for damage prediction, so as to carry out customized product life prediction.

[0019] A further technical solution is that the specific implementation method of step S1 is: establishing a customized product geometric model and completing the establishment of a customized product finite element model.

[0020] A further technical solution is that the specific implementation method of step S2 is:

[0021] S21. Establish a system dynamics simulation model;

[0022] S22. Establish system dynamics simulation.

[0023] A further technical solution is that the specific implementation method of step S3 is: importing a finite element model and establishing a modal analysis.

[0024] A further technical solution is that the specific implementation method of step S4 is:

[0025] S41. Import modal analysis result file;

[0026] S42, establish a rigid-flexible coupling simulation model and analyze it;

[0027] S43, conduct customized product performance evaluation;

[0028] S44. Save the customized product time-load history.

[0029] A further technical solution is that the specific implementation method of step S5 is:

[0030] S51, establishing an external flow field region, which is usually a rectangular parallelepiped;

[0031] S52, hollowing out the customized product model in the outer flow field area;

[0032] S53, setting turbulence model;

[0033] The customized product in this paper adopts the turbulent SST k-ω model, and its control equation is as follows:

[0034]

[0035]

[0036] Where: G k represents the turbulent kinetic energy, G ω is the ω equation, Γ k , Γ ω represent the effective diffusion terms of k and ω respectively, Y k , Y ω Represent the divergence terms of k and ω respectively. ω represents the orthogonal divergence term. S k With S ω is a user-defined term. The SST k-ω model is different from the standard k-ω model in that in the standard k-ω, α ∞ is a constant, while in the SST model, α ∞ The equation is as follows:

[0037] α ∞ =F1α ∞,1 +(1-F1)α ∞,2

[0038] Where:

[0039] SST k-ω is similar to the k-ω model, but with the following improvements:

[0040] (1) The deformation of the SST k-ω model and the k-ω model grows from the addition of the hybrid function and the dual model. The hybrid function is designed for the near-wall region, where the standard k-ω model is valid. There are also free surfaces, which are valid for the deformation of the k-ω model;

[0041] (2) The SST k-ω model incorporates the cross-diffusion from the ω equation;

[0042] (3) Turbulent viscosity takes into account the propagation of turbulent shear stress;

[0043] (4) The model constants are different.

[0044] These modifications make the SST k-ω model more accurate and reliable than the standard k-ω model in a wide range of flow fields.

[0045] S54, set boundary conditions U, P, k, nut, omega, etc.;

[0046] S55. Set the algorithm type to SIMPLE.

[0047] S56, solving the equation by using the NS equation;

[0048] The common form of the NS equation is as follows:

[0049]

[0050] Where: is the volume force; is the pressure differential force; For the viscosity force.

[0051] A further technical solution is that the specific implementation method of step S6 is:

[0052] S61, processing the flow field simulation result data file for subsequent calculation;

[0053] S62, calculating fluid node pressure;

[0054] Since the flow field solver uses the finite volume method for calculation, the calculated pressure is a scalar stored at the center of the volume unit. To perform fluid-solid coupling mapping, the pressure value of the fluid node must be calculated first.

[0055] In the flow field, the flow field nodes are mainly divided into two categories: nodes located on the flow field boundary; nodes located inside the flow field. The calculation of the flow field node pressure value is different for the two types of nodes:

[0056] (1) Nodes located on the flow field boundary:

[0057]

[0058] Where n2 is the node located on the boundary of the flow field; bf1 and bf7 are the values ​​of the boundary surfaces of units c1 and c7 respectively.

[0059] (2) Nodes located inside the flow field:

[0060]

[0061] Where: n1 is the node inside the flow field; c1-c7 is the unit value; V c1 -V c7 is the unit volume.

[0062] S63, searching for the fluid node closest to the point to be interpolated;

[0063] All solid surface nodes are points to be interpolated. Before interpolation, we must first find the fluid node closest to the point to be interpolated, and use the octree method to search.

[0064] The octree method mainly performs search iterations at a given depth, reduces the search range through continuous iterations, and eventually searches for the desired result.

[0065] S64, mapping of points to be interpolated;

[0066] According to the nearest fluid node found in S63, all surface elements containing the node are found, and the point to be interpolated is mapped to the plane to which one of the surface elements belongs:

[0067] Define the space plane equation as:

[0068] Ax+By+Cz+D=0

[0069] Where x, y, z are spatial coordinates, and A, B, C, D are equation coefficients.

[0070] Then substitute the three points of the surface unit into the spatial plane equation and solve the homogeneous linear equations to obtain the equation coefficients:

[0071] Ax1+By1+Cz1+D=0

[0072] Ax2+By2+Cz2+D=0

[0073] Ax3+By3+Cz3+D=0

[0074]

[0075] Where: x1, x2, x3 are the x coordinates of the three points, and the same goes for y and z.

[0076] After obtaining A, B, C, and D, the coordinates of the mapping point of the interpolation point in the plane are obtained according to the coordinate substitution formula:

[0077]

[0078]

[0079]

[0080] Then we need to determine whether the mapping point of the interpolation point is within the surface unit:

[0081] Suppose a point P in space x is a point on the plane where the triangular unit P0P1P2 is located. Since the sides P0P1 and P0P2 of ΔP0P1P2 are not collinear, if the vector and As the basis, we know from linear algebra theory that the vector It can be expressed as:

[0082]

[0083] Where: i (i=1,2) is the unknown coefficient.

[0084] Multiply both sides of the above equation by The equations are:

[0085]

[0086] Solving this equation yields λ i (i=1,2), if λ is satisfied i satisfy:

[0087]

[0088] This indicates that point P x is within ΔP0P1P2, otherwise it indicates that the point is outside ΔP0P1P2.

[0089] S65, shape function interpolation of the point to be interpolated

[0090] According to S64, the mapping point coordinates (x p ,y p ,z p ), according to the invariance of spatial distance, the coordinates of triangles and mapping points are converted into a two-dimensional coordinate system.

[0091] Assume that the three points of the triangle unit are P0, P1, P2, and the mapping point is P x Define P0 as the origin P0'(0,0) in the two-dimensional coordinate system, then the two-dimensional coordinates of P1' are (0, P0P1), and the coordinates of P2' can be defined as follows:

[0092] P2'P1'=P2P1

[0093] P2'P0'=P2P0

[0094] The P2' coordinates are obtained by solving the above equations.

[0095] After completing the conversion of three-dimensional coordinates to two-dimensional coordinates, the triangle shape function is defined as:

[0096]

[0097] a i =x j y m -x m y j

[0098] b i =y j -y m

[0099] c i =x m -x j

[0100] Where: A is the area of ​​the triangle; x j 、x m are the x coordinates of two points other than point i, and y j ,y m is the y coordinate.

[0101] According to the shape function, point P x Pressure is defined as:

[0102]

[0103] Where: n is the number of nodes; N i is the shape function of each point; P i For each point pressure.

[0104] S66, Solving Solid Boundary Loading

[0105] The solid node pressure obtained in S65 cannot be directly applied to the solid boundary as a load and requires some processing.

[0106] In this paper, the solid unit is defined as a hexahedral unit, and each of its surface units is a quadrilateral surface unit. Four Gaussian points are defined in the quadrilateral surface unit. The pressure value at the Gaussian point can be obtained by shape function interpolation. The shape function of the quadrilateral is defined as:

[0107]

[0108]

[0109]

[0110]

[0111] Where: N i(i=1,2,3,4) is the shape function of each point of the quadrilateral; s and r are local coordinates.

[0112] Then define the solid boundary node loading as:

[0113]

[0114] Where: w s 、w r is the Gaussian point weight coefficient, which can be obtained by looking up the Gaussian point table; S is the area of ​​the quadrilateral; N i is the shape function of each point; P i is the pressure at each Gaussian point; J is the Jacobian matrix.

[0115] A further technical solution is that the specific implementation method of step S7 is: outputting the loading result calculated in S6 as a command stream file and calling the calculation.

[0116] A further technical solution is that the specific implementation method of step S8 is: importing the modal result file and establishing the boundary element method for acoustic simulation.

[0117] A further technical solution is that the specific implementation method of step S9 is: importing a finite element model, applying boundary conditions and completing static simulation to obtain strength and displacement results.

[0118] A further technical solution is that the specific implementation method of step S10 is:

[0119] S101, introduce the time-load history in S4 and the static simulation results in S9,

[0120] The four-point cycle counting method uses four consecutive vertices to extract stress cycles. First, four consecutive vertices are taken from the stress time history, and the stress values ​​are (S1, S2, S3, S4) at a time; the stress range of the outer cycle is ΔS0 = |S1-S4|.

[0121] (1) If ΔS1≤ΔS0, extract the inner layer cyclic stress S2→S3; and directly connect the outer layer vertices S1 and S4. Go to step (3).

[0122] (2) If ΔS1≤ΔS0, do not extract the cycle; go to step (3);

[0123] (3) Select 4 consecutive vertices from the remaining vertices and recalculate until all loops are extracted.

[0124] Based on the above principle of rainflow counting method, it is used to extract the amplitude and mean of stress cycle.

[0125] After completing the stress cycle extraction, the load spectrum can be compiled. This article uses four load spectrum compilation methods:

[0126] (1) One-dimensional load spectrum compilation

[0127] First, the load-time history is scanned, and the load cycle and its maximum amplitude are obtained by the rain flow counting method. Since the mean value has little effect on fatigue life, the one-dimensional load spectrum does not consider the influence of the mean value.

[0128] Define the distance between groups: D = (F amax -F amin ) / N

[0129] Where D is the group spacing, Famax and Famin are the maximum and minimum values ​​of the load amplitude, and N is the number of load amplitude groups, which is usually defined as 8, 16, or 32.

[0130] After calculating the group spacing, the upper and lower limits of the amplitude are calculated based on the group spacing:

[0131]

[0132] After determining the upper and lower limits of each level, the load amplitude is graded. The frequency of the group in which the amplitude is located is increased by 1.

[0133] (2) Preparation of two-dimensional load spectrum

[0134] The two-dimensional load spectrum is also called the mean-amplitude spectrum. It introduces the influence of the mean value on fatigue life on the basis of the one-dimensional load spectrum, making the calculation results more accurate.

[0135] The two-dimensional load spectrum compilation introduces the grouping of mean values ​​on the basis of one-dimensional load spectrum compilation, and defines the formula for the mean group spacing:

[0136]

[0137] Calculate the upper and lower limits for each group:

[0138]

[0139] After the calculation is completed, determine which group the load mean and amplitude belong to and add 1 to the frequency of that group.

[0140] (3) Fluctuation center method

[0141] The fluctuation center method simplifies the mean and amplitude of the load-time history of the structure into a univariate random variable and integrates all of the two. The specific method is: average all the load cycle means and use it as the fluctuation center, making it exist as the static part, and then superimpose the amplitude as the dynamic part on the static part. The load data obtained in this way is used to compile the load spectrum to simplify the distribution of the mean.

[0142] The wave center method takes the average of all loads and calculates the total average value:

[0143]

[0144] Among them, σ mi is the mean value of each load cycle, V i is the frequency of each mean.

[0145] (4) Variable Mean Method

[0146] The wave center method is generally applicable to symmetrical load-time histories. For asymmetrical load-time histories, the variable mean method is usually used.

[0147] The variable mean method first groups the load cycles according to the amplitude, and then calculates the average value of each group. The average value calculation formula is:

[0148]

[0149] After the load spectrum is recompiled, fatigue life calculation can be performed.

[0150] S102. Establish multiple damage models and complete life prediction;

[0151] (1) Linear damage accumulation

[0152] The linear damage accumulation model is also known as the Miner-Palmgren theory, or Miner damage theory for short. It refers to the continuous accumulation of fatigue damage to a structure when subjected to a cyclic load of a certain amplitude. The damage generated each time is assumed to be unrelated to the stress received, and other factors are not considered, and a linear relationship is directly accumulated. Under a set of loads with varying amplitudes, the damage caused to the structure by n cyclic loads is:

[0153]

[0154] Where n is the number of independent stress cycles; i is the number of times a single stress acts; N i is the fatigue limit life under a single stress.

[0155] (2) Bilinear damage accumulation

[0156] The bilinear damage accumulation formula proposed by Manson-Grover in 1981 is as follows:

[0157] N Ⅰ =Nexp(ZN φ )

[0158] N Ⅱ =NN Ⅰ

[0159] Among them, z and is the model constant, and the calculation formula is as follows:

[0160]

[0161]

[0162] The damage calculation formulas for each stage under multi-level loads are as follows:

[0163]

[0164] (3) Nonlinear damage accumulation

[0165] The exponential loss model proposed by Dolan and Carten in 1956 requires the following assumptions:

[0166] 1) A crack nucleation period is required to produce permanent fatigue damage;

[0167] 2) As the stress gradually increases, the number of crack cores in the specimen will increase;

[0168] 3) The amount of damage increases with the increase of stress;

[0169] 4) With the increase of stress, the damage growth rate of each cycle increases;

[0170] 5) For all load processes, the total amount of damage that causes fatigue failure of the structure is constant;

[0171] 6) When the stress level is less than the limit value of fatigue damage, the damage continues to accumulate.

[0172] Under a given stress s, the accumulated fatigue damage D is expressed by an exponential relationship:

[0173] D=mrn a

[0174] Wherein, the number of fatigue crack nuclei is represented by m, the damage growth coefficient related to the stress level is represented by r, the number of stress cycles is represented by n, and the exponential constant related to the stress level is represented by a.

[0175] When two levels of stress s1 and s2 of different magnitude act on the specimen in an alternating manner, assuming that stress s1 undergoes an cycles, then (1-a)n is used to represent the number of cycles of stress s2. The estimated formula for fatigue life under two levels of stress is:

[0176]

[0177] Extending the above two-level stress fatigue life estimation formula to multi-level loads, we get:

[0178]

[0179] Where N g is the total fatigue life; N1 is the fatigue life under stress s1; each stress level s i The ratio of the number of cycles under the condition to the total number of cycles is expressed as a i Indicates; the maximum stress amplitude under multi-level alternating loads is represented by s1, r is the stress level, d is the material constant, which is usually determined by a two-level load test. Corten and Dolan recommend taking 4.8 as the normal value.

[0180] After adopting the above technical solution, the beneficial effects of the present invention are as follows: the present invention is based on the multi-domain multi-level collaborative simulation and fatigue life prediction of customized products, and uses a multi-domain solver including fluid domain, solid domain, acoustic field, rigid-flexible coupling and fatigue life prediction to complete the multi-level collaborative simulation of customized products; and considers the mapping law of fluid-solid coupling, as well as various compilation methods of fatigue load spectra and various damage models. The entire multi-domain multi-level collaborative simulation and life prediction are completed in the software. BRIEF DESCRIPTION OF THE DRAWINGS

[0181] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative labor.

[0182] Figure 1 It is a flow chart of a multi-domain and multi-level collaborative simulation and life prediction method for complex customized products of the present invention;

[0183] Figure 2 To customize product models and their finite element models;

[0184] Figure 3 Establishment of external flow field model and its grid schematic diagram;

[0185] Figure 4 This is a schematic diagram of flow field simulation results;

[0186] Figure 5 This is a schematic diagram of flow field node pressure calculation;

[0187] Figure 6 This is the idea diagram of the octree method;

[0188] Figure 7 It is a schematic diagram of converting a three-dimensional plane into a two-dimensional plane;

[0189] Figure 8 A schematic diagram for mapping the points to be interpolated to the interior of the plane;

[0190] Fig. 9 This is the idea diagram of the rain flow counting method;

[0191] Fig.10 Prepare result plots for one-dimensional load spectra;

[0192] Fig.11 Prepare result plots for 2D load spectra;

[0193] Fig.12 This is the idea diagram of the wave center method;

[0194] Fig.13 This is the idea diagram of the variable mean method;

[0195] Fig.14 Customized product life result plots predicted by bilinear damage accumulation model and nonlinear damage accumulation model. DETAILED DESCRIPTION

[0196] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0197] Example:

[0198] like Figures 1 to 14 As shown, the present invention provides a multi-domain and multi-level collaborative simulation and life prediction method for complex customized products, comprising the following steps:

[0199] S1. Establish a finite element model of customized products;

[0200] The specific implementation steps are: establish a customized product model, simplify the model features to improve its simulation calculation efficiency, and establish a customized product finite element model;

[0201] S2, import the model for system dynamics simulation for subsequent model evaluation;

[0202] The specific implementation steps are: establish a system dynamics simulation model and complete the system dynamics simulation calculation.

[0203] S3. Import the model for modal analysis and obtain the modal result related files;

[0204] The specific implementation steps are: import the finite element model in S1, create a modal analysis simulation file and save its modal result file;

[0205] S4. Import the modal results into the dynamics part, establish and analyze the rigid-flexible coupling model of the customized product, and combine the system dynamics simulation results to obtain the performance evaluation of the customized product model and the time-load history used for subsequent life prediction;

[0206] The specific implementation steps are: import the modal analysis result file in S3, generate the rigid-flexible coupling fbi file, and import the rigid-flexible coupling model into the dynamic model;

[0207] S5, import the model and establish the fluid region for flow field analysis;

[0208] The specific implementation steps are:

[0209] S51, establishing an external flow field region, which is usually a rectangular parallelepiped;

[0210] S52, hollowing out the customized product model in the outer flow field area;

[0211] S53, setting turbulence model;

[0212] The customized product in this paper adopts the turbulent SST k-ω model, and its control equation is as follows:

[0213]

[0214]

[0215] Where: G k represents the turbulent kinetic energy, G ω is the ω equation, Γ k , Γ ω represent the effective diffusion terms of k and ω respectively, Y k , Y ω Represent the divergence terms of k and ω respectively. ω represents the orthogonal divergence term. S k With S ω is a user-defined term. The SST k-ω model is different from the standard k-ω model in that in the standard k-ω, α ∞ is a constant, while in the SST model, α ∞ The equation is as follows:

[0216] α ∞ =F1α ∞,1 +(1-F1)α ∞,2

[0217] Where:

[0218] SST k-ω is similar to the k-ω model, but with the following improvements:

[0219] (1) The deformation of the SST k-ω model and the k-ω model grows from the addition of the hybrid function and the dual model. The hybrid function is designed for the near-wall region, where the standard k-ω model is valid. There are also free surfaces, which are valid for the deformation of the k-ω model;

[0220] (2) The SST k-ω model incorporates the cross-diffusion from the ω equation;

[0221] (3) Turbulent viscosity takes into account the propagation of turbulent shear stress;

[0222] (4) The model constants are different.

[0223] These modifications make the SST k-ω model more accurate and reliable than the standard k-ω model in a wide range of flow fields.

[0224] S54, set boundary conditions U, P, k, nut, omega, etc.;

[0225] S55. Set the algorithm type to SIMPLE.

[0226] S56, solving the equation by using the NS equation;

[0227] The common form of the NS equation is as follows:

[0228]

[0229] Where: is the volume force; is the pressure differential force; For the viscosity force.

[0230] S57, completing the solution to obtain the external flow field pressure result of the customized product;

[0231] S6. Obtain the fluid analysis results and combine them with the finite element model of the customized product to calculate the fluid pressure as the boundary condition of the surface of the customized product in the fluid-solid coupling;

[0232] The specific implementation steps are:

[0233] S61, processing the flow field simulation result data file for subsequent calculation;

[0234] S62, calculating fluid node pressure;

[0235] Since the flow field solver uses the finite volume method for calculation, the calculated pressure is a scalar stored at the center of the volume unit. To perform fluid-solid coupling mapping, the pressure value of the fluid node must be calculated first.

[0236] In the flow field, the flow field nodes are mainly divided into two categories: nodes located on the flow field boundary; nodes located inside the flow field. The calculation of the flow field node pressure value is different for the two types of nodes:

[0237] (1) Nodes located on the flow field boundary:

[0238]

[0239] Where n2 is the node located on the boundary of the flow field; bf1 and bf7 are the values ​​of the boundary surfaces of units c1 and c7 respectively.

[0240] (2) Nodes located inside the flow field:

[0241]

[0242] Where: n1 is the node inside the flow field; c1-c7 is the unit value; V c1 -V c7 is the unit volume.

[0243] S63, calculating and obtaining the node pressure of the envelope surface of the solid external flow field;

[0244] S64, searching for the fluid node closest to the point to be interpolated;

[0245] All solid surface nodes are points to be interpolated. Before interpolation, we must first find the fluid node closest to the point to be interpolated, and use the octree method to search.

[0246] The octree method mainly performs search iterations by giving a given depth, and reduces the search range through continuous iterations, and finally searches for the required results;

[0247] S65, mapping of points to be interpolated;

[0248] According to the nearest fluid node found in S64, all surface elements containing the node are found, and the point to be interpolated is mapped to the plane to which one of the surface elements belongs:

[0249] Define the space plane equation as:

[0250] Ax+By+Cz+D=0

[0251] Where x, y, z are spatial coordinates, and A, B, C, D are equation coefficients.

[0252] Then substitute the three points of the surface unit into the spatial plane equation and solve the homogeneous linear equations to obtain the equation coefficients:

[0253] Ax1+By1+Cz1+D=0

[0254] Ax2+By2+Cz2+D=0

[0255] Ax3+By3+Cz3+D=0

[0256]

[0257] Where: x1, x2, x3 are the x coordinates of the three points, and the same goes for y and z.

[0258] After obtaining A, B, C, and D, the coordinates of the mapping point of the interpolation point in the plane are obtained according to the coordinate substitution formula:

[0259]

[0260]

[0261]

[0262] The three-dimensional triangle unit is converted to a two-dimensional coordinate system, and the interpolation points are mapped to the triangle unit plane;

[0263] Then we need to determine whether the mapping point of the interpolation point is within the surface unit:

[0264] Suppose a point P in space x is a point on the plane where the triangular unit P0P1P2 is located. Since the sides P0P1 and P0P2 of ΔP0P1P2 are not collinear, if the vector and As the basis, we know from linear algebra theory that the vector It can be expressed as:

[0265]

[0266] Where: i (i=1,2) is the unknown coefficient.

[0267] Multiply both sides of the above equation by The equations are:

[0268]

[0269] Solving this equation yields λ i (i=1,2), if λ is satisfied i satisfy:

[0270]

[0271] This indicates that point P x is within ΔP0P1P2, otherwise it indicates that the point is outside ΔP0P1P2.

[0272] S66, shape function interpolation of the point to be interpolated

[0273] According to S65, the mapping point coordinates (x p ,y p ,z p ), according to the invariance of spatial distance, the coordinates of triangles and mapping points are converted into a two-dimensional coordinate system.

[0274] Assume that the three points of the triangle unit are P0, P1, P2, and the mapping point is P x Define P0 as the origin P0'(0,0) in the two-dimensional coordinate system, then the two-dimensional coordinates of P1' are (0, P0P1), and the coordinates of P2' can be defined as follows:

[0275] P2'P1'=P2P1

[0276] P2'P0'=P2P0

[0277] The P2' coordinates are obtained by solving the above equations.

[0278] After completing the conversion of three-dimensional coordinates to two-dimensional coordinates, the triangle shape function is defined as:

[0279]

[0280] a i =x j y m -x m y j

[0281] b i =y j -y m

[0282] c i =x m -x j

[0283] Where: A is the area of ​​the triangle; x j 、x m are the x coordinates of two points other than point i, and y j ,y m is the y coordinate.

[0284] According to the shape function, point P x Pressure is defined as:

[0285]

[0286] Where: n is the number of nodes; N i is the shape function of each point; P i For each point pressure.

[0287] S67, Solving Solid Boundary Loading

[0288] The solid node pressure obtained in S66 cannot be directly applied to the solid boundary as a load and requires some processing.

[0289] In this paper, the solid unit is defined as a hexahedral unit, and each of its surface units is a quadrilateral surface unit. Four Gaussian points are defined in the quadrilateral surface unit. The pressure value at the Gaussian point can be obtained by shape function interpolation. The shape function of the quadrilateral is defined as:

[0290]

[0291]

[0292]

[0293]

[0294] Where: N i (i=1,2,3,4) is the shape function of each point of the quadrilateral; s and r are local coordinates.

[0295] Then define the solid boundary node loading as:

[0296]

[0297] Where: w s 、w r is the Gaussian point weight coefficient, which can be obtained by looking up the Gaussian point table; S is the area of ​​the quadrilateral; N i is the shape function of each point; P i is the pressure at each Gaussian point; J is the Jacobian matrix.

[0298] S7, loading fluid-solid coupling boundary conditions to perform fluid-solid coupling analysis;

[0299] The specific implementation steps are: loading the node forces calculated in S6 onto the solid surface nodes, then generating the loading of the solid surface nodes through the program and completing the call calculation.

[0300] S8. Import modal analysis files for acoustic simulation;

[0301] The specific implementation steps are:

[0302] S81, import modal result file;

[0303] S82, creating a customized product model field point;

[0304] S83. Create customized product outer envelope;

[0305] S84, mapping the vibration result of the customized product onto its outer envelope surface;

[0306] S85, establishing boundary element simulation;

[0307] S86. Calculate and obtain the noise field of the customized product;

[0308] S9. Import the model for static analysis to obtain strength, displacement and other results;

[0309] The specific implementation steps are: import the finite element model, apply boundary conditions and complete static simulation to obtain strength and displacement results.

[0310] S10, introduce the time-load history and static strength simulation results, establish the rain flow counting method and various load spectrum compilation methods;

[0311] The specific implementation steps are:

[0312] S101, introduce the time-load history in S4 and the static simulation results in S9,

[0313] The four-point cycle counting method uses four consecutive vertices to extract stress cycles. First, four consecutive vertices are taken from the stress time history, and the stress values ​​are (S1, S2, S3, S4) at a time; the stress range of the outer cycle is ΔS0 = |S1-S4|.

[0314] (1) If ΔS1≤ΔS0, extract the inner layer cyclic stress S2→S3; and directly connect the outer layer vertices S1 and S4. Go to step (3).

[0315] (2) If ΔS1≤ΔS0, do not extract the cycle; go to step (3);

[0316] (3) Select 4 consecutive vertices from the remaining vertices and recalculate until all loops are extracted.

[0317] Based on the above principle of rainflow counting method, it is used to extract the amplitude and mean of stress cycle.

[0318] After completing the stress cycle extraction, the load spectrum can be compiled. This article uses four load spectrum compilation methods:

[0319] (1) One-dimensional load spectrum compilation

[0320] First, the load-time history is scanned, and the load cycle and its maximum amplitude are obtained by the rain flow counting method. Since the mean value has little effect on fatigue life, the one-dimensional load spectrum does not consider the influence of the mean value.

[0321] Define the distance between groups: D = (F amax -F amin ) / N

[0322] Where D is the group spacing, Famax and Famin are the maximum and minimum values ​​of the load amplitude, and N is the number of load amplitude groups, which is usually defined as 8, 16, or 32.

[0323] After calculating the group spacing, the upper and lower limits of the amplitude are calculated based on the group spacing:

[0324]

[0325] After determining the upper and lower limits of each level, the load amplitude is graded. The frequency of the group in which the amplitude is located is increased by 1.

[0326] (2) Preparation of two-dimensional load spectrum

[0327] The two-dimensional load spectrum is also called the mean-amplitude spectrum. It introduces the influence of the mean value on fatigue life on the basis of the one-dimensional load spectrum, making the calculation results more accurate.

[0328] The two-dimensional load spectrum compilation introduces the grouping of mean values ​​on the basis of one-dimensional load spectrum compilation, and defines the formula for the mean group spacing:

[0329]

[0330] Calculate the upper and lower limits for each group:

[0331]

[0332] After the calculation is completed, determine which group the load mean and amplitude belong to and add 1 to the frequency of that group.

[0333] (3) Fluctuation center method

[0334] The fluctuation center method simplifies the mean and amplitude of the load-time history of the structure into a univariate random variable and integrates all of the two. The specific method is: average all the load cycle means and use it as the fluctuation center, making it exist as the static part, and then superimpose the amplitude as the dynamic part on the static part. The load data obtained in this way is used to compile the load spectrum to simplify the distribution of the mean.

[0335] The wave center method takes the average of all loads and calculates the total average value:

[0336]

[0337] Among them, σ mi is the mean value of each load cycle, V i is the frequency of each mean.

[0338] (4) Variable Mean Method

[0339] The wave center method is generally applicable to symmetrical load-time histories. For asymmetrical load-time histories, the variable mean method is usually used.

[0340] The variable mean method first groups the load cycles according to the amplitude, and then calculates the average value of each group. The average value calculation formula is:

[0341]

[0342] After the load spectrum is recompiled, fatigue life calculation can be performed.

[0343] S11. Based on the results of load spectrum compilation, various damage models are established for damage prediction, so as to carry out customized product life prediction.

[0344] S102. Establish multiple damage models and complete life prediction;

[0345] (1) Linear damage accumulation

[0346] The linear damage accumulation model is also known as the Miner-Palmgren theory, or Miner damage theory for short. It refers to the continuous accumulation of fatigue damage to a structure when subjected to a cyclic load of a certain amplitude. The damage generated each time is assumed to be unrelated to the stress received, and other factors are not considered, and a linear relationship is directly accumulated. Under a set of loads with varying amplitudes, the damage caused to the structure by n cyclic loads is:

[0347]

[0348] Where n is the number of independent stress cycles; i is the number of times a single stress acts; N i is the fatigue limit life under a single stress.

[0349] (2) Bilinear damage accumulation

[0350] The bilinear damage accumulation formula proposed by Manson-Grover in 1981 is as follows:

[0351] N Ⅰ =Nexp(ZN φ )

[0352] N Ⅱ =NN Ⅰ

[0353] Among them, z and is the model constant, and the calculation formula is as follows:

[0354]

[0355]

[0356] The damage calculation formulas for each stage under multi-level loads are as follows:

[0357]

[0358] (3) Nonlinear damage accumulation

[0359] The exponential loss model proposed by Dolan and Carten in 1956 requires the following assumptions:

[0360] 1) A crack nucleation period is required to produce permanent fatigue damage;

[0361] 2) As the stress gradually increases, the number of crack cores in the specimen will increase;

[0362] 3) The amount of damage increases with the increase of stress;

[0363] 4) With the increase of stress, the damage growth rate of each cycle increases;

[0364] 5) For all load processes, the total amount of damage that causes fatigue failure of the structure is constant;

[0365] 6) When the stress level is less than the limit value of fatigue damage, the damage continues to accumulate.

[0366] Under a given stress s, the accumulated fatigue damage D is expressed by an exponential relationship:

[0367] D=mrn a

[0368] Wherein, the number of fatigue crack nuclei is represented by m, the damage growth coefficient related to the stress level is represented by r, the number of stress cycles is represented by n, and the exponential constant related to the stress level is represented by a.

[0369] When two levels of stress s1 and s2 of different magnitude act on the specimen in an alternating manner, assuming that stress s1 undergoes an cycles, then (1-a)n is used to represent the number of cycles of stress s2. The estimated formula for fatigue life under two levels of stress is:

[0370]

[0371] Extending the above two-level stress fatigue life estimation formula to multi-level loads, we get:

[0372]

[0373] Where N g is the total fatigue life; N1 is the fatigue life under stress s1; each stress level s i The ratio of the number of cycles under the condition to the total number of cycles is expressed as ai Indicates; the maximum stress amplitude under multi-level alternating loads is represented by s1, r is the stress level, d is the material constant, which is usually determined by a two-level load test. Corten and Dolan recommend taking 4.8 as the normal value.

[0374] The present invention can quickly realize a multi-domain and multi-level collaborative simulation and life prediction method based on customized products. Starting from the operation status of customized products, it comprehensively considers the physical domains and corresponding environmental parameters that need to be calculated for customized products, and can more accurately judge the operation status of customized products in each physical domain.

[0375] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principle of the present invention should be included in the protection scope of the present invention.

Claims

1. A multi-domain and multi-level collaborative simulation and life prediction method for complex customized products, characterized in that: The method comprises the following steps: S1. Establish a finite element model of customized products; S2, import the model for system dynamics simulation for subsequent model evaluation; S3. Import the model for modal analysis and obtain the modal result related files; S4. Import the modal results into the dynamics part, establish and analyze the rigid-flexible coupling model of the customized product, and combine the system dynamics simulation results to obtain the performance evaluation of the customized product model and the time-load history used for subsequent life prediction; S5, import the model and establish the fluid region for flow field analysis; S6. Obtain the fluid analysis results and combine them with the finite element model of the customized product to calculate the fluid pressure as the boundary condition of the surface of the customized product in the fluid-solid coupling; S7, load the fluid-structure coupling boundary conditions and import them into LiToSim for fluid-structure coupling analysis; S8. Import modal analysis files for acoustic simulation; S9. Import the model for static analysis to obtain strength and displacement results; S10, introduce the time-load history and static strength simulation results, establish the rain flow counting method and various load spectrum compilation methods; S11. According to the results of load spectrum compilation, multiple damage models are established for damage prediction, so as to carry out customized product life prediction; The specific implementation method of step S6 is: S61, fluid analysis result data processing; Combined with the fluid analysis result data, processing is performed, including fluid unit connection query, polygonal surface unit splitting into triangle units, and boundary information processing; S62, fluid node pressure calculation; Fluid analysis is based on the finite volume method. The calculated pressure results are stored at the center of the volume unit. To transfer the pressure from the fluid to the solid surface, the pressure at the center of the fluid unit needs to be converted to the fluid node. The specific conversion steps are as follows; (1) First, extract the information of the fluid envelope solid outer surface; (2) Determine whether the fluid nodes on the envelope surface are located at the boundaries of other fluids; (3) If the answer is yes, then traverse all fluid surface units containing the node and determine whether the surface unit is a boundary surface unit; (4) If the answer is yes, a simple average calculation is performed on all boundary surface units containing fluid nodes to obtain the pressure value of the fluid node; (5) If (2) is judged as no, then directly traverse the fluid boundary surface units containing the fluid nodes and perform simple averaging; (6) Interpolation to obtain the solid surface node pressure; The interpolation method uses mapping-shape function interpolation. The mapping-shape function interpolation is based on the boundary information and fluid node pressure obtained in S62. The solid node pressure calculation steps are as follows:

1. Determine the points to be interpolated on the solid surface nodes; Second, use the octree method to traverse all nodes of the fluid envelope surface and find the nearest point to the specific interpolation point; 3. Given the distance error ε, take ε = element_size * 10 -5 , determine the distance between the point to be interpolated and the nearest node and compare it with ε; 4. If the distance in step 3 is less than ε, the pressure at the point to be interpolated is directly assigned to the value of the fluid node; 5. If the distance in step 3 is greater than ε, first search for all fluid surface units to which the fluid node belongs; 6. Select a surface unit and use the spatial plane equation to map the points to be interpolated to the plane where the surface unit is located: The spatial plane equation is: Ax+By+Cz+D=0 Substitute the coordinates of the three nodes of the surface element into the equation to solve the homogeneous linear equation system: Ax1+By1+Cz1+D=0 Ax2+By2+Cz2+D=0 Ax3+By3+Cz3+D=0 After solving A, B, C, and D, the coordinates of the interpolation point in the plane are obtained according to the coordinate substitution formula: In the formula, x p ,y p 、z p is the coordinate of the point to be interpolated in the plane, x0, y0, z0 are the original coordinates of the point to be interpolated; 7. Convert the mapping points and surface units to a two-dimensional coordinate system; The transformation is performed based on the invariance of spatial distance; assuming that the three points of the surface unit are P0, P1, and P2, and P0 is defined as the origin P0'(0, 0) in the two-dimensional coordinate system, then the two-dimensional coordinates of P1' are (0, P0P1), and the coordinates of P2' can be defined as follows: P2'P1'=P2P1 P2'P0'=P2P0 According to the above equations, the coordinates of P2' are obtained; 8. Determine whether the mapping point is within the face unit, if so, continue, otherwise return to step 6 to select the next face unit; 9. Constructing triangle element shape functions in surface elements: to i =x j and m -x m and j b i =and j -and m c i =x m -x j Where A is the area of ​​the triangle; x j 、x m are the x coordinates of two points other than point i, and y j ,y m is the y coordinate; 10. Interpolate the mapping points: According to the shape function, point P x Pressure is defined as: Where: n is the number of nodes; N i is the shape function of each point; P i is the pressure at each point; When solving solid boundary loading, the obtained solid node pressure cannot be directly applied as a load to the solid boundary; If the solid unit is defined as a hexahedral unit, then each of its face units is a quadrilateral face unit. Four Gaussian points are defined in the quadrilateral face unit. The pressure value at the Gaussian point can be obtained by shape function interpolation. The shape function of the quadrilateral is defined as: Where: Ni,i=1,2,3,4 are the shape functions of each point of the quadrilateral; s, r are local coordinates; Then define the solid boundary node loading as: Where: w s 、w r is the Gaussian point weight coefficient, which can be obtained by looking up the Gaussian point table; S is the area of ​​the quadrilateral; N i is the shape function of each point; P i is the pressure at each Gaussian point; J is the Jacobian matrix.

2. The multi-domain and multi-level collaborative simulation and life prediction method for complex customized products according to claim 1, characterized in that: In the step S1, a customized product model is first established to perform customized product segmentation.

3. The multi-domain and multi-level collaborative simulation and life prediction method for complex customized products according to claim 1 is characterized in that: The specific implementation method of step S2 is: establishing a multi-rigid body dynamics model and performing system dynamics simulation.

4. The multi-domain and multi-level collaborative simulation and life prediction method for complex customized products according to claim 1 is characterized in that: The specific implementation method of step S3 is: importing the customized product finite element model and establishing modal analysis.

5. The multi-domain and multi-level collaborative simulation and life prediction method for complex customized products according to claim 1 is characterized in that: The specific implementation method of step S7 is: loading according to the solid surface node force calculated in S6, and forming a command stream file for calculation.

Images