Finite element modelling of Coriolis mass flowmeters with arbitrary pipe geometry and unsteady flow conditions

Flow Measurement and Instrumentation 37 (2014) 119–126

Contents lists available at ScienceDirect

Flow Measurement and Instrumentation journal homepage: www.elsevier.com/locate/flowmeasinst

Finite element modelling of Coriolis mass flowmeters with arbitrary pipe geometry and unsteady flow conditions J. Ruoff a,b,n, M. Hodapp a, H. Kück b a b

Festo AG & Co. KG, Ruiter Strasse 82, 73734 Esslingen, Germany Institute of Micro Assembly Technology, Hahn-Schickard-Gesellschaft, Allmandring 9B, 70569 Stuttgart, Germany

art ic l e i nf o

a b s t r a c t

Article history: Received 25 September 2013 Received in revised form 11 February 2014 Accepted 30 March 2014 Available online 15 April 2014

The behaviour of Coriolis mass flowmeters (CMFs) with arbitrary pipe geometry is investigated in this paper regarding the effects of unsteady flow, geometry and damping. Timoshenko's beam theory, the finite element method (FEM), Mathworks Matlab and Simulink are used for creating a suitable model. Various physical quantities such as Coriolis, centrifugal and gravitational forces, mass inertia, mass flow, unsteady flow forces, damping and material parameters are considered to increase model precision. Stiffening effects due to fluid deflection are partly neglected in the model. The use of exact finite elements and model reduction via modal transformation results in fast and accurate numerical computation which corresponds well with experimental data. As a consequence, this approach allows further optimization of pipe geometries, actuation positions and concepts, actuation und flow control systems and signal processing concepts without co-simulation. & 2014 Elsevier Ltd. All rights reserved.

Keywords: Coriolis flowmeter Finite element Dynamics Fluid structure interaction

1. Introduction Coriolis mass flowmeters are widely used in many industry segments because of their capability for direct mass flow measurement which is independent of pressure, density and temperature with low measurement error and high turn down ratios. The basic measurement principle is based on a flow pipe excited to vibrate in a fundamental mode, usually called drive mode, at the according resonance frequency. The amplitude of the motion is thereby very small, typically a fraction of a millimetre. Mass flow through the pipe then represents a relative movement to the pipe, on which the Coriolis force acts, resulting in a superimposed vibration of another mode, called a Coriolis mode. The amplitude of the Coriolis mode is hereby proportional to mass flow through the pipe and usually two or more magnitudes smaller as the drive mode and can be measured through position sensors as a phase shift of the drive mode. Commercially available metres exist with various pipe shapes [1], diameters ranging from 0.25 to 300 mm [2,3] for flow rates from 100 g/h up to 1500 t/h [2,4] optimized for different kinds of applications regarding sensitivity, operating pressure, pressure drop, usable installation orientation, installation space, piggability, power

n Corresponding author at: Festo AG & Co. KG, Ruiter Strasse 82, 73734 Esslingen, Germany. Tel.: þ49 711 347 57381. E-mail address: [email protected] (J. Ruoff).

http://dx.doi.org/10.1016/j.flowmeasinst.2014.03.010 0955-5986/& 2014 Elsevier Ltd. All rights reserved.

consumption and speed of response. Latest research activities have shown an increased effort in developing smaller micro-machined metres [3,5,6]. Because production and characterization of these metres and their prototypes are expensive and time consuming there are various different approaches to describe their behaviour. Each has different advantages, but also disadvantages like limitations in geometry and modelling precision [7–9], computation time [10] or the use of proprietary FEM codes [11,12], which make data manipulation and controller implementation difficult or co-simulation necessary. Implementing the whole system in Mathworks Matlab allows automatic computation, geometry, filter and controller optimization with the possibility of direct reuse of controller and filter configurations in prototypes. An analytical description of the behaviour of the flow pipe is usually restricted to simple straight geometries with limited accuracy through limitations in the formulation of physical effects [9,13]. Complex geometries usually are modelled using volume or beam elements to discretize the flow pipe. The state of technology modelling fluid conveying pipes is summarized by Ibrahim [14,15]. Through their geometry pipes it can precisely be modelled using beam elements, resulting in a smaller equation system compared to volume elements. Such approaches are shown by Belhadi et al. [12] and Wang et al. [8]. As shown there curved pipe sections can be discretized with a higher number of short straight beam elements resulting in negligible small errors in the modelling of the forces on the tube. Because these short elements are not

120

J. Ruoff et al. / Flow Measurement and Instrumentation 37 (2014) 119–126

necessarily slender enough for the validity of the Euler–Bernoulli beam theory [16], Timoshenko's beam theory [17] is usually used which accounts effects like shear deformation and rotary inertia. Previous studies used stationary flow conditions [7,8,18–20] or neglected several dynamic effects. A general formulation of the dynamics of fluid conveying pipes is described and derived by Laithier and Païdoussis [21] including the effects of flow changes, Coriolis and centrifugal forces as well as rotary inertia. Also differences regarding the derivation using Hamilton's or Newton's principle are discussed. A fully coupled dynamic fluid-structure interaction mechanism including axial, radial and transverse forces is described by Lee et al. [22,23]. In the same way a geometric non-linear formulation was derived by Lee and Chung [24] and Gorman et al. [25]. For the description of a Coriolis mass flowmeter Belhadi et al. [11,12] used a formulation for time-dependent flow by Païdoussis and Issid [26] of a single tube element. They further added the terms for dynamics to the finite element formulation for steady flow of Stack et al. [18]. However Païdoussis and Issid used a pressure of zero at the outlet of the pipe as boundary condition. Using this formulation for a finite element means that pressure at the start of every element needs to equal the inlet pressure of the discretized geometry and that the pressure at the end of every element needs to be zero. This formulation therefore only holds for a geometry discretized by one single finite element. Sadeghi and Karimi-Dona [27] used a formulation derived by Reddy and Wang [28] based on a Hamilton approach for closed systems. However, this approach is invalid through the momentum flow in and out of the considered system by the fluid as shown in [21,29], which is not considered in the boundary conditions. For the unbent tube this formulation proposes an axial acceleration and/or stretching caused by fluid acceleration in the absence of velocity depended pipe friction, where in reality the fluid acceleration is caused by a pressure difference. Some models of Coriolis mass flowmeters like Wang [30] or Belhadi [12] are based on formulations by Païdoussis [16,21]. These models are bound to the boundary conditions and simplifications made by Païdoussis, which are partly based on assumptions by Benjamin [31]. One of these assumptions is that at least one end of the beam is supported or clamped. Therefore axial movement of the pipe with exception of tension is suppressed. It is also assumed that the tension compression through lateral contraction is compensated by the strain through fluid friction. Therefore the whole problem can be considered frictionless [32]. This assumption implies that the upstream pipe end is fixed and that the force on the fixation does not need to be considered [33]. As a finite element needs to be able to describe forces in any direction this assumption is not suitable for a finite element formulation. 2. Strong form of governing equations of fluid conveying pipes In this paper a general finite element formulation based on the Timoshenko beam theory for the fully coupled fluid structure interaction proposed by Lee et al. [22,34] is developed considering dynamic flow conditions. Also suitable boundary conditions and simplifications are proposed. The forces on the pipe and the fluid can be identified in detail considering an infinitesimal element of a bended pipe and the conveyed fluid like shown in Fig. 1. The fluid element in Fig. 1(a) is hereby under the influence of ! the pressure p on the surface A over the length δs. The pressure ! gradient results in a time dependent tangential fluid velocity c ! and the opposing directed friction shear forces γ S on the fluid at ! the perimeter S . The pressure, fluid flow deflection and the mass

γ

γ

+

γ

γ

∂ δ ∂

Fig. 1. Forces acting on a fluid element (a) and the corresponding surrounding pipe element (b).

inertia of the fluid also result in forces directed normal to the fluid ! flow and the pipe which are summarized by N . The gravitational ! forces of the fluid are mf g . In Fig. 1(b) the tensile forces on the ! pipe are summarized by T and the gravitational of the pipe by ! mp g . The masses of the fluid mf and the pipe mp are the masses per length and can be calculated from the respective density and cross sectional area of the fluid Af or the pipe Ap . The damping ! forces on pipe movement are summarized in F Damp . In the following, the properties of the fluid and the pipe are marked with the indices f and p. With these quantities the fluid structure interaction can be described based on the derivation by Lee and Park [34] with the pressure and temperature forces described by Watkins and Anderson [35], Fyrileiv and Collberg [36] and Young and Budynas [37] with the following equations:
 νpRi _ 0 Þ ð1Þ E þ s0 þ  EαΔT Ap u″ þ γ S þ mf g y w0 ¼ mp ðu€ þ g x w0 w d

 
€  ςG I þ E þ s0 EαΔT þ vp 1 þ Ri J pp ϕ″ ¼ 0 ðςρf I pf þ ρp I pp Þϕ f pf 2ð1 þ vÞ d ð2Þ


EI p w″″ þ ðpAf w0 Þ0  s0 þ

νpRi d



EαΔT Ap w″

_ 0 þ c2 w″ þ cc0 w0 Þ þ mf ðc_ w0 þ 2cw € þ mg z þ mf g x w0 ¼ mw

ð3Þ

ðpA_ f Þ þ mf a2 ðc0  2νu_ 0 Þ ¼ 0

ð4Þ


 Af 1 D ν p1  ¼ ðEAp u0 Þ1 Ed Af 2Af

ð5Þ

! _ 0 þ u€ þ cu_ 0 Þ þ γ S ¼ 0 ðpAf Þ0 þ mf ðg z w0 þ g x þ c_ þ cc0 þw0 w

ð6Þ

where E is Young's modulus of the pipe, ν the Poisson ratio of the pipe, ρ the density, I p the polar moment of inertia, J pp the torsional moment of inertia of the pipe, Gf the shear modulus of the fluid, d the pipe wall thickness, Ri the inner pipe radius, Af and Ap the cross section areas, s0 the initial axial stress, α the coefficient of thermal extension of the pipe material, ΔT the temperature change since initial installation, p the pressure on the fluid, the axial displacement u, the transverse displacement w, the torsion ϕ

J. Ruoff et al. / Flow Measurement and Instrumentation 37 (2014) 119–126

and m the fluid and pipe mass combined. The velocity of sound of the fluid a is defined by Lee and Park [34] as a2 ¼

Ef Ed

ρf ðEd þ Ef DÞ

ð7Þ

