Mathematical model of a vibrating level gauge with a resonator in the form of a cantilever elliptical tube
Introduction
In the publication [1], the implementation of the method of measuring the ratio of frequencies using Lissajous figures is considered in detail in Python. As an example, the vibration shapes of a cantilever elliptical tube of a vibrating level gauge were analyzed [2].
Elastically fixed tube of elliptical cross section with the help of excitation systems 5,6,7 performs self-oscillations in one plane, and with the help of systems 8, 9, 10 in another plane perpendicular to the first one. The tube oscillates in two mutually perpendicular planes with different frequencies close to its own. The mass of the tube depends on the level of its filling fluid.
')
With the change in mass, the oscillation frequencies of the tube also change, which are the output signals of the level gauge. Frequencies carry additional information about the multiplicative and additive additional errors that are compensated by the processing of frequencies by the microprocessor 11.
The question of determining the dependence of the oscillation frequency of the tube on the level of the filling fluid remained unresolved, which is the subject of this publication.
Formulation of the problem
Determine the frequency of bending oscillations of the tube in two mutually perpendicular planes by the Rayleigh method using the exact equation of the bending line of the tube from [1].
Using the obtained ratios for frequencies, find the dependences of sensitivity on the level and determine the ranges suitable for controlling the liquid level.
To implement these tasks using Python, consider two methods for solving symbolic and symbolic-numeric ones. Compare specified performance methods
Application of the Rayleigh method for determining the oscillation frequencies of a cantilever tube of elliptical section
The essence of the Rayleigh method is to determine the maximum kinetic Tmax and the maximum potential Pmax of the energies of a conservative oscillatory system, followed by equating the energies Tmax = Pmax.
From the resulting equation, it is possible to find the oscillation frequency of the system, taking into account that all harmonic components can be equated to maximum values, for example, sin (t * w) ** 2 = 1 and cos (t * w) ** 2 = 1.
Python character solution
Apply the following equation for the tube bend line:
z = (sin (k * x) -sinh (k * x) + ((cosh (k) -cos (k)) / (sin (k) -sinh (k))) * (cos (k * x) -cosh (k * x))).
We will write the following program listing for calculating the kinetic energy of an empty tube of length L according to the function –def fp (L), on the first form of oscillations at k = 1.875:
def fp(L): x,A, L, w, t, k =symbols('x AL wtk ') z=(sin(k*x)-sinh(k*x)+((cosh(k) -cos(k))/(sin(k)-sinh(k)))*(cos(k*x)-cosh(k *x))).subs({k:1.875}) dx=z*A*sin(w*t) dx2=(dx.diff(t))**2 return factor(integrate(dx2,(x,0,L)))/(A**2*w**2*cos(t*w)**2) For the kinetic energy of the fluid filling the tube to the level h and oscillating with the tube as a whole, due to the small cross section we obtain the function def fh (h):
def fh(h): x,A, h, w, t, k =symbols(' x A hwtk ') z=(sin(k*x)-sinh(k*x)+((cosh(k) -cos(k))/(sin(k)-sinh(k)))*(cos(k*x)-cosh(k *x))).subs({k:1.875}) dx=z*A*sin(w*t) dx2=(dx.diff(t))**2 return factor(integrate (dx2,(x,0,h)))/(A**2*w**2*cos(t*w)**2) For the potential energy of the tube, the function def pl (L) will have the form:
def pl(L): x,A, J, w, t, k,L =symbols('x AJ wtk L ') z=(sin(k*x)-sinh(k*x)+((cosh(k) -cos(k))/(sin(k)-sinh(k)))*(cos(k*x)-cosh(k *x))).subs({k:1.875}) dx=z*A*sin(w*t) dx2=(dx.diff(x))**2 return factor(integrate (dx2,(x,0,L)))/(A**2*sin(t*w)**2) Now we can make the Rayleigh equation– Tmax = Pmax. To do this, we multiply the functions def fp (L), def fh (h), def pl (L), by the coefficients (A ** 2 * w ** 2 * cos (t * w ) ** 2) and (A ** 2 * sin (t * w) ** 2), into which they were divided in listings. Why divide first, and then multiply, I will explain later when we consider the symbolic-numerical method.
Multiply the functions of kinetic energies by the masses of the unit length of the tube and the liquid — m0 and m. The function of potential energy is multiplied by the stiffness E * J. Solving the Rayleigh equation for circular frequency - w, transforming it into cyclic –f = w / 2 * pi, we get:
f = (0.5 / pi) * ((pl (L)) * (E * J) / (m0 * fp (L) + m * fh (h))) ** 0.5
Having a function for the oscillation frequency, you can get a symbolic solution.
Listing of symbolic solution
from sympy import * from numpy import arange,pi import matplotlib.pyplot as plt import time start = time.time() def fp(L): x,A, L, w, t, k =symbols('x AL wtk ') z=(sin(k*x)-sinh(k*x)+((cosh(k) -cos(k))/(sin(k)-sinh(k)))*(cos(k*x)-cosh(k *x))).subs({k:1.875}) dx=z*A*sin(w*t) dx2=(dx.diff(t))**2 return factor(integrate(dx2,(x,0,L)))/(A**2*w**2*cos(t*w)**2) def fh(h): x,A, h, w, t, k =symbols(' x A hwtk ') z=(sin(k*x)-sinh(k*x)+((cosh(k) -cos(k))/(sin(k)-sinh(k)))*(cos(k*x)-cosh(k *x))).subs({k:1.875}) dx=z*A*sin(w*t) dx2=(dx.diff(t))**2 return factor(integrate (dx2,(x,0,h)))/(A**2*w**2*cos(t*w)**2) def pl(L): x,A, J, w, t, k,L =symbols('x AJ wtk L ') z=(sin(k*x)-sinh(k*x)+((cosh(k) -cos(k))/(sin(k)-sinh(k)))*(cos(k*x)-cosh(k *x))).subs({k:1.875}) dx=z*A*sin(w*t) dx2=(dx.diff(x))**2 return factor(integrate (dx2,(x,0,L)))/(A**2*sin(t*w)**2) """ """ L,h,m,m0,E,J,q =symbols(' L,h m0 m EJ q ') k1=q.subs({q:0.16}) k2=pl(L).subs({L:1}) k3=E.subs({E:196e9})*J.subs({J:2.3e-08}) k4=m0.subs({m0:1.005}) k5=m.subs({m:0.98}) k6=fp(L).subs({L:1}) """ """ x=arange(0.0,1.0,0.01) f=k1*((k2+k3)/(k4*k6+k5*fh(h)))**0.5 y=[k1*((k2+k3)/(k4*k6+k5*fh(h).subs({h:w})))**0.5 for w in arange(0.0,1.0,0.01)] s=f.diff(h) y1=[s.subs({h:w}) for w in arange(0.0,1.0,0.01)] k3=E.subs({E:196e9})*J.subs({J:1.7e-08}) f1=k1*((k2+k3)/(k4*k6+k5*fh(h)))**0.5 y2=[k1*((k2+k3)/(k4*k6+k5*fh(h).subs({h:w})))**0.5 for w in arange(0.0,1.0,0.01)] s1=f1.diff(h) y3=[s1.subs({h:w}) for w in arange(0.0,1.0,0.01)] k7=fh(h).subs({h:0.7}) k3=E.subs({E:196e9})*J.subs({J:2.3e-08}) f1=k1*((k2+k3)/(k4*k6+k5*k7))**0.5 k3=E.subs({E:196e9})*J.subs({J:1.7e-08}) f2=k1*((k2+k3)/(k4*k6+k5*k7))**0.5 x1=[sin(f1*t) for t in arange(0.,2*pi,0.01)] y4=[sin(f2*t) for t in arange(0.,2*pi,0.01)] """ """ plt.subplot(221) plt.plot(x, y, label='f1-zox') plt.plot(x, y2,label='f2- zoy') plt.xlabel(' h') plt.ylabel(' f 1,f2') plt.legend(loc='best') plt.grid(True) plt.subplot(222) plt.plot(x, y1,label='s=df/dh zox') plt.plot(x, y3,label='s1=df1/dh - zoy') plt.ylabel('s') plt.xlabel(' h ') plt.legend(loc='best') plt.grid(True) plt.subplot(223) plt.plot(x1, y4,label='X,Y h=0.7') plt.ylabel(' Y') plt.xlabel('X') plt.grid(True) plt.legend(loc='best') stop = time.time() print (" :",round(stop-start,3)) plt.show() Result:
Time: 165.868
From the graphs obtained it can be seen that the measurement or control can be carried out from the level of 0.75 L tube length for aqueous solutions. Highest sensitivity at the end of the tube.
Taking into account the high frequency stability of the mechanical resonator, of the order of e-6, the measurement can be carried out on a wider section from 0.5 L. And on a section from 0.75 L, microprocessor processing can be applied with compensation for both additive and multiplicative error components.
Conclusion
A symbolic solution gives a good result, but takes a lot of machine time.
Python symbolic number solution
You can get a symbolic solution for the functions def fp (L), def fh (h), def pl (L), so we divided them into (A ** 2 * w ** 2 * cos (t * w) ** 2 ) and (A ** 2 * sin (t * w) ** 2) to reduce the number of variables to two- L, h. And also get a symbolic solution for the sensitivity of the level gauge.
This can be done using the following listing:
Auxiliary listing for symbolic display of functions
from sympy import * from numpy import arange,pi import matplotlib.pyplot as plt import time start = time.time() def fp(L): x,A, L, w, t, k =symbols('x AL wtk ') z=(sin(k*x)-sinh(k*x)+((cosh(k) -cos(k))/(sin(k)-sinh(k)))*(cos(k*x)-cosh(k *x))).subs({k:1.875}) dx=z*A*sin(w*t) dx2=(dx.diff(t))**2 return factor(integrate(dx2,(x,0,L)))/(A**2*w**2*cos(t*w)**2) def fh(h): x,A, h, w, t, k =symbols(' x A hwtk ') z=(sin(k*x)-sinh(k*x)+((cosh(k) -cos(k))/(sin(k)-sinh(k)))*(cos(k*x)-cosh(k *x))).subs({k:1.875}) dx=z*A*sin(w*t) dx2=(dx.diff(t))**2 return factor(integrate (dx2,(x,0,h)))/(A**2*w**2*cos(t*w)**2) def pl(L): x,A, J, w, t, k,L =symbols('x AJ wtk L ') z=(sin(k*x)-sinh(k*x)+((cosh(k) -cos(k))/(sin(k)-sinh(k)))*(cos(k*x)-cosh(k *x))).subs({k:1.875}) dx=z*A*sin(w*t) dx2=(dx.diff(x))**2 return factor(integrate (dx2,(x,0,L)))/(A**2*sin(t*w)**2) L,h,m,m0,E,J,q =symbols(' L,h m0 m EJ q ') f=q*((pl(L)+E*J)/(m0*fp(L)+m*fh(h)))**0.5 s=f.diff(h) print(fp(L)) print(fh(h)) print(pl(L)) print(s) stop = time.time() print (" :",round(stop-start,3)) Time: 6.696
Listing is not loaded with multiple calculations of functions, therefore its execution time is 6.696 s. Moreover, it should be used only once.
Now it remains to copy the results print (fp (L)), print (fh (h)), print (pl (L)), print (s) into the corresponding functions of the next listing in which you can quickly and easily change all the initial data for the calculation. I do not cite the results of print, as in previous listings, so as not to clutter the text.
Listing character - numeric method
from numpy import (pi,cos,cosh,sin,sinh,arange) import matplotlib.pyplot as plt import matplotlib as mpl mpl.rcParams['font.family'] = 'fantasy' mpl.rcParams['font.fantasy'] = 'Comic Sans MS, Arial' import time start = time.time() L=1# . d=1e-3# . e=0.8# .. a0=20e-3# . a=a0+d# . b0=a0*e# . b=b0+d# . E=196e9# /2. rm=7.9e3# /3. rg=1e3# /3. J=(pi/4)*(a**3*b-a0**3*b0)# 4. J1=(pi/4)*(a*b**3-a0*b0**3)# 4. m0=pi*(a*b-a0*b0)*rm# /3 m=pi*(a0*b0)*rg# /3 q=1/(2*pi)# def fp(L):# L return (1.83012983499975*L*sin(1.875*L)**2 + 1.83012983499975*L*cos(1.875*L)**2 - 0.830129834999752*L*sinh(1.875*L)**2 + 0.830129834999752*L*cosh(1.875*L)**2 - 0.869882798951083*sin(1.875*L)**2 + 0.442735911999868*sin(1.875*L)*cos(1.875*L) + 1.73976559790217*sin(1.875*L)*sinh(1.875*L) - 1.9521384906664*sin(1.875*L)*cosh(1.875*L) - 0.885471823999735*cos(1.875*L)*sinh(1.875*L) + 0.976069245333201*sinh(1.875*L)*cosh(1.875*L) - 0.869882798951083*cosh(1.875*L)**2 + 0.869882798951083) def fh(h):# h return (1.83012983499975*h*sin(1.875*h)**2 + 1.14978095631541e-15*h*sin(1.875*h)*sinh(1.875*h) + 1.83012983499975*h*cos(1.875*h)**2 - 1.07791964654569e-15*h*cos(1.875*h)*cosh(1.875*h) - 0.830129834999752*h*sinh(1.875*h)**2 + 0.830129834999752*h*cosh(1.875*h)**2 - 0.869882798951083*sin(1.875*h)**2 + 0.442735911999868*sin(1.875*h)*cos(1.875*h) + 1.73976559790217*sin(1.875*h)*sinh(1.875*h) - 1.9521384906664*sin(1.875*h)*cosh(1.875*h) - 0.885471823999735*cos(1.875*h)*sinh(1.875*h) + 0.976069245333201*sinh(1.875*h)*cosh(1.875*h) - 0.869882798951083*cosh(1.875*h)**2 + 0.869882798951083) def pl(L):# E*J/L**3 return (6.434050201171*L*sin(1.875*L)**2 + 1.48843545191282e-15*L*sin(1.875*L)*sinh(1.875*L) + 6.434050201171*L*cos(1.875*L)**2 - 2.79081647233653e-15*L*cos(1.875*L)*cosh(1.875*L) + 2.918425201171*L*sinh(1.875*L)**2 - 2.918425201171*L*cosh(1.875*L)**2 + 3.0581817150624*sin(1.875*L)**2 - 1.55649344062453*sin(1.875*L)*cos(1.875*L) + 3.11298688124907*sin(1.875*L)*cosh(1.875*L) - 6.86298688124907*cos(1.875*L)*sinh(1.875*L) + 6.1163634301248*cos(1.875*L)*cosh(1.875*L) - 3.0581817150624*sinh(1.875*L)**2 + 3.43149344062453*sinh(1.875*L)*cosh(1.875*L) - 6.1163634301248) def f(m0,q,E,J,h,m):# s=df/dh return -0.5*m0*q*((E*J + 6.434050201171*L*sin(1.875*L)**2 + 6.434050201171*L*cos(1.875*L)**2 + 2.918425201171*L*sinh(1.875*L)**2 - 2.918425201171*L*cosh(1.875*L)**2 + 3.0581817150624*sin(1.875*L)**2 - 1.55649344062453*sin(1.875*L)*cos(1.875*L) + 3.11298688124907*sin(1.875*L)*cosh(1.875*L) - 6.86298688124907*cos(1.875*L)*sinh(1.875*L) + 6.1163634301248*cos(1.875*L)*cosh(1.875*L) - 3.0581817150624*sinh(1.875*L)**2 + 3.43149344062453*sinh(1.875*L)*cosh(1.875*L) - 6.1163634301248)/(m*(1.83012983499975*L*sin(1.875*L)**2 + 1.83012983499975*L*cos(1.875*L)**2 - 0.830129834999752*L*sinh(1.875*L)**2 + 0.830129834999752*L*cosh(1.875*L)**2 - 0.869882798951083*sin(1.875*L)**2 + 0.442735911999868*sin(1.875*L)*cos(1.875*L) + 1.73976559790217*sin(1.875*L)*sinh(1.875*L) - 1.9521384906664*sin(1.875*L)*cosh(1.875*L) - 0.885471823999735*cos(1.875*L)*sinh(1.875*L) - 0.869882798951083*sinh(1.875*L)**2 + 0.976069245333201*sinh(1.875*L)*cosh(1.875*L)) + m0*(1.83012983499975*h*sin(1.875*h)**2 + 1.83012983499975*h*cos(1.875*h)**2 - 0.830129834999752*h*sinh(1.875*h)**2 + 0.830129834999752*h*cosh(1.875*h)**2 - 0.869882798951083*sin(1.875*h)**2 + 0.442735911999868*sin(1.875*h)*cos(1.875*h) + 1.73976559790217*sin(1.875*h)*sinh(1.875*h) - 1.9521384906664*sin(1.875*h)*cosh(1.875*h) - 0.885471823999735*cos(1.875*h)*sinh(1.875*h) - 0.869882798951083*sinh(1.875*h)**2 + 0.976069245333201*sinh(1.875*h)*cosh(1.875*h))))**0.5*(1.0*sin(1.875*h)**2 - 3.26206049606656*sin(1.875*h)*cos(1.875*h) - 2.0*sin(1.875*h)*sinh(1.875*h) + 3.26206049606656*sin(1.875*h)*cosh(1.875*h) + 2.6602596699995*cos(1.875*h)**2 + 3.26206049606656*cos(1.875*h)*sinh(1.875*h) - 5.32051933999901*cos(1.875*h)*cosh(1.875*h) + 1.0*sinh(1.875*h)**2 - 3.26206049606656*sinh(1.875*h)*cosh(1.875*h) + 2.6602596699995*cosh(1.875*h)**2)/(m*(1.83012983499975*L*sin(1.875*L)**2 + 1.83012983499975*L*cos(1.875*L)**2 - 0.830129834999752*L*sinh(1.875*L)**2 + 0.830129834999752*L*cosh(1.875*L)**2 - 0.869882798951083*sin(1.875*L)**2 + 0.442735911999868*sin(1.875*L)*cos(1.875*L) + 1.73976559790217*sin(1.875*L)*sinh(1.875*L) - 1.9521384906664*sin(1.875*L)*cosh(1.875*L) - 0.885471823999735*cos(1.875*L)*sinh(1.875*L) - 0.869882798951083*sinh(1.875*L)**2 + 0.976069245333201*sinh(1.875*L)*cosh(1.875*L)) + m0*(1.83012983499975*h*sin(1.875*h)**2 + 1.83012983499975*h*cos(1.875*h)**2 - 0.830129834999752*h*sinh(1.875*h)**2 + 0.830129834999752*h*cosh(1.875*h)**2 - 0.869882798951083*sin(1.875*h)**2 + 0.442735911999868*sin(1.875*h)*cos(1.875*h) + 1.73976559790217*sin(1.875*h)*sinh(1.875*h) - 1.9521384906664*sin(1.875*h)*cosh(1.875*h) - 0.885471823999735*cos(1.875*h)*sinh(1.875*h) - 0.869882798951083*sinh(1.875*h)**2 + 0.976069245333201*sinh(1.875*h)*cosh(1.875*h))) """ """ x=arange(0.0,1.0,0.01) y1=[q*((pl(L)+E*J)/(m0*fp(L)+m*fh(h)))**0.5 for h in arange(0.0,1.0,0.01)] # ZOX J=(pi/4)*(a*b**3-a0*b0**3)# ZOY y2=[q*((pl(L)+E*J)/(m0*fp(L)+m*fh(h)))**0.5 for h in arange(0.0,1.0,0.01)]# ZOX J=(pi/4)*(a**3*b-a0**3*b0)# ZOX y3=[f(m0,q,E,J,h,m) for h in arange(0.0,1.0,0.01)]# ZOX J=(pi/4)*(a*b**3-a0*b0**3)# ZOY y4=[f(m0,q,E,J,h,m) for h in arange(0.0,1.0,0.01)]# ZOX J=(pi/4)*(a**3*b-a0**3*b0)# ZOX h=0.7 f1=q*((pl(L)+E*J)/(m0*fp(L)+m*fh(h)))**0.5 # h=0.7 ZOX J=(pi/4)*(a*b**3-a0*b0**3)# ZOY f2=q*((pl(L)+E*J)/(m0*fp(L)+m*fh(h)))**0.5 # h=0.7 ZOY x1=[sin(f1*t) for t in arange(0.,2*pi,0.01)]# y5=[sin(f2*t) for t in arange(0.,2*pi,0.01)]# Y """ """ plt.subplot(221) plt.plot(x, y1, label='f-ZOX') plt.plot(x, y2,label='f1-ZOY') plt.xlabel(' h') plt.ylabel(' f 1,f2') plt.legend(loc='best') plt.grid(True) plt.subplot(222) plt.plot(x, y3,label='s1=df1/dh - ZOX') plt.plot(x, y4,label='s2=df2/dh - ZOY') plt.ylabel('s') plt.xlabel(' h ') plt.legend(loc='best') plt.grid(True) plt.subplot(223) plt.plot(x1, y5,label='X,Y h=0.7') plt.ylabel(' Y') plt.xlabel('X') plt.grid(True) plt.legend(loc='best') stop = time.time() print (" :",round(stop-start,3)) plt.show() Result:
Time: 0.447
Thus, we obtained the same result as the symbolic method, but already 165.868 / 0.447 = 371 times faster.
The surface on which the tube oscillates on the first form of its bending vibrations
Listing to build a surface
#!/usr/bin/env python #coding=utf8 import pylab import numpy from numpy import (zeros,arange,ones,pi,sin,cos,sinh,cosh) from matplotlib.colors import LinearSegmentedColormap from matplotlib import cm from mpl_toolkits.mplot3d import Axes3D def trubka(): a=10;b=10;k=q[0] v = arange(0, 2.05*pi, 0.05*pi) u= zeros([len(v),1]) for i in arange(0,len(v)): u[i,0]=[sin(k*w)-sinh(k*w)+((cosh(k) -cos(k))/(sin(k)-sinh(k)))*(cos(k*w)-cosh(k *w) )for w in arange(0, 2.05*pi, 0.05*pi)][i] x=a*u*cos(v) y=b*u*sin(v) z=u*ones(len(v)) return x,y,z x,y,z=trubka() fig = pylab.figure() axes = Axes3D(fig) axes.plot_surface(x, y, z, rstride=4, cstride=4, cmap = cm.jet) pylab.show() Result:
Findings:
The frequencies of bending tube oscillations in two mutually perpendicular planes are determined by the Rayleigh method using the exact equation of the tube's bending line.
Using the obtained relations for frequencies, the dependences of sensitivity on the level were found and ranges suitable for monitoring or measuring the level of the liquid were determined.
To implement these tasks using Python, two methods for solving symbolic and symbolic-numerical methods are considered. The character-numeric method runs 370 times faster than the character method.
References:
1. From two tuning forks from Lissajous experiments to a single elliptical level gauge tube with a step of centuries and all in Python.
2. Vibration level gauge.A.C.№777455
Source: https://habr.com/ru/post/340982/
All Articles