CN109033514B - Method for evaluating elastic instability of plane tube bundle fluid - Google Patents
Method for evaluating elastic instability of plane tube bundle fluid Download PDFInfo
- Publication number
- CN109033514B CN109033514B CN201810626314.0A CN201810626314A CN109033514B CN 109033514 B CN109033514 B CN 109033514B CN 201810626314 A CN201810626314 A CN 201810626314A CN 109033514 B CN109033514 B CN 109033514B
- Authority
- CN
- China
- Prior art keywords
- matrix
- unit
- node
- fluid
- tube
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02E—REDUCTION OF GREENHOUSE GAS [GHG] EMISSIONS, RELATED TO ENERGY GENERATION, TRANSMISSION OR DISTRIBUTION
- Y02E30/00—Energy generation of nuclear origin
- Y02E30/30—Nuclear fission reactors
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Complex Calculations (AREA)
- Investigating Strength Of Materials By Application Of Mechanical Stress (AREA)
Abstract
The invention provides a plane tube bundle fluid elastic instability assessment method which comprises two parts of modal analysis and Connors formula calculation. The modal analysis part is responsible for calculating the natural frequency and the corresponding vibration mode of the heat exchange tube; and calculating the equivalent flow rate and the critical flow rate under different modes according to the Connors formula. The invention has the following advantages: the modal analysis calculation of different outside-tube fluid parameters, nonuniform transverse flow velocity and complex heat exchange tube shape can be realized by utilizing the discrete characteristics of the finite element; on the basis of finite element modal analysis results, a Connors formula is calculated by adopting a Gaussian integration method, so that the weighting effect of the equivalent flow velocity of the vibration mode is realized.
Description
Technical Field
The invention relates to a method for evaluating the elastic instability of a plane tube bundle fluid.
Background
The shell-and-tube heat exchanger is widely applied to the fields of electric power and chemical industry, and mainly utilizes a heat exchange tube bundle to realize large-area heat exchange of fluid media inside and outside the tube. The flow direction of the fluid outside the tubes can be split into longitudinal and transverse directions as it flows through the tube bundle. If the transverse flow velocity is too high, the fluid elasticity instability phenomenon can be generated, so that the heat exchange tube can vibrate greatly and finally be damaged and fail.
Fluid elastic destabilization is a fluid-solid strong coupling problem. The movement of each heat exchange tube affects the fluid force and movement of the other heat exchange tubes and creates self-excitation forces. For the occurrence of destabilization, it is believed that the fluid-solid coupling effects adversely affect structural mass, damping, and effort. Because of the complex mechanism, no very perfect theory exists at present, and most of formulas for predicting the occurrence of instability are empirical formulas obtained through experimental measurement.
Appendix C of heat exchanger standard GB/T151-2014 provides a series of calculation formulas and judgment criteria for fluid induced vibration calculation and check, but the application range is limited to the problem that single medium outside the tube uniformly flows through a straight tube or a U-shaped tube, and the support form of the heat exchange tube has clear requirements. Section N-1330 in ASME BPVC-III applications describes some of the features that result from fluid elastic instability, briefly gives the choice of recommended formulas and empirical constants, but does not give an explicit calculation of specific problems.
The following difficulties are often created using the above criteria:
1) The problem that multiple fluid media outside the pipe cannot be calculated;
2) A problem of inability to calculate non-uniform lateral flow rates;
3) The complex heat exchange tube shape and the variable span support form can not be calculated.
In practical heat exchanger designs, one or more of the above-described problems are often encountered. For example, a serpentine tube high-pressure heater in a conventional thermal power plant, wherein the shell-side outside fluid comprises three types of superheated steam, saturated steam and supercooled water, and belongs to the multi-medium problem; the flow rates of the superheating section, the condensing section and the cooling section are different, and the method belongs to the problem of non-uniform transverse flow rate; the heat exchange tubes are arranged in a serpentine shape, and belong to the problem of complicated heat exchange tube shapes.
The problems can only depend on past experience, and the design is very trouble.
Disclosure of Invention
The invention aims to provide a calculation method for quantitatively analyzing the problem of fluid elastic instability of a planar tube bundle with any shape under the condition of nonuniform transverse flow of various media.
In order to achieve the above purpose, the technical scheme of the invention is to provide a method for evaluating the elastic instability of a plane tube bundle fluid, which is characterized by comprising the following steps:
step 1, simplifying a heat exchange tube into a one-dimensional curve model in a three-dimensional space by utilizing the axis of the heat exchange tube; splitting the one-dimensional curve model into a plurality of line segments, and ensuring that the discrete line segment model is basically consistent with the original model in geometry.
A three-dimensional cartesian coordinate system is defined in space, and node coordinates, cell node numbers, node constraints are listed for subsequent computation. Wherein, define the discrete line segment as the unit, define the terminal point of the discrete line segment as the node. In order to achieve higher computational accuracy under a limited discrete model, a secondary unit is employed here, i.e. the nodes of the unit comprise two end points of the discrete line segments and the midpoint of the line segments, three in total.
The node coordinate list LOC is expressed as:
the unit node number list NOD is expressed as:
the node constraint list UROT is expressed as:
step 2, calculating a unit mass matrix M according to the Timoshenko beam theory, the finite element theory and the numerical integration method e Cell stiffness matrix K e ：
Wherein n is gaussm And n gaussk The number of gaussian integral points for the mass and stiffness arrays, respectively. For a secondary unit with intermediate nodes, n gaussm Taking 3, n gaussk Taking 2. Corresponding Gaussian point position xi i And weight w i The following are provided:
j e ＝l e and/2 is the standard domain ([ -1, 1)]One-dimensional space of (a) to a singleMeta local coordinates ([ -1) e /2，1 e /2]Is a one-dimensional space) Jacobian determinant, l e Is the length of the element line segment:
n (ζ) is a matrix function with respect to the standard domain coordinates ζ:
in the method, in the process of the invention,
is a shape function of the unit node 1;
N 3 ＝N 3 (ξ)＝1-ξ 2 is a shape function of the unit node 3;
matrix C is M e Cross-sectional correlation matrix of (a):
wherein ρ, A, J, I y ，I z The unit density, the sectional area, the torsion moment of inertia, the y-direction bending moment of inertia and the z-direction bending moment of inertia are respectively;
b (ζ) is a matrix function with respect to the standard domain coordinates ζ:
in the method, in the process of the invention,
matrix D is K e Cross-sectional correlation matrix of (a):
in E, K y 、G y 、K z 、G z 、K x G is Young's modulus, y-direction shear correction coefficient, y-direction shear modulus, z-direction shear correction coefficient, z-direction shear modulus, torsion correction coefficient and shear modulus respectively;
the matrix T is a transformation matrix of the local coordinates of the cell ([ -le/2, one-dimensional space of le/2] to global coordinates (three-dimensional space):
in the method, in the process of the invention,
step 3, the unit quality matrix M calculated in the step 2 e Cell stiffness matrix K e Sequentially accumulating the node degrees of freedom numbers into an overall mass matrix M and an overall rigidity K matrix:
of the formula (I)
Rather than summing symbols, it represents the assembly of an overall mass, stiffness matrix;
at the time of assembly, M e And K e All elements in (a) are accumulated in M and K according to the following matrix rank number correspondence:
[(m-1)*3+p-1，(n-1)*3+q-1]
→[(NOD[e，m]-1)*3+p-1，(NOD[e，n]-1)*3+q-1]
where m, n (=1 to 3) are the unit node numbers, and p, q (=1 to 6) are the node degrees of freedom, respectively.
After assembly, according to table UROT, for constrained degrees of freedom, a large number (e.g., by 10 50 )。
Step 4, calculating to obtain the characteristic value and the characteristic vector of the integral mass matrix and the integral rigidity matrix:
the above equation is a generalized eigenvalue equation, where ω 2 Is a characteristic value of the characteristic,
is a feature vector. The calculation amount of equation solving is large, and the equation solving can be carried out through a computer. />
After the solution is completed, the circular frequency omega (or natural frequency) of each order mode can be obtained
) And displacement vectors of all nodes +.>
Step 5, calculating the equivalent flow velocity v under each order mode by adopting Connors formula under the condition that the flow velocity and the density of the along-pipe are different e 2 ：
Wherein ρ is 0 To average out-of-pipe fluid density ρ e For unit e out-of-pipe fluid density, v e For the lateral flow rate of the fluid outside the unit e tube, u j e Is the total displacement of the unit eNode j (extract step 4 results
Corresponding displacement value of (b), N j (ζ) definition in step 3, m 0 Is average mass per unit length, m t e Is the mass per unit length, w, of the unit e i And xi i Weighting and integrating point position for three-point Gaussian product (value synchronization step 3), j e Is Jacobian determinant (value step 3).
Step 6, calculating critical flow velocity v under each order mode by adopting natural frequency result c
Wherein C isIs Connors coefficient, ζ is damping ratio, D o And f is the natural frequency (the calculation result of the extraction step 4) for the outer diameter of the heat exchange tube.
Step 7: calculating the flow rate ratio of each order mode, namely equivalent flow rate/critical flow rate, and if all the flow rate ratios are smaller than 1, checking and passing the fluid elastic instability; if a certain flow rate ratio is greater than or equal to 1, the fluid elasticity instability check is not passed, and the original design can be modified by observing a vibration mode diagram.
Preferably, in the step 2, the calculation formula of the cell density ρ is:
wherein ρ is t Is the density of the heat exchange tube;
m i is the mass of fluid in a unit length of pipe, r i For the inner radius of the heat exchange tube, ρ i Is the fluid density in the tube; />
m o Is the mass of fluid outside the tube per unit length, r o For the outer radius of the heat exchange tube, ρ o For the density of the fluid outside the tube, C M Is an additional quality coefficient.
Preferably, in said step 2, K z 、K y The calculation formula of (2) is as follows:
wherein, v is poisson ratio.
The invention has the following advantages:
1) The modal analysis calculation of different outside-tube fluid parameters, nonuniform transverse flow velocity and complex heat exchange tube shape can be realized by utilizing the discrete characteristics of the finite element;
2) On the basis of finite element modal analysis results, a Connors formula is calculated by adopting a Gaussian integration method, so that the weighting effect of the equivalent flow velocity of the vibration mode is realized.
By combining the advantages, the calculation method can solve the three problems in the background art, and is suitable for calculating the fluid elasticity instability of the conventional shell-and-tube heat exchanger.
Detailed Description
The invention will be further illustrated with reference to specific examples. It is to be understood that these examples are illustrative of the present invention and are not intended to limit the scope of the present invention. Further, it is understood that various changes and modifications may be made by those skilled in the art after reading the teachings of the present invention, and such equivalents are intended to fall within the scope of the claims appended hereto.
The invention provides a plane tube bundle fluid elastic instability assessment method which comprises two parts of modal analysis and Connors formula calculation. The modal analysis part is responsible for calculating the natural frequency and the corresponding vibration mode of the heat exchange tube; and calculating the equivalent flow rate and the critical flow rate under different modes according to the Connors formula.
The modal analysis mechanical model is based on Timoshenko beam theory, and the numerical calculation method adopts a finite element algorithm. Basic assumptions include the theory of small deformation, linear elastic isotropy constitutive relation, timoshenko Liang Jiben assumption (beam cross section remains planar during bending).
In general, the modality analysis section includes the following steps 1 to 4:
step 1: converting the Timoshenko beam equation into a weak form equation, and solving domain dispersion by adopting a linear finite element method to obtain a generalized eigenvalue equation in a matrix form.
Step 2: the discrete units select a secondary beam unit with an intermediate node, and the corresponding unit mass matrix and the corresponding unit stiffness matrix are obtained through numerical integration by utilizing the corresponding shape function and constitutive equation. The numerical integration adopts a Gaussian integration method, the integration points of the mass array are 3, and the integration points of the stiffness array are 2.
Step 3: and assembling all the unit mass matrixes and the unit stiffness matrixes to obtain an overall mass matrix and an overall stiffness matrix.
Step 4: substituting the overall mass matrix and the overall stiffness matrix obtained in the step 3 into a generalized characteristic equation, and solving a characteristic vector (vibration mode) and a characteristic value (circular frequency).
The Connors formula calculation part includes the following steps 5 to 7:
step 5: and (3) adopting a Connors equivalent flow rate formula and a vibration mode result, and obtaining the equivalent flow rate after weighted average by numerical integration (Gaussian integration).
Step 6: and calculating the critical flow rate by adopting a Connors critical flow rate formula and a natural frequency result.
Step 7: the flow rate ratio (equivalent flow rate/critical flow rate) of each order mode is calculated. If all flow ratios are less than 1, the fluid elastic instability check passes. If a certain flow rate ratio is greater than or equal to 1, the fluid elastic instability check does not pass; at this time, the original design can be modified by observing the vibration pattern diagram.
The calculation step comprises a large number of iteration and matrix operation, and particularly, the algorithm of solving the generalized eigenvalue of the sparse matrix is needed in the step 4, and the calculation can be performed through a computer.
The modal analysis mainly converts the Tie-Skoch beam equation into a weak form equation, solves domain dispersion, and ignores external force, so that an equation in a matrix form can be obtained:
Mü+Ku＝0
where u is the displacement vector, u is the acceleration vector, and the general solution of the equation is:
in the method, in the process of the invention,
is a displacement characteristic vector, ω is a circular frequency, and t is time. Substituting the general solution form into an equation to obtain a generalized eigenvalue equation:
wherein the matrices M and K are as follows,
representing the assembly of the overall mass matrix and the overall stiffness matrix:
specifically, the modal analysis includes the following calculation steps:
step 1, simplifying the heat exchange tube into a one-dimensional curve model in a three-dimensional space by utilizing the axis of the heat exchange tube (section parameters are considered in unit calculation in step 2); splitting the one-dimensional curve model into a plurality of line segments, and ensuring that the discrete line segment model is basically consistent with the original model in geometry.
A three-dimensional cartesian coordinate system is defined in space, and node coordinates, cell node numbers, node constraints are listed for subsequent computation. Wherein, define the discrete line segment as the unit, define the terminal point of the discrete line segment as the node. In order to achieve higher computational accuracy under a limited discrete model, a secondary unit is employed here, i.e. the nodes of the unit comprise two end points of the discrete line segments and the midpoint of the line segments, three in total.
The node coordinate list LOC is expressed as;
the unit node number list NOD is expressed as:
the node constraint list UROT is expressed as:
Step 2, solving an integral expression of the unit mass/stiffness matrix by generally adopting Gaussian product. The main reason is that: 1) The manual deduction integral operation is avoided, so that the solving process of the unit rigidity/quality matrix is more universal; 2) The effect of volume/shear locking can be reduced by reducing the integration.
The gaussian product will also be used in the subsequent Connors algorithm, the correlation calculation being given in detail here.
For a general function f (x) defined over the standard domain Γ, its integral over the standard domain Γ may be approximately solved by a gaussian product equation:
in xi i Is the Gaussian point location on the standard domain Γ, W i Is a corresponding weighting.
There is a definition in the cell domain Ω e General function g above, define x: Γ.fwdarw.Ω e Then:
in the method, in the process of the invention,
is Jacobian determinant of transform x, j i ＝j(ξi)，g i ＝g(x(ξ i ))，x(ξ i ) For the transformation of the actual coordinates with the standard domain, g (x (ζ) i ) A relation of the function g with respect to the standard domain coordinates ζ.
Applying gaussian product to the integral expression of the mass/stiffness matrix can result in:
wherein n is gaussm And n gaussk The number of gaussian integral points for the mass and stiffness arrays, respectively. For a secondary unit with intermediate nodes, n gaussm Taking 3, n gaussk Taking 2. Corresponding Gaussian point position xi i And weight w i The following are provided:
j e ＝l e and/2 is the standard domain ([ -1, 1)]Is transformed to unit local coordinates ([ -1) e /2，1 e /2]Is a one-dimensional space) Jacobian determinant, l e Is the length of the element line segment:
n (ζ) is a matrix function with respect to the standard domain coordinates ζ:
in the method, in the process of the invention,
is a shape function of the unit node 1; />
N 3 ＝N 3 (ξ)＝1-ξ 2 is a shape function of the unit node 3;
matrix C is M e Cross-sectional correlation matrix of (a):
wherein ρ, A, J, I y ，I z The cell density, the cross-sectional area, the torsional moment of inertia, the y-bending moment of inertia, and the z-bending moment of inertia are respectively defined. The influence of the additional mass needs to be considered when the pipe vibrates in the fluid, and the calculation method is as follows:
wherein ρ is t Is the density of the heat exchange tube;
m i is the mass of fluid in a unit length of pipe, r i For the inner radius of the heat exchange tube, ρ i Is the fluid density in the tube; />
m o Is the mass of fluid outside the tube per unit length, r o For the outer radius of the heat exchange tube, ρ o For the density of the fluid outside the tube, C M Is an additional quality coefficient.
B (ζ) is a matrix function with respect to the standard domain coordinates ζ:
in the method, in the process of the invention,
matrix D is K e Cross-sectional correlation matrix of (a):
in E, K y 、G y 、K z 、G z 、K x G is Young's modulus, y-direction shear correction coefficient, y-direction shear modulus, z-direction shear correction coefficient, z-direction shear modulus, torsion correction coefficient, and shear modulus, respectively. For round tubes, torsion does not generate warpage, thus K x ＝1。
Wherein v is poisson ratio, r i R is the inner radius of the heat exchange tube o Is the outer radius of the heat exchange tube.
The matrix T is a transformation matrix of the local coordinates of the cell ([ -le/2, one-dimensional space of le/2] to global coordinates (three-dimensional space):
in the method, in the process of the invention,
step 3, the unit quality matrix M calculated in the step 2 e Cell stiffness matrix K e Sequentially accumulating the node degrees of freedom numbers into an overall mass matrix M and an overall rigidity K matrix:
of the formula (I)
Rather than summing symbols, it represents the assembly of an overall mass, stiffness matrix;
at the time of assembly, M e And K e All elements in (a) are accumulated in M and K according to the following matrix rank number correspondence:
[(m-1)*3+p-1，(n-1)*3+q-1]
→[(NOD[e，m]-1)*3+p-1，(NOD[e，n]-1)*3+q-1]
where m, n (=1 to 3) are the unit node numbers, and p, q (=1 to 6) are the node degrees of freedom, respectively.
After assembly, according to table UROT, for constrained degrees of freedom, a large number (e.g., by 10 50 )。
Step 4, calculating to obtain the characteristic value and the characteristic vector of the integral mass matrix and the integral rigidity matrix:
the above equation is a generalized eigenvalue equation, where ω 2 Is a characteristic value of the characteristic,
is a feature vector. The calculation amount of equation solving is large, and the equation solving can be carried out through a computer.
After the solution is completed, the circular frequency omega (or natural frequency) of each order mode can be obtained
) And displacement vectors of all nodes +.>
The Connors formula calculation includes steps 5 to 7:
step 5, calculating the equivalent flow velocity v under each order mode by adopting Connors formula under the condition that the flow velocity and the density of the along-pipe are different e 2 ：
Wherein ρ (x) is the density of the fluid outside the tube, v (x) is the lateral flow velocity outside the tube, mt (x) is the mass per unit length, ψ (x) is the displacement under the corresponding mode, and L is the total length of the heat exchange tube.
Because the formula includes integral calculations, the calculation is still performed using a gaussian product formula. The number of gaussian product points is taken to be 3. Zeta type toy i And W is i Local coordinates and weights of Gaussian integral points are respectively given; j (j) e Jacobian determinant for element e, of size L e /2。
Wherein ρ is 0 To average out-of-pipe fluid density ρ e For unit e out-of-pipe fluid density, v e For the lateral flow rate of the fluid outside the unit e tube, u j e Is the total displacement of the unit eNode j (extract step 4 results
Corresponding displacement value of (b), N j (ζ) definition in step 3, m 0 Is average mass per unit length, m t e Is the mass per unit length of unit e.
Step 6, calculating critical flow velocity V under each order mode by adopting natural frequency result c
Wherein C is a Connors coefficient, ζ is a damping ratio，D o And f is the natural frequency (the calculation result of the extraction step 4) for the outer diameter of the heat exchange tube.
Step 7: calculating the flow rate ratio of each order mode, namely equivalent flow rate/critical flow rate, and if all the flow rate ratios are smaller than 1, checking and passing the fluid elastic instability; if a certain flow rate ratio is greater than or equal to 1, the fluid elasticity instability check is not passed, and the original design can be modified by observing a vibration mode diagram.
Calculation example 1 (conventional problem): subject matter of input condition reference heat exchanger Standard GB/T151-2014 appendix C.7
Comparing the calculation result of the algorithm with the calculation result of GB/T151-2014 as follows:
the data show that the invention is consistent with the calculation result of GB/T151-2014 for the problem of fluid elasticity instability of the conventional heat exchange tube in which the medium uniformly flows outside the single tube.
For calculation of the problem of fluid elastic instability of heat exchange tubes of complex shape in heterogeneous flows of various media, reference can be made to calculation example 2.
Calculation example 2: and calculating the fluid elastic instability problem of the serpentine tube bundle under the heterogeneous flow of various out-of-tube media.
The calculation result of the algorithm is as follows:
in addition, the mode shape graphs of different modes can be obtained through mode analysis. For the mode with the flow rate ratio larger than 1, the position of the occurrence of the fluid elastic instability can be judged by observing the vibration mode diagram, and then the design modification is carried out on the weak area.
Claims (3)
1. A method for assessing the elastic instability of a fluid in a planar tube bundle, comprising the following steps:
step 1, simplifying a heat exchange tube into a one-dimensional curve model in a three-dimensional space by utilizing the axis of the heat exchange tube; splitting the one-dimensional curve model into a plurality of line segments, wherein the discrete line segment model is basically consistent with the original model in geometry;
defining a three-dimensional Cartesian coordinate system in a space, and listing node coordinates, unit node numbers and node constraints for subsequent calculation, wherein a discrete line segment is defined as a unit, the end point of the discrete line segment is defined as a node, and the node of the unit comprises two end points of the discrete line segment and a line segment midpoint, wherein the three end points are in total:
the node coordinate list LOC is expressed as:
the unit node number list NOD is expressed as:
the node constraint list UROT is expressed as:
step 2, calculating a unit mass matrix M according to the Timoshenko beam theory, the finite element theory and the numerical integration method e Unit and method for manufacturing sameStiffness matrix K e ：
Wherein n is gaussm And n gaussk The number of Gaussian integral points, ζ, of the mass array and the stiffness array respectively i Is Gaussian point position, w i Is weighted; j (j) e ＝l e 2 is Jacobian determinant of standard domain transformation to cell local coordinates, l e Is the length of the element line segment:
n (ζ) is a matrix function with respect to the standard domain coordinates ζ:
in the method, in the process of the invention,
is a shape function of the unit node 1;
N 3 ＝N 3 (ξ)＝1-ξ 2 is a shape function of the unit node 3;
matrix C is M e Cross-sectional correlation matrix of (a):
wherein ρ, A, J, I y ，I z The unit density, the sectional area, the torsion moment of inertia, the y-direction bending moment of inertia and the z-direction bending moment of inertia are respectively;
b (ζ) is a matrix function with respect to the standard domain coordinates ζ:
in the method, in the process of the invention,
matrix D is K e Cross-sectional correlation matrix of (a):
in E, K y 、G y 、K z 、G z 、K x G is Young's modulus, y-direction shear correction coefficient, y-direction shear modulus, z-direction shear correction coefficient, z-direction shear modulus, torsion correction coefficient and shear modulus respectively;
the matrix T is a transformation matrix of the local coordinates of the cell to global coordinates:
in the method, in the process of the invention,
step 3, the unit quality matrix M calculated in the step 2 e Cell stiffness matrix K e Sequentially accumulating the node degrees of freedom numbers into an overall mass matrix M and an overall rigidity matrix K:
of the formula (I)
Rather than summing symbols, it represents the assembly of an overall mass, stiffness matrix;
at the time of assembly, M e And K e All elements in (a) are accumulated in M and K according to the following matrix rank number correspondence:
[(m-1)*3+-1,(n-1)*3+-1]
→[(NOD[e,m]-1)*3+-1,(NOD[e,n]-1)*3+-1]
wherein m, n=1, 2,3 are the unit node numbers, p, q=1, 2,3, …,6 are the node degrees of freedom, respectively;
after the assembly is completed, setting a large number of main elements corresponding to M and K for the constrained degree of freedom according to a table UROT;
step 4, calculating to obtain the characteristic value and the characteristic vector of the integral mass matrix and the integral rigidity matrix:
the above equation is a generalized eigenvalue equation, where ω 2 Is a characteristic value of the characteristic,
is a characteristic vector, and the circular frequency omega or the natural frequency +.>
And displacement vectors of all nodes +.>
Step 5, calculating the equivalent flow velocity v under each order mode by adopting Connors formula under the condition that the flow velocity and the density of the along-pipe are different e 2 ：
Wherein: ρ 0 Is the average out-of-tube fluid density; ρ e The density of the fluid outside the unit e pipe; v e A lateral flow rate of fluid outside the cell e-tube; u (u) j e Extracting the result of the step 4 for the total displacement of the unit eNode j
Corresponding displacement values of (a); n (N) j (ζ) the definition in step 3; m is m 0 Is the average mass per unit length; m is m t e Mass per unit length of unit e; w (w) i And xi i A step 3 of weighting and integrating point position for three-point Gaussian product and value synchronization; j (j) e Taking a Jacobian determinant as a value in the step 3;
step 6, calculating critical flow velocity v under each order mode by adopting natural frequency result c
Wherein C is Connors coefficient, ζ is damping ratio, D o The outer diameter of the heat exchange tube is f, and the natural frequency is f;
step 7: calculating the flow rate ratio of each order mode, namely equivalent flow rate/critical flow rate, and if all the flow rate ratios are smaller than 1, checking and passing the fluid elastic instability; if a certain flow rate ratio is greater than or equal to 1, the fluid elasticity instability check is not passed, and the original design can be modified by observing a vibration mode diagram.
2. The method for evaluating the elastic destabilization of a fluid in a flat tube bundle according to claim 1, wherein in the step 2, a calculation formula of the cell density ρ is:
wherein ρ is t Is the density of the heat exchange tube;
m i is the mass of fluid in a unit length of pipe, r i For the inner radius of the heat exchange tube, ρ i Is the fluid density in the tube; />
m o Is the mass of fluid outside the tube per unit length, r o For the outer radius of the heat exchange tube, ρ o For the density of the fluid outside the tube, C M Is an additional quality coefficient.
3. A method for assessing the elastohydrodynamic instability of a tube bundle according to claim 2, wherein in said step 2, K z 、K y The calculation formula of (2) is as follows:
wherein, v is poisson ratio.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810626314.0A CN109033514B (en) | 2018-06-15 | 2018-06-15 | Method for evaluating elastic instability of plane tube bundle fluid |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810626314.0A CN109033514B (en) | 2018-06-15 | 2018-06-15 | Method for evaluating elastic instability of plane tube bundle fluid |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN109033514A CN109033514A (en) | 2018-12-18 |
| CN109033514B true CN109033514B (en) | 2023-04-28 |
Family
ID=64609499
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201810626314.0A Active CN109033514B (en) | 2018-06-15 | 2018-06-15 | Method for evaluating elastic instability of plane tube bundle fluid |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN109033514B (en) |
Families Citing this family (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111680378B (en) * | 2020-07-17 | 2022-09-16 | 天华化工机械及自动化研究设计院有限公司 | ANSYS-based heat exchanger tube bundle modal analysis method in liquid filling state |
| CN115114872B (en) * | 2022-07-20 | 2023-09-26 | 中国核动力研究设计院 | Parameter identification method and system for predicting tube bundle fluid bullet instability |
| CN115238494B (en) * | 2022-07-21 | 2023-10-20 | 中国核动力研究设计院 | Component position identification method for pipeline flow bullet instability |
Citations (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2006099549A (en) * | 2004-09-30 | 2006-04-13 | Keio Gijuku | Numerical calculation method for higher order mixed element poisson solver based on gsmac finite element method |
| CN101839470A (en) * | 2010-03-16 | 2010-09-22 | 上海电气电站设备有限公司 | U-shaped tube high pressure heater for nuclear power generating unit |
| CN102679772A (en) * | 2012-06-15 | 2012-09-19 | 张家港市江南锅炉压力容器有限公司 | Tube plate type heat exchanger |
| KR101192335B1 (en) * | 2011-05-16 | 2012-10-26 | 한국과학기술연구원 | Method for simulating fluid flow and recording medium for performing the method |
| CN104899387A (en) * | 2015-06-17 | 2015-09-09 | 西南石油大学 | Stream generator tube bundle two-phase transverse fluid elastic force unstability analysis method |
| WO2016069313A1 (en) * | 2014-10-27 | 2016-05-06 | Ebullient, Llc | Two-phase cooling system component |
| CN107025315A (en) * | 2016-02-02 | 2017-08-08 | 上海核工程研究设计院 | A kind of U-shaped heat-transfer pipe Flow vibration of nuclear power station steam generator and fretting wear coupling analysis computational methods |
| CN107480322A (en) * | 2017-06-23 | 2017-12-15 | 中国工程物理研究院总体工程研究所 | Free body multiple spot correlation pulse pressure random vibration analysis computational methods |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101256502B (en) * | 2007-02-27 | 2011-02-09 | 国际商业机器公司 | System and method for simulating multiprocessor system |
-
2018
- 2018-06-15 CN CN201810626314.0A patent/CN109033514B/en active Active
Patent Citations (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2006099549A (en) * | 2004-09-30 | 2006-04-13 | Keio Gijuku | Numerical calculation method for higher order mixed element poisson solver based on gsmac finite element method |
| CN101839470A (en) * | 2010-03-16 | 2010-09-22 | 上海电气电站设备有限公司 | U-shaped tube high pressure heater for nuclear power generating unit |
| KR101192335B1 (en) * | 2011-05-16 | 2012-10-26 | 한국과학기술연구원 | Method for simulating fluid flow and recording medium for performing the method |
| CN102679772A (en) * | 2012-06-15 | 2012-09-19 | 张家港市江南锅炉压力容器有限公司 | Tube plate type heat exchanger |
| WO2016069313A1 (en) * | 2014-10-27 | 2016-05-06 | Ebullient, Llc | Two-phase cooling system component |
| CN104899387A (en) * | 2015-06-17 | 2015-09-09 | 西南石油大学 | Stream generator tube bundle two-phase transverse fluid elastic force unstability analysis method |
| CN107025315A (en) * | 2016-02-02 | 2017-08-08 | 上海核工程研究设计院 | A kind of U-shaped heat-transfer pipe Flow vibration of nuclear power station steam generator and fretting wear coupling analysis computational methods |
| CN107480322A (en) * | 2017-06-23 | 2017-12-15 | 中国工程物理研究院总体工程研究所 | Free body multiple spot correlation pulse pressure random vibration analysis computational methods |
Non-Patent Citations (10)
| Title |
|---|
| 三段式高压加热器管束流体弹性失稳临界流速的计算;王莉华等;《力学季刊》;20100915(第03期);全文 * |
| 两相横向流诱发管束弹性不稳定性的试验研究进展;张亚楠等;《轻工机械》;20171231(第03期);全文 * |
| 影响管阵流体弹性不稳定性因素的研究;黄平等;《核动力工程》;19911028(第05期);全文 * |
| 换热器管束中的流体弹性不稳定性;聂清德，郭宝玉，丁学仁，金楠，陈旭;《力学学报》;19960319(第02期);全文 * |
| 换热器管束流体弹性不稳定性的预测;张明贤，聂清德，侯曾炎，王献旗;《石油化工设备》;19940425(第04期);全文 * |
| 管束的流体弹性不稳定性的研究;聂清德等;《振动工程学报》;19930630(第02期);全文 * |
| 管阵流体弹性不稳定性分析;吴一红等;《应用力学学报》;19930630(第02期);全文 * |
| 蒸汽发生器传热管流弹失稳计算;韩同行等;《压力容器》;20141130(第11期);全文 * |
| 螺旋管束流体弹性不稳定;薄涵亮，马昌文，佟允宪，姚梅生;《核科学与工程》;19950330(第01期);全文 * |
| 高压加热器管束振动机理简析;牛忠华等;《电站辅机》;20110330(第01期);全文 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN109033514A (en) | 2018-12-18 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109033514B (en) | Method for evaluating elastic instability of plane tube bundle fluid | |
| Mirzaeifar et al. | Coupled thermo-mechanical analysis of shape memory alloy circular bars in pure torsion | |
| CN112580236B (en) | Rapid analysis method for nonlinear dynamic response of thermal protection connection structure | |
| El-Bakari et al. | Assessing impact force localization by using a particle swarm optimization algorithm | |
| Gao et al. | Inverse identification of the mechanical parameters of a pipeline hoop and analysis of the effect of preload | |
| Zhang et al. | Weak form quadrature element analysis of planar slender beams based on geometrically exact beam theory | |
| Ji et al. | Numerical analysis of shell-side flow-induced vibration of elastic tube bundle in heat exchanger | |
| Albanesi et al. | Inverse finite element method for large‐displacement beams | |
| JP2019091316A (en) | Heat exchanger analysis method | |
| Di et al. | Eigen-parameter decomposition of element matrices for structural damage detection | |
| CN110807225B (en) | Transfer matrix calculation stability optimization method based on dimensionless analysis | |
| Chen et al. | Eigensolution reanalysis of modified structures using epsilon‐algorithm | |
| JP6497588B2 (en) | Deformation stress characteristic acquisition method and seismic evaluation method of heat exchanger | |
| Yan et al. | A comprehensive comparison on vibration and heat transfer of two elastic heat transfer tube bundles | |
| CN116306178A (en) | Structural strain inversion method based on self-adaptive shape function and equivalent neutral layer | |
| US20190293482A1 (en) | Self-excited vibration evaluation method | |
| Qimin et al. | The optimal design and simulation of helical spring based on particle swarm algorithm and MATLAB | |
| CN113408040A (en) | Analog data correction method and system in civil engineering | |
| JP4731233B2 (en) | Computer-implemented method for determining structural dimensions | |
| CN115048846B (en) | Model order reduction and stability judgment method and system for tube bundle fluid bullet system | |
| JP3946699B2 (en) | Region shape optimization method | |
| Tang et al. | Flow-induced vibration of curved cylinder arrays subject to loose support | |
| Elbanhawy et al. | Simulations of Fully-Flexible Fuel Bundle Response due to Turbulence Excitation | |
| CN112179169B (en) | Temperature control heat exchange method of three-fluid heat exchanger | |
| Ji et al. | Numerical Calculation Method |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |