A vibration simulation method for ceramic matrix composites considering heat and damage comprehensively

Abstract

The application discloses a ceramic matrix composite material vibration simulation method considering heat and damage comprehensively, uses a finite element analysis method to perform dynamic analysis on CMC material, considers damage based on a woven ceramic matrix composite material arbitrary loading and unloading constitutive model, discretizes a temperature field and solves each node temperature by using an interpolation method to obtain each point thermal strain, solves a dynamic equation by using a central difference method in a time domain, can calculate the response of CMC material at each moment under different excitation frequencies, and obtains an amplitude-frequency curve, so that accurate and efficient simulation measurement of CMC high-temperature nonlinear dynamic response is realized.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of CMC material dynamics simulation, and particularly relates to a ceramic matrix composite material vibration simulation method considering heat and damage comprehensively. BACKGROUND

[0002] As one of the most important core components of an airplane, most components of an aero-engine are subjected to huge thermal load and force load in a harsh environment of high temperature and high frequency vibration for a long time. Future aerospace vehicles will fly at more than 5 Mach, and aerodynamic heating will quickly heat the vehicle body to more than 1000 DEG C, and the local working temperature of the engine will exceed 1700 DEG C, so a new type of high-temperature-resistant material is urgently needed to meet the requirements of the new generation of engines. Composite material vibration over-limit, i.e. resonance, flutter and other undesirable vibrations, will cause engine structure failure. Therefore, under ultra-high temperature conditions, whether the vibration performance of the new high-temperature-resistant composite material can meet the requirements of the new generation of aero-engines is one of the key factors to determine whether the new high-temperature-resistant composite material is qualified. Therefore, it is necessary to simulate the vibration performance of the composite material under ultra-high temperature conditions, which is of great significance to the further design of high-performance aero-engines.

[0003] At present, most of the simulation of the vibration characteristics of CMC materials is based on linear theory to calculate the natural frequency and damping of CMC materials, and there is little simulation of the nonlinear dynamics response of CMC materials. Due to the complexity of the weaving damage mechanism, the dynamics response of CMC materials lacks damage calculation, and the influence of thermal strain on the dynamics characteristics of CMC materials under high temperature conditions is also rarely considered.

[0004] Therefore, it is necessary to design a CMC vibration simulation method considering heat and damage to fill the gap in the field of vibration simulation of CMC composite materials under ultra-high temperature conditions. SUMMARY

[0005] The purpose of the present application is to provide a ceramic matrix composite material vibration simulation method considering heat and damage comprehensively, which can use the loading and unloading constitutive model of CMC material to consider damage in the dynamics calculation of CMC material, change the different temperature field provided by the user at a finite number of time into a temperature field meeting the time interval Δt in the central difference method, and calculate the thermal strain caused by temperature at different times, so as to consider heat field and damage in the dynamics calculation at the same time.

[0006] To achieve the above technical purpose, the technical scheme adopted by the present application is as follows:

[0007] A ceramic matrix composite material vibration simulation method considering heat and damage comprehensively, the method comprising the following steps:

[0008] S1, using ANSYS software to establish a finite element model of ceramic matrix composites, output unit and node information; calculate the unit geometry matrix B, mass matrix M and damping matrix C that do not change with time;

[0009] S2, input node temperature information, use the quadratic difference method to average the provided temperature field to the temperature field that satisfies the time interval Δt;

[0010] S3, calculate the mass matrix M and damping matrix C based on the element geometry matrix B;

[0011] S4, calculate the stress load K at time t t ;

[0012] S5, input external load information and calculate the external load matrix Q t ;

[0013] S6, calculate the node displacement vector a at the next moment t+Δt ;

[0014] S7, calculate the total strain vector of the 8 Gaussian integration points in each element

[0015] S8, calculate the thermal strain vector of 8 Gaussian integration points in each unit

[0016] S9, calculate the elastic strain vectors of 8 Gaussian integration points in each unit

[0017] S10, based on the CMCs constitutive model, calculate the stress vector σ at the next moment at the eight Gaussian integration points of each unit t+Δt ;

[0018] S11, determine whether the number of calculations has been reached, if so, proceed to step S12, otherwise set t = t + Δt and proceed to step S4;

[0019] S12, calculate the displacement response of the CMC material at each moment under different excitation frequencies, obtain the amplitude-frequency curve, and write the displacement, strain, and stress files at the corresponding moment.

