Математическая модель вибрационного уровнемера с резонатором в виде консольной эллиптической трубки


    Введение


    В публикации [1] подробно рассмотрена реализация на Python метода измерения отношения частот с использованием фигур Лиссажу. В качестве примера были проанализированы формы колебаний консольной эллиптической трубки вибрационного уровнемера [2].



    Упруго закреплённая трубка эллиптического сечения с помощью систем возбуждения 5,6,7 совершает автоколебания в одной плоскости, а с помощью систем 8, 9, 10 в другой плоскости перпендикулярной первой. Трубка колеблется в двух взаимно перпендикулярных плоскостях с разными частотами близкими к собственным. Масса трубки зависит от уровня заполняющей её жидкости.

    С изменением массы меняются и частоты колебаний трубки, которые и являются выходными сигналами уровнемера. Частоты несут дополнительную информацию о мультипликативных и аддитивных дополнительных погрешностях, компенсируемых при обработке частот микропроцессором 11.

    Остался не решённым вопрос определения зависимости частот колебаний трубки от уровня заполняющей жидкости что и является предметом данной публикации.

    Постановка задачи


    Определить частоты изгибных колебаний трубки в двух взаимно перпендикулярных плоскостях методом Релея с использованием точного уравнения изгибной линии трубки из публикации [1].

    С использованием полученных соотношений для частот найти зависимости чувствительности от уровня и определить диапазоны пригодные для контроля уровня жидкости.

    Для реализации указанных задач средствами Python рассмотреть два метода решения символьный и символьно-численный. Сравнить указанные методы по производительности

    Применение метода Релея для определения частот колебаний консольной трубки эллиптического сечения


    Суть метода Релея состоит в определении максимальной кинетической Tmax и максимальной потенциальной Pmax энергий консервативной колебательной системы с последующим приравниванием энергий Tmax= Pmax.

    Из полученного уравнения можно найти частоту колебаний системы с учётом того, что все гармонические составляющие можно приравнять к максимальным значениям, например sin(t*w)**2=1 и cos(t*w)**2=1.

    Символьное решение на Python


    Применим следующее уравнение изгибной линии трубки:
    z=(sin(k*x)-sinh(k*x)+((cosh(k) -cos(k))/(sin(k)-sinh(k)))*(cos(k*x)-cosh(k *x))).
    Запишем следующий листинг программы для расчёта кинетической энергии пустой трубки длинной L по функции –def fp(L), на первой форме колебаний при k=1.875:

    def fp(L):
             x,A, L, w, t, k  =symbols('x A L w t k ')
             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)

    Для кинетической энергии жидкости, заполняющей трубку до уровня h и колеблющейся с трубкой как единое целое, вследствие малого сечения получим функцию def fh(h):

    def fh(h):
             x,A, h, w, t, k  =symbols(' x A h w t k ')
             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) будет иметь вид:

    def pl(L):
             x,A, J, w, t, k,L  =symbols('x A J w t k 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)

    Теперь можно составить уравнение Релея– Tmax= Pmax.Для этого умножим функции def fp(L), def fh(h), def pl(L), на коэффициенты (A**2*w**2*cos(t*w) **2) и (A**2*sin(t*w) **2), на которые их делили в листингах. Зачем сначала делить, а потом умножать поясню позже, когда будем рассматривать символьно-численный метод.

    Умножим функции кинетических энергий на массы единицы длины трубки и жидкости – m0 и m. Функцию потенциальной энергии умножим на жёсткость E*J. Решая уравнение Релея для круговой частоты – w, преобразовав её в циклическую –f=w/2*pi, получим:

    f=(0.5/pi)*((pl(L))*(E*J)/(m0*fp(L)+m*fh(h)))**0.5

    Имея функцию для частоты колебаний можно получить символьное решение.

    Листинг символьного решения
    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 A L w t k ')
             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 h w t k ')
             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 A J w t k 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 E J 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()


    Результат:


    Время: 165.868

    Из полученных графиков видно, что измерение или контроль можно осуществлять с уровня 0.75 длинны L трубки для водных растворов. Самая высокая чувствительность на конце трубки.

    Учитывая высокую стабильность частоты механического резонатора, порядка e-6, измерение можно осуществлять и на более широком участке от 0.5 L. А на участке от 0.75 L можно применить микропроцессорную обработку с компенсацией как аддитивной, так и мультипликативной составляющих погрешности.

    Вывод


    Символьное решение даёт хороший результат, но занимает много машинного времени.

    Символьно -численное решение на Python


    Можно получить символьное решение для функций def fp(L), def fh(h), def pl(L), поэтому мы и делили их на (A**2*w**2*cos(t*w) **2) и (A**2*sin(t*w) **2), чтобы сократить количество переменных до двух- L,h. А так же получить символьное решение для чувствительности уровнемера.

    Это можно сделать с помощью следующего листинга:

    Вспомогательный листинг для символьного отображения функций
    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 A L w t k ')
             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 h w t k ')
             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 A J w t k 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 E J 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))


    Время: 6.696

    Листинг не нагружен многократными вычислениями функций, поэтому время его выполнения 6.696 с. Тем более его нужно использовать всего один раз.

    Теперь остаётся скопировать результаты print(fp(L)), print(fh(h)), print(pl(L)), print(s) в соответствующие функции следующего листинга в котором удобно и быстро можно менять все исходные данные для расчёта. Результаты print не привожу, как и в предыдущих листингах, чтобы не загромождать текст.

    Листинг символьно – числового метода
    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()
    



    Результат:

    Время: 0.447

    Таким образом получили тот же результат что и символьным методом, но уже в 165.868/0.447=371 раз быстрее.

    Поверхность по которой колеблется трубка на первой форме её изгибных колебаний



    Листинг для построения поверхности
    #!/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()  



    Результат:



    Выводы:


    Определены частоты изгибных колебаний трубки в двух взаимно перпендикулярных плоскостях методом Релея с использованием точного уравнения изгибной линии трубки.

    С использованием полученных соотношений для частот найдены зависимости чувствительности от уровня и определены диапазоны пригодные для контроля или измерения уровня жидкости.

    Для реализации указанных задач средствами Python рассмотрены два метода решения символьный и символьно-численный. Символьно-численный метод выполняется в 370 раз быстрее чем символьный.

    Ссылки:

    1. От двух камертонов из опытов Лиссажу к одной эллиптической уровнемерной трубке с шагом в столетия и всё на Python.
    2. Вибрационный уровнемер.А.С.№777455
    Поделиться публикацией

    Комментарии 8

      –2

      "Имя функцию для частоты колебаний можно получить символьное решение".

        –1
        С теоретической точки зрения здорово. Но возникает ряд вопросов. Я так понял что стальная трубка закреплена снизу и трясут её сверху электромагнитами.
        1. Не понятно что мешает аппроксимировать график h(f1) и h(f2). Вместо вычисления кучи экспонент умноженных на 10^-15. Это позволит еще ускорить вычисления.
        2. Частоты 5 и 6 Гц как бы намекают на времена измерения n/(f2-f1) это десятки секунд на измерение. Как-то долго. Да и точность измерения вызывает вопросы?
        3. Как влияют внешние загрязнения и вибрации например насосов?
        4. Как много энергии необходимо подводить для поддержания автоколебаний.
        5. Может проще в трубке стоячие волны возбуждать, они гораздо проще связаны с уровнем.
        6. В целом метод выглядит довольно не практично. Как-то сложно, дорого еще и гудит.
          –1
          Отвечаю по порядку вопросов.
          1. Если аппроксимировать функции def fp(L), def fh(h), def pl(L) утрачивается основное преимущество метода Релея — точность определения частоты при точной упругой линии.
          Кроме того что для Вас «куча экспонент» для компьютера доли секунды и ускорять уже бессмысленно.
          2. В статье не предлагается вычитать частоты. Цифровая обработка осуществляется в реальном масштабе времени для каждой частоты отдельно.Точность измерения определяется нестабильностью частоты а она составляет e-6 или чтобы Вам было понятно 10^-6 (для механических резонаторов).
          3. В публикации я рассматривал схему устройства. В реальной конструкции используется механический фильтр в виде массивного основания и защитный кожух — как и для любого другого технологического оборудования.
          4. Резонатор закрепляется в узле изгибных колебаний при этом добротность достигает нескольких тысяч, а при резонансе энергия возбуждения равна энергии потерь на внутреннее трение в материале резонатора, а эти потери при малых амплитудах минимальны.
          5. Метод стоячих волн основан на разнице акустических импедансов сред поэтому имеет ограниченное применение и его реализация значительно сложнее чем вибрационно-частотного с механическим резонатором.
          6. Ведущие мировые фирмы выпускают вибрационные сигнализаторы уровня, поэтому утверждение о не практичности метода голословно.Интересно Ваше утверждение о гудении — резонатор колеблется на амплитудах мене 10^-3 м., и не генерирует звук.
            0
            1. Какая точность измерения высоты уровня? (Для датчиков обычно применяют микроконтроллеры, а не компьютеры)
            2. У вас в примере частоты 5 и 6 Гц. Сколько вам понадобится времени что бы измерить эту частоту с точностью 1%?
            3. Я не видел подобного оборудования. Но в других бывали резонансные частоты кожухов на 120Гц например.
            4. Смущает именно низкие частоты собственных колебаний. И геморой с калибровкой.
            5. Стоячие волны был как пример. Есть просто по отраженной волне и куча других (еще). Но все они используют высокие частоты или вообще не звуковое излучение.
            6. А они точно измеряют уровень или просто крайние точки?

            Какие преимущества есть у такого датчика перед другими? Что сподвигнет проектировщика к этому датчику?
              0
              Отвечаю.
              Для сигнализатора или уровнемера с узким диапазоном короткая трубка из кварцевого стекла и частота в кгц. Преимущество простота и и частотный выход. Возможен и погружной вариант в виде сплошной камертонной вилки — приведенный расчёт вполне пригоден и для такого варианта. О други мотивах поинтересуйтесь у производителей life-prog.ru/view_msinv.php?id=77.
                0
                Да но это тоже пороговый датчик и глубину он не измеряет. Он измеряет срыв колебаний. И частота 260Гц то есть он воет — жуть.
                  –1
                  В статье рассмотрен датчик с непрерывными автоколебаниями частота которых зависит от массы заполняющей жидкости а следовательно от её уровня. Резонатор не генерирует звук. При необходимости начальная частота может быть увеличена уменьшением длины и увеличением модуля упругости резонатор. Девиацию частоты увеличивают уменьшением плотности материала резонатора по отношению к плотности жидкости. При закреплении обеих концов трубки резонатор в форме эпилептической трубки заменяет уровне мерное стекло с отсчётом с уровня 0.75L.
                    +1
                    Из статьи не понятно, какая точность измерения и какие требования к датчику. А пример выглядит крайне не практичным. И из-за закрепления снизу точность измерения низких уровней (<0.25L) просто никакая. А как и сколько по времени вы собираетесь измерять разницу например между 5.21 Гц и 5.22 Гц?

        Только полноправные пользователи могут оставлять комментарии. Войдите, пожалуйста.

        Самое читаемое