with the bulk modulus of the fluid Ef and the fluid density ρf . The bulk modulus of the fluid is a material property characterizing its compressibility. The derivations with respect to time are marked as ð U Þ ¼ ∂=∂t, with respect to the location as ð0 Þ ¼ ∂=∂x and the substantial derivate as ð1Þ ¼ ∂=∂x þ c∂=∂t. Eq. (1) describes the dynamic axial behaviour of the pipe. The first term models the influence of internal pressure, temperature and initial tension as a modification of the elastic modulus. The following terms are for the fluid friction and gravitational influence on a bended pipe. On the right side are the terms for axial acceleration, axial acting gravitational and centrifugal forces. Eq. (2) describes the torsional behaviour of the pipe. The first term describes the mass moment of inertia in dependence of the polar moment of inertia of the pipe, the fluid and the rotation participation factor ς. The rotation participation factor describes how much the fluid follows the rotation of the tube and depends on the fluid viscosity, the pipe surface roughness, the speed of tube rotation and the tube diameter. It can be derived from the Navier–Stokes-equations based on the work of Drahm et al. [38]. For gases and low viscous liquids these terms only have very little influence. The second part of the equation describes the torsional stiffness. The influence of the pressure forces, like shown by Young and Budynas [37], and stresses which act through circumferential tension for isotropic materials on the torsion stiffness can be derived from Hook's Law in cylinder coordinates: 2 3 2 32 s 3 εrr rr 1 ν ν 6 7 7 6 ν 76 6 εϕϕ 7 6 sϕϕ 7 1  ν 0 6 7 6 7 6 7 6 7 6 εzz 7 szz 7 76 ν ν 1 6 7 16 7 76 6 7¼ 6 6 ð8Þ 76 s 7 6 εϕz 7 E 6 7 ν 1 þ ϕ z 6 7 6 7 7 6 76 6 ε 7 6 7 6 7 szr 7 1þν 0 6 zr 7 4 56 4 5 4 5 εr ϕ sr ϕ 1þν The influence of temperature and initial tension acting as axial stresses can be calculated as a modification of the shear modulus. The pressure force occurs as circumferential tension and can directly be added. The transverse dynamic of the pipe is described by Eq. (3) with the pressure and centrifugal forces of the fluid on curved paths [39,40], the Coriolis and the gravitational forces. The continuity in Eq. (4) of the fluid describes the correlation among pressure changes, fluid compression and axial strain. Eq. (5) contains the Poisson coupling between the radial and axial deformations of the pipe through Poisson coupling under the influence of pressure. The fluid momentum in Eq. (6) describes the relations between the pressure and the gravitational forces on the fluid and the resulting fluid acceleration and opposing friction forces. The resultant ! tangential shear force γ S caused by friction only equals the pressure drop if no other force like gravitational or centrifugal forces are acting on the fluid. The friction forces between the fluid and the pipe can be modelled for laminar and turbulent flow in a circular cross-section using the Darcy–Weisbach-equation [41] as !

γ S ¼ fD

! mf jcj c 2D

ð9Þ

with the Darcy friction factor f D . With the Eqs. (1)–(6) the variables pressure p, transverse displacement w, torsion ϕ , axial displacement u, the fluid velocity c and the radial displacement through the inner diameter D, which

121

is derived from the cross-sectional area Af , can be determined from given boundary conditions. If the radial deformation is not of interest, Eq. (5) can be neglected. The remaining system of five coupled equations can then be solved numerically using methods like the finite element method.

3. Simplified governing equations of fluid conveying pipes For most problems of interest, certain conditions allow the simplification of these equations resulting in faster computation with negligible error. In order to model Coriolis flowmeters for liquids and small pipe diameters the following assumptions and simplifications are made in this paper:

 The fluid is nearly incompressible. Therefore the wavelength is





  

 

much larger as the pipe length of the metre. The fluid velocity is therefore nearly constant over the metre length ðc0 ¼ 0Þ. The effects of a velocity profile are also neglected. The fluid pressure is assumed to be low. The resultant radial deformation in Eq. (5) has negligible influence on the dynamic behaviour of the metre. Therefore the cross section is constant over the whole geometry length ðA1 ¼ 0Þ. The arbitrary pipe geometry of the metre is modelled with the help of straight pipe elements which are treated as Timoshenko beams, including transverse shear deformation, rotary inertia and transverse acceleration. The radii of curvature of the modelled pipe geometry are much larger than the tube diameter. The pressure loss through deflection is therefore neglected. The pressure varies linear over a pipe element. The pipe diameter is small and the influence of gravitational forces on the fluid acceleration is compensated by capillary forces, fluid surface tension and viscosity. Therefore the gravitational forces on the fluid act in the same manner as the gravitational forces on the pipe. The deformations of the flow pipe of the Coriolis metre are small compared to the metre geometry. Therefore all nonlinear geometrical effects are neglected. The centrifugal forces on the fluid are neglected. As the in- and outlet of the pipe geometry are fixed, the centrifugal forces only lead to local pressure changes. Over the whole pipe length, these changes add up to zero. For this reason, there is no resulting global force and therefore no fluid acceleration.

It should be noted that modelling axial stresses as modification of the stiffness like done in this paper or by Lee and Park [34] and Wang et al. [8] are an approximation of effects like stress stiffening and axial tension through nonlinear geometrical behaviour. This formulation only modifies the stiffness and does not describe the real occurring deformations. The axial forces through temperature changes and internal pressure are usually causing a mixture of stresses and strains depending on the boundary conditions and the pipe geometry. Stress stiffening and spin softening could be described using a non-linear formulation or a non-linear preceding calculation for the determination of the deformations, the stresses and the resulting stiffness modifications and a proceeding calculation on the modified system. As the considered axial effects only have minor influence on the dynamic behaviour though their assumed small value range the modification of stiffness is an appropriate approximation of these effects. With these assumptions the equations of motion can be simplified to
 νpRi E þ s0 þ  EαΔT Ap u″ þ γ S ¼ mp u€ þ mg x ð10Þ d

122