[0020] Furthermore, the geometry matrix B, mass matrix M, damping matrix C, stress load matrix K t and the external load matrix Q t The expression is:

[0021]

[0022] Where, Ω e is the unit calculation domain, Γe is the unit external force domain, is the unit stress at the current moment, q tis the external force load at the current time, c is the damping coefficient, N e is the shape function of the i-th node of the element, ρ is the material density, G e is the conversion matrix of the element node degrees of freedom and the global degrees of freedom, B e is the element geometry matrix, the superscript T is the transpose symbol, the subscript t represents the current time, and the superscript e represents the element.

[0023] Further, in step S7, the total strain vector of the eight Gauss integral points in each element is calculated by the following formula

[0024]

[0025] In the formula, is the shape function of the i-th node of the element, is the displacement vector of the i-th node of the element, B e is the element geometry matrix, n is the total number of element nodes, is the element displacement vector, is the derivative with respect to the x direction, is the derivative with respect to the y direction, the superscript e represents the element, and the subscript t+Δt represents the next time.

[0026] Further, in step S8, the thermal strain vector of the eight Gauss integral points in each element is calculated by the following formula

[0027]

[0028] In the formula, is the instantaneous thermal expansion coefficient of the element, is the temperature change of the element, the superscript e represents the element, and the subscript t+Δt represents the next time.

[0029] Further, in step S9, the elastic strain vector of the eight Gauss integral points in each element is calculated by the following formula

[0030]

[0031] In the formula, B e is the element geometry matrix, is the temperature change of the element, is the thermal expansion coefficient of the element in each direction, is the temperature change of the element, the superscript e represents the element, and the subscript t+Δt represents the next time.

[0032] Further, the CMCs constitutive model is:

[0033]

[0034] where p u1 、p u2 、p u3 、p u4 、p r1 、p r2 、p r3 and p r4 All are about the maximum strain ε max function; ε is the strain; is the strain rate, which is used to judge whether it is a loading or unloading process; σ uploud is the stress during the loading process, σ reloud is the stress during the unloading process.

[0035] Furthermore, the central difference method is used to calculate the node displacement vector a at the next moment. t+Δt :

[0036] a t+Δt =M -1 (K t -Q t )Δt 2 +2a t -a t-Δt

[0037] Where a is the displacement, M is the mass matrix, K t is the stress load, Q t is the external load matrix, where the starting condition of the central difference method is:

[0038]

[0039] Where a0 is the initial displacement condition, is the initial velocity condition, is the initial acceleration condition, and the convergence conditions are as follows:

[0040]

[0041] Where Δt is the time step, Δt cr is the minimum time step, ω n is the highest order natural frequency of the vibration system, T n is the minimum period of the vibration system; T n The calculation formula is:

[0042]

[0043] Where, is the minimum element size, E is the elastic modulus of the material, and ρ is the material density.

[0044] Compared with the prior art, the present invention has the following beneficial effects:

[0045] First, the present invention's vibration simulation method for ceramic-based composites, which comprehensively considers heat and damage, incorporates the damage of CMC materials into dynamic calculations and substitutes the damage of CMC into the simulation calculations from a macroscopic phenomenological loading and unloading constitutive model, which is more consistent with experimental data. Compared with the previous method of quantifying damage per unit time, the use of a fitted macroscopic phenomenological loading and unloading constitutive model to consider its damage is closer to reality and more consistent with experimental data.

[0046] Second, the present invention's vibration simulation method for ceramic-based composites, which comprehensively considers heat and damage, incorporates thermal strain into the calculation of the dynamic response of CMC materials. To cope with the different temperature fields at different times provided by the user, the quadratic difference method is used to decompose the temperature field into temperature fields that satisfy the time interval Δt, so that in each iteration, the dynamic response calculation considering the thermal field under transient conditions can be analyzed. BRIEF DESCRIPTION OF THE DRAWINGS

[0047] Figure 1 This is a flow chart of the ceramic matrix composite material vibration simulation method that comprehensively considers heat and damage of the present invention;

[0048] Figure 2 Schematic diagram of the constitutive model of ceramic matrix composites at room temperature;

[0049] Figure 3 Schematic diagram of time domain displacement obtained from dynamic frequency sweep experiment on ceramic matrix composites;

