Введение
В публикации [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