J. Ruoff et al. / Flow Measurement and Instrumentation 37 (2014) 119–126



 
€  ςG I þ E þ s0  EαΔT þ vp 1 þ Ri I pp ϕ″ ¼ 0 ðςρf I f þ ρt I p Þϕ f pf 2ð1 þ vÞ d ð11Þ


 νpRi _ 0 þ c2 w″Þ  EαΔT Ap w″ þ mf ðc_ w0 þ2cw EI p w″″ þ ðpAf w0 Þ0  s0 þ d € þ mg z ¼ mw

ð12Þ

! p0 Af þ mf c_ þ γ S ¼ 0

ð13Þ

0

þ ðp0 Af þ mf c_ Þw0  ks Gp Ap ðw″  θ Þ ¼  mg z

ð14Þ


 νpRi  EαΔT I p θ″ ¼ 0 ðρf I f þ ρp I p Þθ€  ks Gp Ap ðw0  θÞ  E þ s0 þ d ð15Þ with the polar moments of inertia of the pipe I, the shear modulus Gp and the shear correction factor [42] ð7 þ 6vÞð1 þ sÞ2 þ ð20 þ12vÞs2

; s¼

di do

ð16Þ

for a round pipe with inner diameter di , outer diameter do and Poisson's ratiov. To solve the differential equations on a one dimensional discretized geometry the following cubic hermit shape functions for bending [43,44]: 2 3 ð1  3η2 þ 2η3 þ ψ ð1  ηÞÞ 6 7 Lðη  2η2 þ η3 þ ψ ðη  η2 Þ=2Þ 7 1 6 6 7 T ½N w  ¼ ð17Þ 6 7 7 ð3η2  2η3 þ ψηÞ 1þψ6 4 5 2 3 2 Lð  η þ η þ ψ ðη  ηÞ=2Þ quadratic shape functions for shear deformation 2 3 6ð  η þ η2 Þ=L 6 7 ð1  4η þ 3η2 þ ψ ð1  ηÞÞ 7 1 6 6 7 ½N θ T ¼ 6 7 7 1þψ6 6ðη  η2 Þ=L 4 5 ð  2η þ 3η2 þ ψηÞ and linear shape functions for torsion and axial tension " # 1η ½N L T ¼

η

ð18Þ

ð19Þ

with

ψ¼

12EI ks GAL2

; η¼

x L

ð20Þ

are used for consistent interpolation without shear locking [45]. The length of the element is herby described by L. Based on the Eqs. (10), (11), (14), (15) and the corresponding shape functions, the terms for the element matrices of the finite element formulation for bending can be derived as Z L M 44 ¼ ðmp þ mr ÞN Tw N w þ ðρf I f þ ρp I p ÞN Tθ N θ dx ð21Þ 0

0

Z K 44 ¼

L

2mf cðNTw N0w  N w N0T w Þdx



s0 þ

νpRi

d 0 þ ðp0 Af þmf c_ ÞðNTw N 0w  Nw N w'T^ Þ
 νpRi 0  EαΔT N 0T  E þ s0 þ θ Nθ d

0 0T þ ks GAðN θ N Tθ N 0w N Tθ  N θ N 0T w  N w N w Þdx

Z K 22 ¼

L




0

E þ s0 þ

νpRi d

 0 EαΔT Ap N 0T L N L dx

and torsion Z L ðςρf I pf þ ρt I pp ÞN TL N L dx M 22 ¼ 0

Z K 22 ¼

ð22Þ

  0  EαΔT Ap þ pAf þ mf c2 N 0T w Nw

0

The behaviour of fluid conveying pipes in a plane can be described using Timoshenko's beam theory. The governing equations, for the transverse deflection w at a specific point x at the pipe at time t can be derived from Eq. (12) for the bending and shear deformation as

  νpRi € þ 2mf cw _ 0 mw  EαΔT Ap þ pAf þ mf c2 w″ s0 þ d

6ð1 þ vÞð1 þ sÞ2

L

and analogue for axial tension Z L ðmp þ mr ÞN TL NL dx M 22 ¼

4. Finite element formulation

ks ¼

Z B44 ¼

L 0




ςGf Ipf þ

ð23Þ

ð24Þ ð25Þ

ð26Þ



  E þ s0  EαΔT R 0 þ vp 1 þ i J pp N 0T L N L dx 2ð1 þ vÞ d ð27Þ

with the torsional moment of inertia J pp of the pipe. Because the formulation for the transverse deflection in both transverse plans w and v is equal for symmetric cross sections, the terms in the elements matrices are identical too. Therefore the well known 12  12 element matrices for mass M e , damping Be and stiffness K e can be determinated from the terms . From these element matrices the dynamic equation of motion of the metre's pipe geometry can be assembled and written in standard matrix form for dynamic problems as 2

M

! d q dq þ BðcÞ þKðc; c_ ; c2 ; pÞq ¼ Fðt; c; pÞ þ F D ðc; psÞ dt dt 2

ð28Þ

where q is the vector of nodal displacements. The position, velocity and acceleration independent forces can be added to the model on the right side as shown by Friedman [43] via Z L " ½N  #T " Fðt; c; pÞ # w F¼ dx ð29Þ ½N θ  Mðt; c; pÞ 0 as consistent nodal forces F and nodal moments M. In addition to the Coriolis force, a speed dependent Rayleigh Damping is added to the velocity dependent damping matrix B to model material and air damping. More advanced approaches like squeeze film damping [46], damping through use of combination of viscoelastic material of Kelvin–Voigt type [47] and Navier–Stokes based air damping [48] also lead to acceptable results but are not as suitable as modal damping parameters fitted from measurement data for available metre geometries. The damping factors can easily be identified from experimental data, wisely from the Coriolis and the drive mode, and are usually very small through the high quality factor of the oscillating system. Damping has therefore only very little influence on the Coriolis amplitude or the phase shifts as long as actuation and Coriolis mode frequency are separated.