[0050] Figure 4 for Figure 3 A magnified image of the resonance peak in the mesodynamic sweep experiment;

[0051] Figure 5 This is a time domain diagram of the displacement solution for the simulation at 0.5s around the resonance peak;

[0052] Figure 6 Schematic diagram of the displacement time domain curve obtained by simulation when the frequency sweep starts at 0.1s (150Hz-150.1Hz);

[0053] Figure 7 Schematic diagram of the displacement time domain curve obtained by simulation when the frequency sweep ends at 0.1s (239.9Hz-240Hz). DETAILED DESCRIPTION

[0054] The embodiments of the present invention are described in further detail below with reference to the accompanying drawings.

[0055] The present invention discloses a vibration simulation method for ceramic matrix composites that comprehensively considers heat and damage. The finite element analysis method is used to perform dynamic analysis on CMC materials. Damage is considered based on the arbitrary loading and unloading constitutive model of woven ceramic matrix composites. The temperature field is discretized and the interpolation method is used to solve the thermal strain at each point where the temperature of each node is solved. The central difference method is used in the time domain to solve the dynamic equation. The response of the CMC material at each moment under different excitation frequencies can be calculated, and the amplitude-frequency curve can be obtained, thereby realizing accurate and efficient simulation measurement of the high-temperature nonlinear dynamic response of the CMC.

[0056] The simulation calculation steps of the ceramic matrix composite material vibration simulation method of the present invention are as follows: Figure 1 The detailed description is as follows:

[0057] S1, using ANSYS software to establish a finite element model of ceramic matrix composites, output unit and node information; calculate the unit geometry matrix B, mass matrix M and damping matrix C that do not change with time;

[0058] S2, input node temperature information, use the quadratic difference method to average the provided temperature field to the temperature field that satisfies the time interval Δt;

[0059] S3, calculate the mass matrix M and damping matrix C based on the element geometry matrix B;

[0060] S4, calculate the stress load K at time t t ;

[0061] S5, input external load information and calculate the external load matrix Q t ;

[0062] As a preferred example, the geometry matrix B, mass matrix M, damping matrix C, stress load matrix K t and the external load matrix Q t The expression is:

[0063]

[0064] Where, Ω e is the unit calculation domain, Γ e is the unit external force domain, is the unit stress, q t is the external load, c is the damping coefficient, N e The shape function of the i-th node of the unit, ρ is the material density, G e B is the conversion matrix between unit node degrees of freedom and overall degrees of freedom, e is the unit geometry matrix, the superscript T is the transpose symbol, the superscript e represents the unit, and the subscript t+Δt represents the next moment.

[0065] S6, calculate the node displacement vector a at the next moment t+Δt ;

[0066] S7, calculate the total strain vector of the 8 Gaussian integration points in each element

[0067] In step S7, the total strain vector of the eight Gaussian integration points in each unit is calculated using the following formula:

[0068]

[0069] Where, is the shape function of the i-th node of the unit, is the displacement vector of the i-th node of the unit, B e is the unit geometry matrix, n is the total number of unit nodes, is the unit displacement vector, To find the derivative in the x direction, To find the derivative in the y direction, the superscript e represents the unit, and the subscript t+Δt represents the next moment.

[0070] S8, calculate the thermal strain vector of 8 Gaussian integration points in each unit

[0071] In step S8, the following formula is used to calculate the thermal strain vector of the eight Gaussian integral points in each unit:

[0072]

[0073] Where, is the unit instantaneous thermal expansion coefficient, is the unit temperature change, the superscript e represents the unit, and the subscript t+Δt represents the next moment.

[0074] S9, calculate the elastic strain vectors of 8 Gaussian integration points in each unit

[0075] In step S9, the elastic strain vectors of the eight Gaussian integration points in each unit are calculated using the following formula:

[0076]

[0077] Where B e is the unit geometry matrix, is the unit temperature change, is the thermal expansion coefficient of the unit in all directions, is the unit temperature change, the superscript e represents the unit, and the subscript t+Δt represents the next moment.

[0078] S10, calculating the stress vector of each unit at the next moment at 8 Gauss integral points based on the CMCs constitutive model t+Δt ;

[0079] The CMCs constitutive model is as follows:

[0080]

