CN109033514A - A kind of flat tube beam fluid elastic instability assessment method - Google Patents
A kind of flat tube beam fluid elastic instability assessment method Download PDFInfo
- Publication number
- CN109033514A CN109033514A CN201810626314.0A CN201810626314A CN109033514A CN 109033514 A CN109033514 A CN 109033514A CN 201810626314 A CN201810626314 A CN 201810626314A CN 109033514 A CN109033514 A CN 109033514A
- Authority
- CN
- China
- Prior art keywords
- formula
- matrix
- node
- unit
- fluid
- 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.)
- Granted
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
A kind of flat tube beam fluid elastic instability assessment method provided by the invention includes that model analysis and Connors formula calculate two parts.It is responsible for calculating the intrinsic frequency and the corresponding vibration shape of heat exchanger tube in model analysis part；Connors formula calculates the equivalent flow velocity and critical flow velocity then calculated according to modal analysis result under different modalities.The present invention has the advantage that the model analysis that different extratubal fluid parameters, non-homogeneous transverse flow speed, Complex Heat tube shape may be implemented calculates using the discrete feature of finite element itself；On the basis of finite element modal analysis result, Connors formula is calculated using Gauss integration method, realizes the vibration shape to the weighting effect of equivalent flow velocity.
Description
Technical field
The present invention relates to the methods that a kind of pair of flat tube beam fluid elastic instability is evaluated.
Background technique
Shell-and-tube heat exchanger is widely used in electric power, chemical field, mainly realizes that outer fluid is situated between in managing using heat-exchanging tube bundle
The large area heat exchange of matter.For extratubal fluid when flowing through tube bank, flow direction can be analyzed to vertical and horizontal.If transverse flow
Speed is excessive, then can generate elasticity of fluid unstable phenomenon, and heat exchanger tube is caused substantially to vibrate and final damage inactivation.
Elasticity of fluid unstability is a kind of solid close coupling problem of stream.The movement of every heat exchanger tube will affect fluid force and its
The movement of his heat exchanger tube, and generate self excitation force.Generation for unstable phenomenon, it is considered that be fluid structure interaction to structure matter
Amount, damping and active force produce adverse effect.It is most of pre- currently without very perfect theory due to its complicated mechanism
The formula that dendrometry surely occurs all is the empirical equation obtained by test measurement.
The appendix C of heat exchanger standard GB/T 151-2014 gives a system to the calculating and check of fluid-induced vibration
Column count formula and decision criteria, but use scope is only limitted to the problem of outer Single Medium of pipe flows uniformly through straight tube or U-tube, it is right
The supporting form of heat exchanger tube clearly requires.N-1330 chapters and sections in ASME BPVC-III Appendices describe fluid
Some features that elastic instability generates, briefly give the selection of recommended formula and empirical, but not to particular problem
Provide specific calculation method.
Following difficulty is often generated using above-mentioned standard:
1) can not computer tube outer multiple fluid medium the problem of；
2) the problem of non-homogeneous transverse flow speed can not be calculated；
3) Complex Heat tube shape can not be calculated, become the problem of span supporting form.
In actual design of heat exchanger, one or more above problems are often encountered.Such as in conventional thermal power plant
Coiled pipe high-pressure heater, the extratubal fluid of shell-side includes superheated steam, saturated vapor and three kinds of subcooled water, belongs to multimedium
Problem；Superheat section, condense section and dredge cold section flow velocity it is different, belong to non-homogeneous transverse flow speed problem；Heat exchanger tube is snakelike
Arrangement, belong to Complex Heat tube shape problem.
Many puzzlements can only be caused in design by past use experience by encountering the above problem.
Summary of the invention
The object of the present invention is to provide a kind of tube banks of quantitative analysis arbitrary shape plane in the non-homogeneous transverse flow of a variety of media
The calculation method of elasticity of fluid destabilization problems under dynamic.
In order to achieve the above object, the technical solution of the present invention is to provide a kind of evaluations of flat tube beam fluid elastic instability
Method, which comprises the following steps:
Step 1, the axis using heat exchanger tube, one-dimensional curve model heat exchanger tube being reduced in three-dimensional space；It will be one-dimensional
If curve model is split as main section, discrete lines segment model and master mould must be guaranteed geometrically almost the same.
Three dimensional cartesian coordinates system is defined in space, and node coordinate, cell node number, joint constraint are arranged
Table, so as to subsequent calculating.Wherein, discrete line segment is defined as unit, the endpoint of discrete line segment is defined as node.In order to have
Reach higher computational accuracy under the discrete model of limit, uses two sub-cells here, i.e. the node of unit includes discrete line segment
Two endpoints and line segment midpoint, totally three.
Node coordinate list LOC table is shown as:
In formula, xi、yi、ziThe x-axis of respectively i-th node, y-axis, z-axis coordinate；
Cell node numbered list NOD is indicated are as follows:
In formula,The number of 1st, 2,3 node of respectively e-th unit, the 1st, 2,3 of e-th of unit
Node is two endpoints and an intermediate node of e-th of unit；
Joint constraint list UROT is indicated are as follows:
In formula, uxi、uyi、uzi、rotxi、rotyi、rotzi6 freedom degrees of respectively i-th node constrain；
Step 2, according to Timoshenko beam theory, finite element theory and numerical integration method, computing unit moment of mass
Battle array MeAnd element stiffness matrix Ke:
In formula, ngaussmAnd ngausskThe respectively Gauss integration point quantity of Mass matrix and Stiffness Matrix.For band intermediate node
Two sub-cells, ngaussmTake 3, ngausskTake 2.Corresponding Gauss point position ξiWith weighting wiIt is as follows:
je=le/ 2 are converted into unit local coordinate ([- 1 for canonical domain (one-dimensional space of [- 1,1])e/ 2,1e/ 2] one
Dimension space) Jacobian determinant, leIt is the length of unit line segment:
N (ξ) is the matrix function about canonical domain coordinate ξ:
In formula,For the shape function of cell node 1；
For the shape function of cell node 2；
N3=N3(ξ)=1- ξ2For the shape function of cell node 3；
Matrix C is MeSection correlation matrix:
In formula, ρ, A, J, Iy, IzRespectively cell density, sectional area, torsional moment of inertia, y are to bending the moment of inertia, z to bending
The moment of inertia；
B (ξ) is the matrix function about canonical domain coordinate ξ:
In formula,
Matrix D is KeSection correlation matrix:
In formula, E, Ky、Gy、Kz、Gz、Kx, G be respectively Young's modulus, y to Splice variant, y to modulus of shearing, z to cutting
Cut correction factor, z to modulus of shearing, torsion correction factor, modulus of shearing；
Matrix T is unit local coordinate (one-dimensional space of [- le/2, le/2]) to the transformation of world coordinates (three-dimensional space)
Matrix:
In formula,
Step 3, the element mass matrix M that will be calculated in step 2eAnd element stiffness matrix KeAccording to the volume of the degree of freedom on a node basis
It number successively adds up into total quality matrix M and overall stiffness K matrix:
In formulaIt not sums symbol, but represents the assembling of total quality, stiffness matrix；
In assembling, by MeAnd KeIn all elements be added in M and K according to following matrix position No. corresponding relationship:
[(m-1) * 3+p-1, (n-1) * 3+q-1]
→ [(NOD [e, m] -1) * 3+p-1, (NOD [e, n] -1) * 3+q-1]
Wherein m, n (=1~3) are respectively cell node number, and p, q (=1~6) are respectively the degree of freedom on a node basis.
After being completed, according to table UROT, for restrained freedom degree, the corresponding pivot of M and K set big number (such as
Multiply 1050)。
Step 4, the characteristic value and feature vector that total quality matrix and Bulk stiffness matrix is calculated:
Above formula is generalized eigenvalue equation, wherein ω2It is characteristic value,It is feature vector.The calculation amount of equation solution compared with
Greatly, it can be solved by computer.
Circular frequency ω (or the intrinsic frequency of every first-order modal can be obtained after the completion of solving) and all nodes
Motion vector
Step 5, in the case of different from density along pipe flow velocity, calculated under every first-order modal using Connors formula
Equivalent flow velocity ve 2:
In formula, ρ0For average extratubal fluid density, ρeFor unit e extratubal fluid density, veFor the cross of unit e extratubal fluid
To flow velocity, uj eFor total displacement (4 result of extraction step of unit e node jIn corresponding displacement value), Nj(ξ) is the same as in step 3
Definition, m0For mean unit linear mass, mt eFor the linear mass of unit e, wiAnd ξiFor 3 Gaussian quadratures weighting with
Point position (the same step 3) of value, jeFor Jacobian determinant (the same step 3) of value.
Step 6 calculates critical flow velocity v under every first-order modal using intrinsic frequency resultc
In formula, C is Connors coefficient, and ζ is damping ratio, DoFor exchange heat pipe outside diameter, f be intrinsic frequency (extraction step 4
Calculated result).
Step 7: calculate the velocity ratio of every first-order modal, i.e., equivalent flow velocity/critical flow velocity, if all velocity ratios less than 1,
Then elasticity of fluid unstability check passes through；It is more than or equal to 1 if there is some velocity ratio, then elasticity of fluid unstability is checked obstructed
It crosses, can be modified at this time by observing bending vibation mode picture to original design.
Preferably, in the step 2, the calculation formula of cell density ρ are as follows:
In formula, ρtFor heat exchanger tube density；miFor unit length tube fluid quality, riIt is in heat exchanger tube half
Diameter, ρiFor tube fluid density；moFor unit length extratubal fluid quality, roFor heat exchanger tube outer radius, ρo
For extratubal fluid density, CMFor mass coefficient.
Preferably, in the step 2, Kz、KyCalculation formula are as follows:
In formula, υ is Poisson's ratio.
The present invention has the advantage that
1) the discrete feature for utilizing finite element itself may be implemented different extratubal fluid parameters, non-homogeneous transverse flow speed, answer
The model analysis of miscellaneous heat exchange tube shape calculates；
2) on the basis of finite element modal analysis result, Connors formula is calculated using Gauss integration method,
The vibration shape is realized to the weighting effect of equivalent flow velocity.
In conjunction with above-mentioned advantage, so that this calculation method is able to solve the three classes problem mentioned in background technique, it is suitable for mesh
Preceding common shell-and-tube heat exchanger elasticity of fluid unstability calculates.
Specific embodiment
Present invention will be further explained below with reference to specific examples.It should be understood that these embodiments are merely to illustrate the present invention
Rather than it limits the scope of the invention.In addition, it should also be understood that, after reading the content taught by the present invention, those skilled in the art
Member can make various changes or modifications the present invention, and such equivalent forms equally fall within the application the appended claims and limited
Range.
A kind of flat tube beam fluid elastic instability assessment method provided by the invention includes model analysis and Connors formula
Calculate two parts.It is responsible for calculating the intrinsic frequency and the corresponding vibration shape of heat exchanger tube in model analysis part；Connors formula calculates then
The equivalent flow velocity and critical flow velocity under different modalities are calculated according to modal analysis result.
Model analysis mechanical model is based on Timoshenko beam theory, and numerical computation method uses finite element algorithm.Substantially
Assuming that including small deformation theory, linear elasticity isotropism constitutive relation, (beam section had been bent Timoshenko beam basic assumption
Plane is still kept in journey).
In general, model analysis part includes the following steps 1 to step 4:
Step 1: converting weak form equation for Timoshenko Beam equation, using linear finite method, domain will be solved
It is discrete, obtain the generalized eigenvalue equation of matrix form.
Step 2: discrete unit selects the secondary beam element with intermediate node, using corresponding shape function and constitutive equation,
Numerical integration obtains corresponding element mass matrix and element stiffness matrix.Numerical integration uses Gaussian quadrature method, Mass matrix
Point be 3, the point of Stiffness Matrix is 2.
Step 3: all element mass matrix and element stiffness matrix being assembled, total quality matrix and whole is obtained
Body stiffness matrix.
Step 4: the total quality matrix and Bulk stiffness matrix that step 3 is acquired substitute into Generalized Characteristic Equation, solve special
Levy vector (vibration shape) and characteristic value (circular frequency).
Connors formula calculating section then includes the following steps 5 to step 7:
Step 5: using the equivalent velocity formula of Connors and the vibration shape as a result, numerical integration (Gauss integration) to obtain weighting flat
Equivalent flow velocity after.
Step 6: using Connors critical flow velocity formula and intrinsic frequency as a result, calculating critical flow velocity.
Step 7: calculating the velocity ratio (equivalent flow velocity/critical flow velocity) of every first-order modal.If all velocity ratios less than 1,
Then elasticity of fluid unstability check passes through.It is more than or equal to 1 if there is some velocity ratio, then elasticity of fluid unstability is checked obstructed
It crosses；It can be modified at this time by observing bending vibation mode picture to original design.
Above-mentioned calculating step includes that a large amount of iteration and matrix operation, especially step 4 are needed using sparse matrix broad sense spy
The algorithm that value indicative solves, can be calculated by computer.
Model analysis mainly utilizes timoshenko beam is equations turned for weak form equation, it is discrete will to solve domain, and ignore
External force, then the equation of available matrix form:
M ü+Ku=0
In formula, u is motion vector, and ü is vector acceleration, the general solution form of equation are as follows:
In formula,For displacement characteristic vector, ω is circular frequency, and t is the time.General solution form is substituted into equation, obtains broad sense
Eigenvalue equation:
Wherein matrix M and K are as follows,Represent the assembling of total quality matrix and Bulk stiffness matrix:
Specifically, model analysis includes following calculating step:
Heat exchanger tube is reduced to the one-dimensional curve model in three-dimensional space and (joined in section by step 1, the axis using heat exchanger tube
Number will account in the calculating of step 2 unit)；If one-dimensional curve model is split as main section, it must guarantee discrete line segment mould
Type and master mould are geometrically almost the same.
Three dimensional cartesian coordinates system is defined in space, and node coordinate, cell node number, joint constraint are arranged
Table, so as to subsequent calculating.Wherein, discrete line segment is defined as unit, the endpoint of discrete line segment is defined as node.In order to have
Reach higher computational accuracy under the discrete model of limit, uses two sub-cells here, i.e. the node of unit includes discrete line segment
Two endpoints and line segment midpoint, totally three.
Node coordinate list LOC table is shown as；
In formula, xi、yi、ziThe x-axis of respectively i-th node, y-axis, z-axis coordinate；
Cell node numbered list NOD is indicated are as follows:
In formula,The number of 1st, 2,3 node of respectively e-th unit, the 1st, 2,3 of e-th of unit
Node is two endpoints and an intermediate node of e-th of unit；
Joint constraint list UROT is indicated are as follows:
In formula, uxi、uyi、uzi、rotxi、rotyi、rotzi6 freedom degrees of respectively i-th node constrain.
Step 2 generallys use Gaussian quadrature to solve element quality/stiffness matrix integral expression.Main cause exists
In: 1) it avoids deriving integral operation manually, keeps element stiffness/mass matrix solution procedure more pervasive；2) contracting can be passed through
The mode of debulk point reduces the influence of volume/shear locking.
Gaussian quadrature can be also used in subsequent Connors algorithm, provides relevant calculation in detail here.
For the generic function f (x) being defined on canonical domain Γ, the integral on canonical domain Γ can be asked by Gauss
Product formula approximate solution:
In formula, ξiIt is the Gauss point position on canonical domain Γ, WiIt is to weight accordingly.
Existing one is defined on unit domain ΩeOn generic function g, define x: Γ → Ωe, then:
In formula,It is the Jacobian determinant for converting x, ji=j (ξ i), gi=g (x (ξi)), x
(ξi) be actual coordinate and canonical domain transformation relation, g (x (ξi)) it is relationship of the function g about canonical domain coordinate ξ.
Gaussian quadrature is applied to quality/stiffness matrix integral expression, available:
In formula, ngaussmAnd ngausskThe respectively Gauss integration point quantity of Mass matrix and Stiffness Matrix.For band intermediate node
Two sub-cells, ngaussmTake 3, ngausskTake 2.Corresponding Gauss point position ξiWith weighting wiIt is as follows:
je=le/ 2 are converted into unit local coordinate ([- 1 for canonical domain (one-dimensional space of [- 1,1])e/ 2,1e/ 2] one
Dimension space) Jacobian determinant, leIt is the length of unit line segment:
N (ξ) is the matrix function about canonical domain coordinate ξ:
In formula,For the shape function of cell node 1；
For the shape function of cell node 2；
N3=N3(ξ)=1- ξ2For the shape function of cell node 3；
Matrix C is MeSection correlation matrix:
In formula, ρ, A, J, Iy, IzRespectively cell density, sectional area, torsional moment of inertia, y are to bending the moment of inertia, z to bending
The moment of inertia.Need to consider the influence of additional mass when pipe vibrates in a fluid, calculation method is as follows:
In formula, ρtFor heat exchanger tube density；miFor unit length tube fluid quality, riIt is in heat exchanger tube half
Diameter, ρiFor tube fluid density；moFor unit length extratubal fluid quality, roFor heat exchanger tube outer radius,
ρoFor extratubal fluid density, CMFor mass coefficient.
B (ξ) is the matrix function about canonical domain coordinate ξ:
In formula,
Matrix D is KeSection correlation matrix:
In formula, E, Ky、Gy、Kz、Gz、Kx, G be respectively Young's modulus, y to Splice variant, y to modulus of shearing, z to cutting
Cut correction factor, z to modulus of shearing, torsion correction factor, modulus of shearing.For round tube, torsion can't generate warpage, therefore
Kx=1.
In formula, υ is Poisson's ratio, riFor heat exchanger tube inside radius, roFor heat exchanger tube outer radius.
Matrix T is unit local coordinate (one-dimensional space of [- le/2, le/2]) to the transformation of world coordinates (three-dimensional space)
Matrix:
In formula,
Step 3, the element mass matrix M that will be calculated in step 2eAnd element stiffness matrix KeAccording to the volume of the degree of freedom on a node basis
It number successively adds up into total quality matrix M and overall stiffness K matrix:
In formulaIt not sums symbol, but represents the assembling of total quality, stiffness matrix；
In assembling, by MeAnd KeIn all elements be added in M and K according to following matrix position No. corresponding relationship:
[(m-1) * 3+p-1, (n-1) * 3+q-1]
→ [(NOD [e, m] -1) * 3+p-1, (NOD [e, n] -1) * 3+q-1]
Wherein m, n (=1~3) are respectively cell node number, and p, q (=1~6) are respectively the degree of freedom on a node basis.
After being completed, according to table UROT, for restrained freedom degree, the corresponding pivot of M and K set big number (such as
Multiply 1050)。
Step 4, the characteristic value and feature vector that total quality matrix and Bulk stiffness matrix is calculated:
Above formula is generalized eigenvalue equation, wherein ω2It is characteristic value,It is feature vector.The calculation amount of equation solution compared with
Greatly, it can be solved by computer.
Circular frequency ω (or the intrinsic frequency of every first-order modal can be obtained after the completion of solving) and all nodes
Motion vector
Connors formula is calculated including step 5 to step 7:
Step 5, in the case of different from density along pipe flow velocity, calculated under every first-order modal using Connors formula
Equivalent flow velocity ve 2:
In formula, ρ (x) is extratubal fluid density, and v (x) is to manage outer transverse flow speed, and mt (x) is unit linear mass, and ψ (x) is
Displacement under corresponding mode, L are heat exchanger tube overall length.
Due to containing integral calculation in formula, still calculated here using Gauss quadrature formula.Gaussian quadrature point
Quantity take 3.ξiAnd WiThe respectively local coordinate and weighting of Gauss integration point；jeFor the Jacobian determinant of unit e,
Size is Le/2。
In formula, ρ0For average extratubal fluid density, ρeFor unit e extratubal fluid density, veFor the cross of unit e extratubal fluid
To flow velocity, uj eFor total displacement (4 result of extraction step of unit e node jIn corresponding displacement value), Nj(ξ) is the same as in step 3
Definition, m0For mean unit linear mass, mt eFor the linear mass of unit e.
Step 6 calculates critical flow velocity V under every first-order modal using intrinsic frequency resultc
In formula, C is Connors coefficient, and ζ is damping ratio, DoFor exchange heat pipe outside diameter, f be intrinsic frequency (extraction step 4
Calculated result).
Step 7: calculate the velocity ratio of every first-order modal, i.e., equivalent flow velocity/critical flow velocity, if all velocity ratios less than 1,
Then elasticity of fluid unstability check passes through；It is more than or equal to 1 if there is some velocity ratio, then elasticity of fluid unstability is checked obstructed
It crosses, can be modified at this time by observing bending vibation mode picture to original design.
Example 1 (general issues): input condition refers to the example of heat exchanger standard GB/T 151-2014 appendix C .7
This algorithm calculated result and GB/T 151-2014 calculated result are compared as follows:
It is above-mentioned statistics indicate that, for conventional heat transfer pipe outside single pipe homogeneous media flow elasticity of fluid destabilization problems,
The present invention is consistent with the calculated result of GB/T 151-2014.
For complicated shape heat exchanger tube the elasticity of fluid destabilization problems of a variety of medium Non-Uniform Flows calculating, reference can be made to
Example 2.
Example 2: elasticity of fluid destabilization problems of the coiled pipe tube bank outside a variety of pipes under medium Non-Uniform Flow are calculated.
This algorithm calculated result is as follows:
In addition, the bending vibation mode picture of different modalities can be obtained by model analysis.It is greater than 1 mode for velocity ratio, it can be with
By observing bending vibation mode picture, judge the position that elasticity of fluid unstability occurs, to be then designed modification to weakness zone.
Claims (3)
1. a kind of flat tube beam fluid elastic instability assessment method, which comprises the following steps:
Step 1, the axis using heat exchanger tube, one-dimensional curve model heat exchanger tube being reduced in three-dimensional space；By one-dimensional curve
If model is split as main section, discrete lines segment model and master mould must be guaranteed geometrically almost the same；
Three dimensional cartesian coordinates system is defined in space, and list is carried out to node coordinate, cell node number, joint constraint, with
Continue after an action of the bowels and calculate, wherein discrete line segment is defined as unit, the endpoint of discrete line segment is defined as node, the node packet of unit
Include two endpoints and the line segment midpoint of discrete line segment, totally three:
Node coordinate list LOC table is shown as:
In formula, xi、yi、ziThe x-axis of respectively i-th node, y-axis, z-axis coordinate；
Cell node numbered list NOD is indicated are as follows:
In formula,The number of 1st, 2,3 node of respectively e-th unit, the 1st, 2,3 section of e-th of unit
Point is two endpoints and an intermediate node of e-th of unit；
Joint constraint list UROT is indicated are as follows:
In formula, uxi、uyi、uzi、rotxi、rotyi、rotzi6 freedom degrees of respectively i-th node constrain；
Step 2, according to Timoshenko beam theory, finite element theory and numerical integration method, computing unit mass matrix MeAnd
Element stiffness matrix Ke:
In formula, ngaussmAnd ngausskThe respectively Gauss integration point quantity of Mass matrix and Stiffness Matrix, ξiFor Gauss point position, wiFor
Weighting；je=le/ 2 are converted into the Jacobian determinant of unit local coordinate, l for canonical domaineIt is the length of unit line segment:
N (ξ) is the matrix function about canonical domain coordinate ξ:
In formula,For the shape function of cell node 1；
For the shape function of cell node 2；
N3=N3(ξ)=1- ξ2For the shape function of cell node 3；
Matrix C is MeSection correlation matrix:
In formula, ρ, A, J, Iy, IzRespectively cell density, sectional area, torsional moment of inertia, y to bending the moment of inertia, z to bending inertia
Square；
B (ξ) is the matrix function about canonical domain coordinate ξ:
In formula,
Matrix D is KeSection correlation matrix:
In formula, E, Ky、Gy、Kz、Gz、Kx, G be respectively that Young's modulus, y are repaired to Splice variant, y to modulus of shearing, z to shearing
Positive coefficient, z to modulus of shearing, torsion correction factor, modulus of shearing；
Matrix T is transformation matrix of the unit local coordinate to world coordinates:
In formula,
Step 3, the element mass matrix M that will be calculated in step 2eAnd element stiffness matrix KeAccording to the degree of freedom on a node basis number according to
It is secondary to add up into total quality matrix M and Bulk stiffness matrix K:
In formulaIt not sums symbol, but represents the assembling of total quality, stiffness matrix；
In assembling, by MeAnd KeIn all elements be added in M and K according to following matrix position No. corresponding relationship:
[(m-1) * 3+p-1, (n-1) * 3+q-1]
→ [(NOD [e, m] -1) * 3+p-1, (NOD [e, n] -1) * 3+q-1]
Wherein m, n (=1~3) are respectively cell node number, and p, q (=1~6) are respectively the degree of freedom on a node basis.
After being completed, according to table UROT, for restrained freedom degree, big number is set in the corresponding pivot of M and K；
Step 4, the characteristic value and feature vector that total quality matrix and Bulk stiffness matrix is calculated:
Above formula is generalized eigenvalue equation, wherein ω2It is characteristic value,It is feature vector, can be obtained after the completion of equation solution every
The circular frequency ω or intrinsic frequency of first-order modalAnd the motion vector of all nodes
Step 5, in the case of different from density along pipe flow velocity, using Connors formula calculate under every first-order modal etc.
Imitate flow velocity ve 2:
In formula, ρ0For average extratubal fluid density, ρeFor unit e extratubal fluid density, veFor the transverse flow of unit e extratubal fluid
Speed, uj eFor total displacement (4 result of extraction step of unit e node jIn corresponding displacement value), Nj(ξ) with the definition in step 3,
m0For mean unit linear mass, mt eFor the linear mass of unit e, wiAnd ξiWeighting and integral for 3 Gaussian quadratures
Point position (the same step 3) of value, jeFor Jacobian determinant (the same step 3) of value；
Step 6 calculates critical flow velocity V under every first-order modal using intrinsic frequency resultc
In formula, C is Connors coefficient, and ζ is damping ratio, DoFor the pipe outside diameter that exchanges heat, f is intrinsic frequency；
Step 7: calculating the velocity ratio of every first-order modal, i.e., equivalent flow velocity/critical flow velocity, if all velocity ratios are flowed less than 1
The check of body elasticity unstability passes through；Being more than or equal to 1 if there is some velocity ratio, then elasticity of fluid unstability check does not pass through, this
When can by observe bending vibation mode picture, to original design modify.
2. a kind of flat tube beam fluid elastic instability assessment method as described in claim 1, which is characterized in that in the step
In 2, the calculation formula of cell density ρ are as follows:
In formula, ρtFor heat exchanger tube density；miFor unit length tube fluid quality, riFor heat exchanger tube inside radius, ρi
For tube fluid density；moFor unit length extratubal fluid quality, roFor heat exchanger tube outer radius, ρoFor pipe
Outflow volume density, CMFor mass coefficient.
3. a kind of flat tube beam fluid elastic instability assessment method as claimed in claim 2, which is characterized in that in the step
In 2, Kz、KyCalculation formula are as follows:
In formula, υ is Poisson's 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 true CN109033514A (en) | 2018-12-18 |
| CN109033514B 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) |
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111680378A (en) * | 2020-07-17 | 2020-09-18 | 天华化工机械及自动化研究设计院有限公司 | ANSYS-based heat exchanger tube bundle modal analysis method in liquid filling state |
| CN115114872A (en) * | 2022-07-20 | 2022-09-27 | 中国核动力研究设计院 | Parameter identification method and system for predicting tube bundle fluid bomb instability |
| CN115238494A (en) * | 2022-07-21 | 2022-10-25 | 中国核动力研究设计院 | Method for identifying position of part of pipeline fluid bomb instability |
Citations (9)
| 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 |
| US20080208558A1 (en) * | 2007-02-27 | 2008-08-28 | Huayong Wang | System and method for simulating a multiprocessor system |
| 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 |
-
2018
- 2018-06-15 CN CN201810626314.0A patent/CN109033514B/en active Active
Patent Citations (9)
| 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 |
| US20080208558A1 (en) * | 2007-02-27 | 2008-08-28 | Huayong Wang | System and method for simulating a multiprocessor system |
| 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 |
|---|
| 吴一红等: "管阵流体弹性不稳定性分析", 《应用力学学报》 * |
| 张亚楠等: "两相横向流诱发管束弹性不稳定性的试验研究进展", 《轻工机械》 * |
| 张明贤，聂清德，侯曾炎，王献旗: "换热器管束流体弹性不稳定性的预测", 《石油化工设备》 * |
| 牛忠华等: "高压加热器管束振动机理简析", 《电站辅机》 * |
| 王莉华等: "三段式高压加热器管束流体弹性失稳临界流速的计算", 《力学季刊》 * |
| 聂清德，郭宝玉，丁学仁，金楠，陈旭: "换热器管束中的流体弹性不稳定性", 《力学学报》 * |
| 聂清德等: "管束的流体弹性不稳定性的研究", 《振动工程学报》 * |
| 薄涵亮，马昌文，佟允宪，姚梅生: "螺旋管束流体弹性不稳定", 《核科学与工程》 * |
| 韩同行等: "蒸汽发生器传热管流弹失稳计算", 《压力容器》 * |
| 黄平等: "影响管阵流体弹性不稳定性因素的研究", 《核动力工程》 * |
Cited By (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111680378A (en) * | 2020-07-17 | 2020-09-18 | 天华化工机械及自动化研究设计院有限公司 | ANSYS-based heat exchanger tube bundle modal analysis method in liquid filling state |
| CN111680378B (en) * | 2020-07-17 | 2022-09-16 | 天华化工机械及自动化研究设计院有限公司 | ANSYS-based heat exchanger tube bundle modal analysis method in liquid filling state |
| CN115114872A (en) * | 2022-07-20 | 2022-09-27 | 中国核动力研究设计院 | Parameter identification method and system for predicting tube bundle fluid bomb instability |
| CN115114872B (en) * | 2022-07-20 | 2023-09-26 | 中国核动力研究设计院 | Parameter identification method and system for predicting tube bundle fluid bullet instability |
| CN115238494A (en) * | 2022-07-21 | 2022-10-25 | 中国核动力研究设计院 | Method for identifying position of part of pipeline fluid bomb instability |
| CN115238494B (en) * | 2022-07-21 | 2023-10-20 | 中国核动力研究设计院 | Component position identification method for pipeline flow bullet instability |
Also Published As
| Publication number | Publication date |
|---|---|
| CN109033514B (en) | 2023-04-28 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Zhou et al. | Nonlinear vibration control of a cantilevered fluid-conveying pipe using the idea of nonlinear energy sink | |
| CN109033514A (en) | A kind of flat tube beam fluid elastic instability assessment method | |
| Shen et al. | Implicit WENO scheme and high order viscous formulas for compressible flows | |
| Nagaoka et al. | Accurate fluid-structure interaction computations using elements without mid-side nodes | |
| Reigstad | Numerical network models and entropy principles for isothermal junction flow | |
| Guo et al. | Uncertain frequency responses of clamp-pipeline systems using an interval-based method | |
| WO2011155539A1 (en) | Numerical analysis device, element generation program, and numerical analysis method | |
| CN105260806B (en) | A kind of fluid structurecoupling kinetic characteristics prediction technique of pipe-line system | |
| Yamashita et al. | Numerical convergence of finite element solutions of nonrational B-spline element and absolute nodal coordinate formulation | |
| Sonawane et al. | High resolution incompressible flow computations over unstructured mesh using SDWLS gradients | |
| Ji et al. | Numerical analysis of shell-side flow-induced vibration of elastic tube bundle in heat exchanger | |
| WO2024051724A1 (en) | Flow-induced vibration test apparatus and method, computer device, storage medium, and product | |
| Scholten et al. | An uncoupled method for fluid-structure interaction analysis with application to aerostructural design | |
| Xie et al. | A numerical simulation of VIV on a flexible circular cylinder | |
| Anderson et al. | High-Order Stabilized Finite Elements on Dynamic Meshes | |
| Cao et al. | Dynamic modeling and experimental verification of clamp–pipeline system with soft nonlinearity | |
| Righi et al. | ROM-based uncertainties quantification of flutter speed prediction of the BSCW wing | |
| Inci et al. | The effect and selection of solution sequence in co-simulation | |
| JP6747957B2 (en) | Method of analyzing damping structure | |
| Patalano et al. | A critical exposition of model order reduction techniques: application to a slewing flexible beam | |
| CN109726454B (en) | Fluid-solid coupling modeling method and device for pipeline system | |
| Moosavian et al. | Unified matrix frameworks for water hammer analysis in pipe networks | |
| Dou et al. | A novel retaining clip for vibration reduction of fluid-conveying pipes by piecewise constraints | |
| CN113050274A (en) | Triangular lattice phononic crystal band gap design method based on wavelet boundary element model | |
| Vourganti et al. | Effect of nonlinear cladding stiffness on the stability and Hopf bifurcation of a heat-exchanger tube subject to cross-flow |
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 |