5. Element transitions As mentioned before curved pipe sections can be discretized with a higher number of short straight pipe elements with unsteady transitions like shown in Fig. 2. As the length of the elbow is short compared with the length of the straight elements the elbow parts between the straight

J. Ruoff et al. / Flow Measurement and Instrumentation 37 (2014) 119–126

123

dominant modes. For these reasons the modal reduction technique is used to accelerate computation. To allow direct evaluation in modal space the mode shapes are normalized to maximum translative amplitude. Pure torsion mode shapes are normalized to maximum angle of 2π.

7. Numerical computation

Fig. 2. Forces from flow deflection and pressure at unsteady element transitions.

elements are idealized as massless rigid parts. The pressure pi at distance x from the element start of each straight pipe element can then be approximated as L þx pi ðxÞ ¼ pI  ðpI pO Þ ai Lt

ð30Þ

with the inlet pressure pI , the outlet pressure pO , the distance from the inlet Lai and the total length of pipe geometry Lt . For small deflections the distance from the inlet Lai and the total length of pipe geometry Lt can be calculated from the distance of the undeformed geometry from the inlet plus the sum of all axial elongations due to strains between inlet and the desired position. The direction of the resulting force due to the fluid deflection ! e D is then ! ! ! e Di þ 1 ¼ e i  e i þ 1

ð31Þ ! with the unit direction vectors e i of the correspondent elements at this node. These vectors need to be calculated from the initial position of each element combined with the nodal displacements of the elements. With these quantities the influence of the change of direction and the influence of pressure at the discontinuous intersection of two tube elements can be described as ! ! F Di ¼ ðAf pi þ mf c2 Þ e Di

The computations were performed on a standard workstation using a single core at 4.5 GHz using below 1 GB of memory. To speed up computation matrix vectors were divided into constant elements with variable scalar factors. For time transient simulation the first 30 modes were considered resulting in a ratio of simulated time to simulation time greater one. For the excitation of the drive mode in time transient simulation a phased lock loop algorithm in combination with a proportional-integral controller was used.

8. Validation of the model The impetus behind this work is to achieve a reliable model of a Coriolis mass flowmeter. The described theoretical model is verified on the pipe geometry shown in Fig. 3. The actuation of the third eigenmode (EM) of this pipe geometry results in a fluid velocity dependent Coriolis force which excites the first eigenmode with the frequency of the actuation. Table 1 Comparison of different model reduction techniques ( þ þ very good,   inappropriate). Criteria

Modal

Balanced

Krylov

Preservation of stability Stationary accuracy Global error boundary Efficient computation Quality estimation

þþ þþ  0 þþ

þþ þþ þþ þ þ

 þ  þþ þ

ð32Þ

with pi the pressure at the inlet of the element i.

6. Model reduction The so established system of differential equations has six degrees of freedom of motion for every node of the pipe geometry. Usually some motions are more dominant than others and others are nearly direct coupled. For fast computation it is desirable to reduce the degrees of freedom. Different model reduction techniques like order reduction using Krylov subspace methods, Balanced Truncation, Proper Orthogonal Decomposition and modal order reduction after Litz have been investigated and are compared in Table 1. For a detailed explanation see [49]. Modal order reduction has the advantage of transparency through the physical background for mechanical systems. Also dominant terms can be identified after Litz [50] and stationary accuracy can be ensured after Guth [51]. Because Coriolis flowmeters usually are excited in one of their basic eigenmodes and higher order eigenmodes experience higher damping it is sufficient to only consider modes up to two magnitudes higher frequency than the exciting and the Coriolis mode. Alternatively the modal participation factors for each mode can be calculated in the unreduced system for a step to identify the

Eigenmode 1

Eigenmode 2

Eigenmode 3

Fig. 3. First three eigenmodes of the investigated pipe geometry.

Table 2 Comparison of computed eigenfrequencies with different software tools for model verification. Source

Element

EM 1

EM 2

EM 3

Measurement Mathworks Matlab 2012b Abaqus 6.12 Comsol multiphysics 4.1

– Beam elements Volume elements Volume elements

Ref.  0.45%  0.22% þ 0.56%

Ref.  0.62%  0.46% þ 0.61%

Ref. þ 0.25% þ 0.12%  0.19%

Normalized Drive Frequency

1.014 1.012 1.01 1.008 1.006 1.004

Measurement

1.002 1

Simulation

0.998 810

850

890

930

970

Density [kg /m3] Fig. 4. Third eigenfrequency as function of fluid density.

1010

1 0.8 0.6 0.4 0.2 Pipe 1

0 0

Pipe 2

Simulation

Theory

0.2 0.4 0.6 0.8 Internal fluid velocity / critical velocity

1

Fig. 5. Eigenfrequency as function of flow velocity.

Frequency(c) / Frequency(c=0)