[0081] Wherein p u1 , p u2 , p u3 , p u4 , p r1 , p r2 , p r3 and p r4 are functions of the maximum strain epsilon max ; epsilon is a strain; is a strain rate, used to determine whether it is a loading or unloading process; sigma uploud is a stress in the loading process, sigma reloud is a stress in the unloading process. As shown in formula (1), it is a constitutive model of ceramic matrix composites at room temperature. Figure 2

[0082] S11, determining whether the calculation number is reached, if yes, turning to step S12, otherwise t=t+Delta t and turning to step S4;

[0083] S12, calculating the displacement response of the CMC material at each moment under different excitation frequencies, and obtaining an amplitude-frequency curve, and writing a displacement, strain and stress file corresponding to the moment.

[0084] In the present application, the method for calculating t=t+Delta t, i.e. the stress and strain vector of the next moment, is the central difference method, and the iterative formula is as follows:

[0085] a t+Δt =M -1 (K t -Q t )Delta t 2 +2a t -a t-Δt

[0086] The starting condition of the above central difference method is as follows:

[0087]

[0088] The time step Delta t of the central difference method determines the stability of the calculation. When Delta t is less than the critical Delta t cr , the calculation result converges, otherwise it diverges, so the convergence condition of the central difference method is as follows:

[0089] ​

[0090] ω n is the highest order natural frequency of the vibration system, T n is the minimum period of the vibration system. n The calculation formula is:

[0091]

[0092] The minimum unit size is .

[0093] The dynamic sweep frequency experiment (150Hz-240Hz) of the ceramic matrix composite flat piece under the acceleration excitation of 10m / s 2 at high temperature 200℃ was carried out, and the corresponding time domain displacement graph was obtained, as shown in Figure 3 , the peak part was amplified to obtain Figure 4 .

[0094] According to the above theoretical calculation method, the time domain displacement graph of the ceramic matrix composite under the acceleration excitation of 10m / s 2 at high temperature 200℃ was calculated by programming with c++. Due to the large amount of transient dynamics calculation, the calculation time is limited, and the displacement time domain graph of 0.5s around the resonance peak was calculated, as shown in Figure 5 . The displacement time domain curves at the beginning and end of the sweep frequency were also calculated, i.e. the sweep frequency of 150-150.1Hz and 239.9-240Hz in the experiment, as shown in Figure 6 and Figure 7 .

[0095] Those skilled in the art should understand that the embodiments of the present application can be provided as a method, a system, or a computer program product. Therefore, the present application can adopt a complete hardware embodiment, a complete software embodiment, or an embodiment combining software and hardware aspects. Moreover, the present application can adopt the form of a computer program product implemented on one or more computer usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) containing computer usable program code. The solutions in the embodiments of the present application can be implemented in various computer languages, such as object-oriented programming language Java and interpreted scripting language JavaScript.

[0096] The present application is described with reference to the flowcharts and / or block diagrams of the methods, devices (systems), and computer program products according to the embodiments of the present application. It should be understood that each process and / or box in the flowchart and / or block diagram, as well as the combination of the processes and / or boxes in the flowchart and / or block diagram, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, a special-purpose computer, an embedded processor, or other programmable data processing device to produce a machine, so that the instructions executed by the processor of the computer or other programmable data processing device generate instructions for implementing the steps in the process. Figure 1 a process or multiple processes and / or boxes Figure 1 A device that provides the functions specified in a block or multiple blocks.

[0097] These computer program instructions may also be stored in a computer readable memory that can direct a computer or other programmable data processing device to work in a specific manner, so that the instructions stored in the computer readable memory produce an article of manufacture comprising an instruction device, which implements the process Figure 1 a process or multiple processes and / or boxes Figure 1 The function specified in one or more boxes.

[0098] These computer program instructions can also be loaded onto a computer or other programmable data processing device so that a series of operating steps are executed on the computer or other programmable device to produce a computer-implemented process, thereby providing the instructions for executing the steps in the process. Figure 1 a process or multiple processes and / or boxes Figure 1 A step that specifies a function in one or more boxes.

[0099] Although the preferred embodiments of the present application have been described, those skilled in the art may make additional changes and modifications to these embodiments once they have learned the basic creative concept. Therefore, the appended claims are intended to be interpreted as including the preferred embodiments and all changes and modifications that fall within the scope of the present application.

[0100] Obviously, those skilled in the art may make various changes and modifications to this application without departing from the spirit and scope of this application. Thus, if these modifications and variations of this application fall within the scope of the claims of this application and their equivalents, this application is intended to include these modifications and variations.

Claims

1. A vibration simulation method for ceramic matrix composite materials that comprehensively considers heat and damage, characterized in that: The method The following steps are involved: S1, using ANSYS software to establish a finite element model of ceramic matrix composites, output unit and node information; calculate the unit geometry matrix B, mass matrix M and damping matrix C that do not change with time; S2, input node temperature information, use the quadratic difference method to average the provided temperature field to the temperature field that satisfies the time interval Δt; S3, calculate the mass matrix M and damping matrix C based on the element geometry matrix B; S4, calculate the stress load K at time t t ; S5, input external load information and calculate the external load matrix Q t ; S6, calculate the node displacement vector a at the next moment t+Δt ; S7, calculate the total strain vector of the 8 Gaussian integration points in each element S8, calculate the thermal strain vector of 8 Gaussian integration points in each unit S9, calculate the elastic strain vectors of 8 Gaussian integration points in each unit S10, based on the CMCs constitutive model, calculate the stress vector σ at the next moment at the eight Gaussian integration points of each unit t+Δt ; S11, determine whether the number of calculations has been reached, if so, proceed to step S12, otherwise set t = t + Δt and proceed to step S4; S12, calculate the displacement response of the CMC material at each moment under different excitation frequencies, obtain the amplitude-frequency curve, and write the displacement, strain, and stress files at the corresponding moment; In step S7, the total strain vector of the eight Gaussian integration points in each unit is calculated using the following formula: Where, is the shape function of the i-th node of the unit, is the displacement vector of the i-th node of the unit, B e is the unit geometry matrix, n is the total number of unit nodes, is the unit displacement vector, To find the derivative in the x direction, To find the derivative in the y direction, the superscript e represents the unit and the subscript t+Δt represents the next moment; In step S8, the following formula is used to calculate the thermal strain vector of the eight Gaussian integral points in each unit: Where, is the unit instantaneous thermal expansion coefficient, is the unit change temperature, the superscript e represents the unit, and the subscript t+Δt represents the next moment; In step S9, the elastic strain vectors of the eight Gaussian integration points in each unit are calculated using the following formula: Where B e is the unit geometry matrix, is the unit temperature change, is the thermal expansion coefficient of the unit in all directions, is the unit temperature change, the superscript e represents the unit, and the subscript t+Δt represents the next moment.

2. The ceramic matrix composite material vibration simulation method considering heat and damage according to claim 1 is characterized in that: Geometry matrix B, mass matrix M, damping matrix C, stress load matrix K t and the external load matrix Q t The expression is: Where, Ω e is the unit calculation domain, Γ e is the unit external force domain, is the unit stress at the current moment, q t is the external force load at the current moment, c is the damping coefficient, N e The shape function of the i-th node of the unit, ρ is the material density, G e B is the conversion matrix between unit node degrees of freedom and overall degrees of freedom, e is the unit geometry matrix, the superscript T is the transpose symbol, the subscript t represents the current moment, and the superscript e represents the unit.

3. The ceramic matrix composite material vibration simulation method considering heat and damage according to claim 1, characterized in that: The constitutive model of CMCs is: where p u1 、p u2 、p u3 、p u4 、p r1 、p r2 、p r3 and p r4 All are about the maximum strain ε max function; ε is the strain; is the strain rate, which is used to determine whether it is a loading or unloading process; σ uploud is the stress during the loading process, σ reloud is the stress during the unloading process.

4. The ceramic matrix composite material vibration simulation method considering heat and damage according to claim 1, characterized in that: Use the central difference method to calculate the node displacement vector a at the next moment t+Δt : a t+Δt =M -1 (K t -Q t )Δt 2 +2a t -a t-Δt Where a is the displacement, M is the mass matrix, K t is the stress load, Q t is the external load matrix, where the starting condition of the central difference method is: Where a0 is the initial displacement condition, is the initial velocity condition, is the initial acceleration condition, and the convergence conditions are as follows: Where Δt is the time step, Δt cr is the minimum time step, ω n is the highest order natural frequency of the vibration system, T n is the minimum period of the vibration system; T n The calculation formula is: Where, is the minimum element size, E is the elastic modulus of the material, and ρ is the material density.

Images