To verify the basic behaviour through mass and stiffness distributions the eigenfrequencies of the system without mass flow are compared with simulated results from simulation programs in Table 2. The calculation in Mathworks Matlab is based on the beam formulation shown in this paper. The geometry model in Abaqus and the model in Comsol Multiphysics consist of volume elements using a geometric non-linear formulation. In these models the pipe wall thickness is modelled by three layers of elements to receive most realistic behaviour. This results in an over four magnitudes greater equation system. The reference measurements were done using a SIOS SP-S 120 laser vibrometer with a measurement uncertainty of the displacement amplitude of 0.3 nm and a uncertainty of the frequency measurement of 0.2 Hz in the investigated measurement range. The measured displacement amplitude and frequency were evaluated over 4096 samples resulting in a measurement uncertainty of 0.0048 nm and 0.003 Hz. Table 2 shows that all results are very close to the reference. The differences between the computational results are due to the different element formulations and computation algorithms. As the result of the simplified beam formulation lies between the more realistic non-linear geometric formulations, the simplification to use beam theory is valid. For further testing of the model the ability to measure the fluid density in the metre can be investigated. As mass, stiffness and internal pipe volume of the metre are constant and the stiffness of the fluid is usually negligible low, the eigenfrequencies of the metre shift in dependence of the fluid density. Ethanol water solutions with different densities have been used for the verification of the simulation. The results are shown in Fig. 4. In companion to the experiment, the simulation shows an offset in normalized drive frequency and also a deviation in the gradient. This offset is due to inaccurate damping parameters, differences in the material properties, production tolerances and the not ideally stiff fixation. The different gradients with an error of 4.57% are probably caused by a difference in the internal volume of the pipe to the assumed one of the metre. The influence of the fluid flow on the behaviour of the pipe is tested by an eigenfrequency calculation of a straight pipe, which is simply supported at its ends. The configuration of the simulation corresponds to the setup of the measurements performed by Dodds and Runyan [52]. The theoretical results for this setup are described by Blevins [53]. The variation in frequency of the bending vibrations as a function of fluid velocity is described by Housner's equation [54]. This equation also shows that at a critical fluid velocity a statically unstable condition exists, similar to the Euler buckling formula for beams. The results of the finite element beam formulation, the analytical solution of Housner's equation and measurement on two identical straight pipes are compared in Fig. 5. The frequency of the first bending vibration mode is normalized relative to the eigenfrequency with fluid present but without

Frequency(c) / Frequency(c=0)

J. Ruoff et al. / Flow Measurement and Instrumentation 37 (2014) 119–126

1 0.8 0.6 0.4

conventional inextensible theory (Chen) conventional ti l iinextensible t ibl th theory (Mi (Misra)) (Misra)) modified difi d iinextensible t ibl th theory (Mi M tl b Matlab

0.2 0 0

0.2

0.4

0.6

1

0.8

Internal fluid velocity / critical velocity Fig. 6. Eigenfrequency of the first eigenmode of a semi-circular flow pipe as function of flow velocity in comparison with different theories [16].

Coriolis amplitude [µm]

124

4 3 2

Simulation

1

Measurement

0 0

200

400 600 Mass flow [g/h]

800

1000

Fig. 7. Comparison of Coriolis amplitude for simulation and measurement.

flow. Fig. 5 shows a good correlation among theory, numerical approximation and experiment. As predicted by the theory, the system becomes unstable at the critical fluid velocity which results in a plastic deformation of the pipe. For a semi-circular flow pipe the eigenfrequency of the first eigenmode is shown in Fig. 6 in dependence of the flow velocity. It can be seen that the proposed formulation of the fluid velocity dependent eigenfrequencies for curved flow pipe geometries is in good agreement with other conventional inextensible theories. The disagreement to the modified inextensible or extensible 0 theories [16] is a result of the way the steady-state forces Π [16] on the element transitions are accounted for. They are modelled as fluid velocity dependent forces on the finite elements and have therefore no influence on the eigenfrequencies. The fluid velocity dependent forces on the element transitions are nevertheless formulated precisely and can describe fluid velocity dependent deformations of the flow pipe geometry. However, since the maximum flow velocity used in normal operation is usually very small compared to the critical fluid velocity this behaviour only leads to minor inaccuracies within the considered boundary conditions. For the investigated flow pipe geometry (Fig. 3) with a ratio of 0.1 for internal fluid velocity to critical velocity at nominal flow velocity, the frequency error is below 0.05%.

J. Ruoff et al. / Flow Measurement and Instrumentation 37 (2014) 119–126

125

2 1 0 -1 -2 -3

Mode 1

-4

Mode 4

Mode5

Mode 7

Mode 8

Time [s] 1.0

P sure [kP Press Pa]

70 60

0.8

50 0.6

40 30

0.4

20

0.2

10 0

Fluid velocity Time [s]

Pressure

F d vellocityy [m//s] Fluid

M Mode amp plitud de [µ µm]

3

0.0

Fig. 8. Mode amplitudes, fluid velocity and pressure difference over time for a flow change. The excitation mode is not shown in the upper diagram.

The complete model of the mass flowmeter with the geometry shown in Fig. 3 was then verified for different mass flows by the determination of the Coriolis amplitude. The dynamic behaviour mainly consists of a superposition of third eigenmode where the actuation takes place and the first eigenmode which is mass flow depended. The Coriolis amplitude can be measured at the centre at the bottom of the geometry where actuation amplitude is always zero (see Fig. 3, eigenmode 3). The measured and the simulated Coriolis amplitudes are shown in Fig. 7. The simulation and the measurement both show a linear dependence of the mass flow. The gradients differ by 2.24% which is the result of the previous differences regarding eigenfrequencies, mass distribution, stiffness and geometry errors.

9. Prediction of unsteady flow behaviour Experimental investigation of time dependent modal behaviour of the flow pipe implies high requirements on space and time measurement. Also the areas of the measurements points are very small and can move on curved paths, which makes contactless measurement difficult. Contactless measurement is necessary, as sensor elements attached to the pipe would influence the behaviour of the metre. Contactless excitation of the flow pipe can be realized by using a magnetic flow pipe and electromagnets or by inducing a current on the flow and (electro-) magnets. Alternatively the excitation can be realized by a small alternating force on the fixation on the flow pipe. If contactless measurement is not possible, attachments to the flow pipe like sensors and exciters can be modelled as an additional lumped mass by increasing the pipe density of the corresponding finite element if they are close the flow pipe. If they are further away from the flow pipe they need to be modelled by additional finite elements without fluid describing their additional mass inertia and position. Instead of these complex measurements the described and verified FEM model can be used to predict the dynamic behaviour of the metre to a mass flow change. The time transient simulation of unsteady flow for a metre with a pipe geometry in Fig. 3 is shown in Fig. 8. The modes four and eight are mainly influenced by gravity. The first mode is actuated by the Coriolis force. The modes five and seven are influenced by fluid flow forces and mode coupling. Modes which are not shown, with exception of the excitation mode, have an amplitude close to zero.

10. Conclusion In this study, a finite element beam model is developed for arbitrarily shaped three dimensional pipes for Coriolis mass flowmeters at internal unsteady flow. In addition to the usually presented equations for bending also terms for axial and torsion behaviour are included. Therefore a set of coupled fluid dynamic equations with consideration of boundary conditions and Timoshenko's beam theory are derived, in which the necessary physicals effects are considered. Next, the equations are adapted under certain boundary conditions and formulated using Timoshenko's beam theory. For different pipe geometries models using the presented finite element formulation, conventional FE models, analytical solutions and measurement data are compared. The results show the high accuracy of the model proposed in this paper. In order to investigate the dynamic characteristics of a pipe geometry carrying an unsteady liquid flow, a transient finite element analysis has been carried out. References [1] Baker RC. Flow Measurement Handbook: Industrial Designs, Operating Principles, Performance, and Applications. Cambridge University Press; 2005. [2] Bronkhorst Cori-Tech B.V.. Datasheet M12 Coriolis Mass Flowmeter for Liquids and Gases; 2012. [3] Clark C, Wang S, Cheesewright R. The perfomance characteristics of a micromachined Coriolis flowmeter: an evaluation by simulation. Flow Meas Instrum 2006;17:325–33. [4] Rheonik, RHM 160—the world's largest 12" Coriolis mass flowmeter, Rheonik RHM 160 Datasheet - v6, 2006. [5] Smith R, Sparks D, Riley D, Najafi N. A MEMS-based Coriolis mass flow sensor for industrial applications. IEEE Trans Ind Electron 2009;56(4):1066–71. [6] Frahnow, R. Massendurchfluss- und Dichtemessung mit einer resonanten Messzelle in Volumenmikromechanik [Dissertation]. Siemens und Technische Universität Chemnitz; 2007. [7] Samera G, Fana S. Modeling of Coriolis mass flowmeter of a general planeshape pipe. Flow Meas Instrum 2010;21(1):40–7. [8] Wang T, Baker RC, Hussain Y. An advanced numerical model for single straight tube Coriolis flowmeters. Trans ASME: J Fluids Eng 2006;128:1346–50. [9] Kutin J, Bajsic I. An analytical estimation of the Coriolis meter's characteristics based on modal superposition. Flow Meas Instrum 2002;12(5–6):345–51. [10] Mole N, Bobovnik G, Kutin J, Štok B, Bajsić I. An improved three-dimensional coupled fluid-structure model for Coriolis flowmeters. J Fluids Struct 2008;24 (4):559–75. [11] Wang S, Clark C, Cheesewright R. Virtual Coriolis flowmeter: a tool for simulation and design. Proc Inst Mech Eng Part C: J Mech Eng Sci 2006;220 (6):817–35. [12] Belhadi A, Cheesewright R, Clark C. The simulation of Coriolis meter response to pulsating flow using a general purpose F.E. Code. J Fluids Struct 2000;14:613–34. [13] Kazahaya M. A mathematical model and error analysis of Coriolis mass flowmeters. IEEE Trans Instrum Meas 2011;60(4):1163–74.

126

J. Ruoff et al. / Flow Measurement and Instrumentation 37 (2014) 119–126

[14] Ibrahim RA. Overview of mechanics of pipes conveying fluids—part I: fundamental studies. J Press Vessel Technol 2010;132(3):034001–32. [15] Ibrahim RA. Mechanics of pipes conveying fluids—part II: applications and fluidelastic problems. Trans ASME: J Press Vessel Technol 2011;133 (2):0240011–30. [16] Païdoussis MP. Fluid–Structure Interactions: Slender Structures and Axial Flow, 1. Academic Press, San Diego; 1998. [17] Timoshenko SP. On the correction for shear of the differential equation for transverse vibrations of prismatic bars. Philos Mag 1921;41(6):742–74. [18] Stack, CP, Garnett, RB and Pawlas, GE. A finite element for the vibration analysis of a fluid-conveying Timoshenko beam, AIAA/ASME/ASCE/AHS/ASC 34th structures, structural dynamics, and materials conference. vol. 1; 1993, p. 2120–2129. [19] Xia W, Wang L. Microfluid-induced vibration and stability of structures modeled as microscale pipes conveying fluid based on non-classical Timoshenko beam theory. Microfluid Nanofluid 2010;9(4):955–62. [20] Wang T, Baker RC. Manufacturing variation of the measuring tube in a Coriolis flowmeter. IEE Proc: Sci Meas Technol 2004;151(3):201–4. [21] Laithier BE, Païdoussis MP. The equations of motion of initially stressed Timoshenko tubular beams conveying fluid. J Sound Vib 1981;79(2):175–95. [22] Lee U, Pak CH, Hong SC. The dynamics of a piping system with internal unsteady flow. J Sound Vib 1995;180(2):297–311. [23] Lee U, Kim J. Dynamics of branched pipeline systems conveying internal unsteady flow. ASME J Vib Acoust 1999;121:114–22. [24] Lee SI, Chung J. New non-linear modelling for vibration analysis of a straight pipe conveying fluid. J Sound Vib 2002;254(2):313–25. [25] Gorman DG, Reese JM, Zhang YL. Vibration of a flexible pipe conveying viscous pulsating fluid flow. J Sound Vib 2000;230(2):379–92. [26] Païdoussis MP, Issid NT. Dynamic stability of pipes conveying fluid. J Sound Vib 1974;33(3):267–94. [27] Sadeghi MH, Karimi-Dona MH. Dynamic behavior of a fluid conveying pipe subjected to a moving sprung mass—an FEM-state space approach. Int J Press Vessels Pip 2011;88(4):123–31. [28] Reddy JN, Wang CM. Dynamics of Fluid-Conveying Beams, Centre of Offshore Research and Engineering. Singapore: National University of Singapore; 2004. [29] Benjamin TB. Dynamics of a system of articulated pipes conveying fluid. I. Theory. Proc R Soc Lond 1961;261(1307):457–86. [30] Wang T, Hussain Y. Pressure effects on Coriolis mass flowmeters. Flow Meas Instrum 2010;21(4):504–10. [31] Benjamin TB. Dynamics of a system of articulated pipes conveying fluid. I. Theory. Proc R Soc Lond 1961;261(1307):457–86. [32] Païdoussis MP, Li GX. Pipes conveying fluid: a model dynamical problem. J Fluids Struct 1993;7(2):137–204. [33] Öz HR. Non-linear vibrations and stability analysis of tensioned pipes conveying fluid with variable velocity. Int J Non-Linear Mech 2001;36:1031–9. [34] Lee U, Park J. Spectral element modelling and analysis of a pipeline conveying internal unsteady fluid. J Fluids Struct 2006;22(2):273–92.

[35] Watkins RK, Anderson LR. Structural Mechanics of Buried Pipes. CRC Press, Boca Raton; 2000. [36] Fyrileiv, O and Collberg, L. Influence of pressure in pipeline design–effective axial force. In: Proceedings of the 24th international conference on offshore mechanics and arctic engineering, (OMAE2005-67502), Halikidiki, Greece; 2005. [37] Young WC, Budynas RG. Roark's Formulas for Stress and Strain. 7th editionMcgraw-Hill Professional, New York; 2001. [38] Drahm W, Fuchs M, Bjö nnes H, Matt C, Wenger A, Sprich H-J. Coriolis-massflow-meter with direct viscosity measurement. Technisches Messen 2004, 148–153;71. [39] Catinaccio, A. Pipes under internal pressure and bending, ph-ep Tech-Note2009-004, Cern, Geneva, Switzerland. [40] Pak CH, Hong SC, Yun YS. On the vibrations of three-dimensional angled piping systems conveying fluid. KSME J 1991;5(2):86–92. [41] White FM. Fluid Mechanics. 4th editionMcGraw-Hill Higher Education, Boston; 1998. [42] Cowper GR. The shear coefficient in Timoshenko's beam theory. J Appl Mech 1966;33:335–40. [43] Friedman Z, Kosmatka JB. An impoved tow-node Timoshenko beam finite element. Comput Struct 1993;47(3):473–81. [44] Luo Y. An efficient 3D Timoshenko beam element with consistent shape functions. Adv Theory Appl Mech 2008;1(3):95–106. [45] Dhondt G. The Finite Element Method for Three-Dimensional Thermomechanical Applications. John Wiley & Sons, Chichester; 2004. [46] Bao M, Yang H. Squeeze film air damping in MEMS. Sens Actuators A: Phys 2007;136(1):3–27. [47] Meyers MA, Chawla KK. Mechanical Behavior of Materials. 2nd edition, Cambridge University Press, Cambridge; 2009. [48] Berli C, Cardona A. On the calculation of viscous damping of microbeam resonators in air. J Sound Vib 2009;327:249–53. [49] Antoulas Athanasios C. Approximation of Large-Scale Dynamical Systems (Advances in Design and Control). SIAM, Philadelphia; 2005. [50] Litz L. Reduktion der Ordnung linearer Zustandsraummodelle mittels modaler Verfahren. Hochschulsammlung Ingenieurwissenschaft Datenverarbeitung. Stuttgart: HochschulVerlag; 1979. [51] Guth R. Stationär genaue Ordnungsreduktion balancierter Zustandsraummodelle. Automatisierungstechnik 1991;39:286–90. [52] Dodds HL, Runyan HL. Effect of high-velocity fluid flow on the bending vibrations and static divergence of a simply supported pipe. 1965 (National Aeronautics and Space Administration Report NASA TN D-2870). [53] Blevins RD. Flow-Induced Vibration Second Edition. Krieger Publishing Company, Malabar, Florida; 2001. [54] Xi-de Z, Tao D, Wen Z, Wen-jun S. Correction for housner's equation of bending vibration of a pipe line containing flowing fluid. Appl Math Mech 1993;14(2):159–61.