От двух камертонов из опытов Лиссажу к одной эллиптической уровнемерной трубке с шагом в столетия и всё на Python




    Картинки из сети, качество желает лучшего, но они достаточно точно отражают суть опыта по визуализации фигур. Зри в корень – основа мудрости поколений.

    Немного истории


    Ещё в школе на уроках физики я вглядывался в осциллограф, на экране которого, сменяя друг друга, появлялись разные фигуры: сначала простые – линия, парабола, круг, эллипс, потом фигуры становились всё более насыщенные непрерывными волнообразными линиями, напоминающие мне кружева. Автором этого кружевного дива был Жюль Антуан Лиссажу французский физик, член — корреспондент Парижской АН (1879) [1]. Сами фигуры — это замкнутые траектории, прочерчиваемые точкой, совершающей одновременно два гармонических колебания в двух взаимно перпендикулярных направлениях [2]. Думаю, что в те далёкие от современности годы основной заслугой Жюля, кроме конечно накопленных опытом знаний математики и физики, была простая механическая визуализация этих фигур подручными средствами. Захотелось конструировать подобно Жулю максимально просто и наглядно, реализовать его идеи применительно к современной задаче линейных измерений. Но сделать это путём математического моделирования с графической визуализацией его результатов на Python. Но сначала рассмотрим классический вариант [3] построения фигур.

    Какими должны быть фигуры Лиссажу


    Для этого воспользуемся системой уравнений, описывающих фигуры:



    x(t), y(t) в общем случае зависящие от времени гармонические колебания вдоль взаимно перпендикулярных плоскостей, частоты b, a и начальная фаза d. Для анализа фигур в вычислениях принимают постоянным модуль разности частот |b — a| = 1. Будем рассматривать отношение круговых частот b / a и начальную фазу d. Имеем для линии A = B d = 0, окружности , и параболы . Основные отношения частот, удовлетворяющие условию, занесём во вложенный список m=[[0],[2,2],[2,1],[1,2],[3,2],[3,4],[5,4],[5,6],[9,8]].

    Код для построения графиков каждой из фигур на отдельных графиках
    #!/usr/bin/env python
    #coding=utf8
    import numpy as np
    from numpy import sin,pi
    import matplotlib.pyplot as plt
    m=[[0],[2,2],[2,1],[1,2],[3,2],[3,4],[5,4],[5,6],[9,8]]# отношение круговых частот  
    for i in m:          
                    if i[0]==0: 
                               a=1
                               x=[sin(a*t) for t in np.arange(0.,2*pi,0.01)]
                               y=[sin(a*t) for t in np.arange(0.,2*pi,0.01)]
                               plt.plot(x, y, 'r')# график для линии
                               plt.grid(True)
                               plt.show()
                    else:
                                    a=i[0]
                                    b=i[1]
                                     d=0.5*pi
                                    x=[sin(a*t+d) for t in np.arange(0.,2*pi,0.01)]
                                    y=[sin(b*t) for t in np.arange(0.,2*pi,0.01)]
                                    plt.plot(x, y, 'r') # график для различных отношений  a/b
                                    #круговых частот
                                    plt.grid (True)
                                    plt.show()
    

    Результат не привожу, отдельные фигуры не впечатляют. Хочу коллаж из «кружев».

    Код программы для построения на одной форме графиков для четырёх фигур при m= [3,4], [5,4],[5,6],[9,8]]
    #!/usr/bin/env python
    #coding=utf8
    import numpy as np
    from numpy import sin,pi
    import matplotlib.pyplot as plt
    m=[[3,4],[5,4],[5,6],[9,8]] # отношение круговых частот
    plt.figure(1)
    for i in m:         
             a=i[0]
             b=i[1]
             d=0.5*pi
             x=[sin(a*t+d) for t in np.arange(0.,2*pi,0.01)]
             y=[sin(b*t) for t in np.arange(0.,2*pi,0.01)]        
             if m.index(i)==0:
                      plt.subplot(221)
                      plt.plot(x, y, 'k') # график для различных отношений  a/b круговых частот                                
                      plt.grid(True)
             elif m.index(i)==1:
                      plt.subplot(222)
                      plt.plot(x, y, 'g')               
                      plt.grid(True)
             elif m.index(i)==2:
                      plt.subplot(223)
                      plt.plot(x, y, 'b') 
                      plt.grid(True)
             else:
                      plt.subplot(224)
                      plt.plot(x, y, 'r')                  
                      plt.grid(True )    
    plt.show()
    


    И вот они «кружева».


    Что нельзя отнести к фигурам Лиссажу по определению о их замкнутости


    Зачем нам |b — a| = 1, “за флажки!” попробуем например так m=[[1,3],[1,5],[1,7],[1,9]]



    На втором графике при m=0,2 получена незамкнутая траектория, которая по определению не является фигурой Лbссажу.

    В поисках механических аналогов


    Поищем аналогии фигур в измерительной технике и вот вибрационный уровнемер с резонатором в виде эллиптической трубки [4].

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

    Условия адекватного моделирования



    Для более-менее корректной привязки фигур Лиссажу к работе упомянутого уровнемера, следует учесть следующие обстоятельства. Во-первых, закреплённая одним концом трубка эллиптического сечения — это колебательная система с распределёнными параметрами, что сильно усложняет анализ её колебаний. Во-вторых, отношение частот колебаний трубки не может изменяться произвольно, оно зависит от эллипсности сечения и допустимых зазоров в системе возбуждения колебаний. Для отношения частот можно получить простое соотношение.



    К чему принадлежат переменные, a, b, a0, b0 ясно из рисунка и кроме того формула для циклической частоты осциллятора известна из школьного курса физики. Для «реализации на Python в последнее отношение введём толщину стенки и показатель эллипсности внутреннего сечения трубки, тогда вместо четырёх переменных получим три.



    Код программы для определения. допустимого изменения отношения частот
    #!/usr/bin/env python
    #coding=utf8
    import numpy as np
    from numpy import sqrt
    import matplotlib.pyplot as plt
    import matplotlib as mpl
    mpl.rcParams['font.family'] = 'fantasy'
    mpl.rcParams['font.fantasy'] = 'Comic Sans MS, Arial'
    d=0.5
    a=9
    x=[w for w in  np.linspace(0.8,0.95,15)]
    y=[sqrt(x**-2*((1+d/a)**3*(1+d/(x*a))-1)/((1+d/a)*(1+d/(x*a))**3-1)) for x in  np.linspace(0.8,0.95,15)]
    plt.plot(x, y, 'r', label='Толщина стенки трубки в мм. --  %s' %str(d))
    d=0.7
    y=[sqrt(x**-2*((1+d/a)**3*(1+d/(x*a))-1)/((1+d/a)*(1+d/(x*a))**3-1)) for x in  np.linspace(0.8,0.95,15)]
    plt.plot(x, y, 'b',label='Толщина стенки трубки в мм.--  %s' %str(d))
    d=1.0
    y=[sqrt(x**-2*((1+d/a)**3*(1+d/(x*a))-1)/((1+d/a)*(1+d/(x*a))**3-1)) for x in  np.linspace(0.8,0.95,15)]
    plt.plot(x, y, 'g', label='Толщина стенки трубки в мм.--  %s' %str(d))
    plt.ylabel('Отношение частот колебаний эллиптической трубки')
    plt.xlabel('Отношение длин малой и большой полуосей')
    plt.title('Определение допустимого диапазона для отношения частот')
    plt.legend(loc='best')
    plt.grid(True)
    plt.show()
    


    В результате работы программы получим график.



    График построен для малой внутренней полуоси в 9 мм. Для конструктивно допустимого отношения малой к большой полуоси сечения в диапазоне от 0.8 до 0.95. Это основной фактор влияния на отношение частот, которое изменяется от 1.18 до 1.04. Толщина стенки влияет незначительно. Теперь у нас есть диапазон отношений и ним можно воспользоваться для дальнейшего моделирования.

    Формы колебаний вертикальной оси трубки


    Что касается распределённых механических параметров консольной трубки, то они при помощи равенства собственных частот и импеданса могут быть приведены к сосредоточенной массе жёсткости и демпфированию. Кроме того, для определения форм изгибных колебаний консольной трубки можно получить выражение для распределённых параметров. Уравнение для форм – балочные функции имеет вид:


    где — корни уравнения:


    Следует отметить что, не смотря на большое количество публикаций о формах и частотах колебаний консольного стержня, балки или трубки уравнения (4) нигде не приводяться, только рисунки без координат. Поэтому уравнение (4), я вывел через условия на концах и балочные функции, проверил по корням (5) и расположению узлов. Однако это тривиальное уравнение, о котором просто забыли.

    Код программы для численного определения корней уравнения 1.1 и построения трёх форм изгибных колебаний оси трубки
    1.1 —
    #!/usr/bin/env python
    #coding=utf8
    from scipy.optimize import *
    import numpy as np
    from numpy import pi,cos,cosh,sin,sinh
    import matplotlib.pyplot as plt
    import matplotlib as mpl
    mpl.rcParams['font.family'] = 'fantasy'
    mpl.rcParams['font.fantasy'] = 'Comic Sans MS, Arial'
    d=[]
    for i in range(0,4):
             x=brentq(lambda x:cosh(x)*cos(x)+1,0+pi*i,pi+pi*i)
             p=round(x,3)
             if p not in d:
                      d.append(p)
    x=[w for w in np.linspace(0,1,100)]
    k=d[0]
    z=[sin(k*x)-sinh(k*x)+((cosh(k) -cos(k))/(sin(k)-sinh(k)))*(cos(k*x)-cosh(k *x) )for x in np.linspace(0,1,100)]
    plt.plot(z, x, 'g', label='Первая форма для корня -  %s' %str(k))    
    k=d[1]
    z=[sin(k*x)-sinh(k*x)+((cosh(k) -cos(k))/(sin(k)-sinh(k)))*(cos(k*x)-cosh(k *x)) for x in np.linspace(0,1,100)]
    plt.plot(z, x, 'b', label='Вторая форма для корня -  %s' %str(k))   
    k=d[2]
    z=[sin(k*x)-sinh(k*x)+((cosh(k) -cos(k))/(sin(k)-sinh(k)))*(cos(k*x)-cosh(k *x)) for x in np.linspace(0,1,100)]
    plt.plot(z, x, 'r', label='Третья форма для корня -  %s' %str(k))
    plt.title('Первые три формы изгибных колебаний осевой линии трубки')
    plt.xlabel(' Координата вдоль оси OX ')
    plt.ylabel(' Координата положения осевой линии трубки вдоль оси OZ ')
    plt.legend(loc='best')
    plt.grid(True)
    plt.show()
    

    В результате работы программы получим график построенный с учётом вертикального положения трубки.



    На графике координата осевой линии приведена к длине трубки, а амплитуда нормирована. Положение узлов колебаний трубки относительно места её крепления в точности соответствует теории колебаний.

    По каким траекториям движется конец трубки


    Последнее препятствие — сложность получения осмысленного численного решения дифференциальных уравнений колебаний, при условии варьирования несколькими параметрами одновременно. Тут на помощь пришли две мои статьи о колебательном звене на Python [5,6], в которых приведена методика получения точных символьных решений дифференциальных уравнений.

    Запишем два условно независимых уравнения для колебаний трубки в плоскости OX и OY с разными частотами a и b отношение между которыми выбрано из ранее установленного диапазона. Остальные параметры выбраны во правильной взаимосвязи, но произвольно для лучшей демонстрации результата.



    Здесь введены следующие обозначения (для упрощения без индексов).

    ─ приведенная амплитуда силы, ─ коэффициент затухания, ─ собственная частота колебаний системы, m ─ сосредоточенная масса одинаковая для обоих уравнений, ─ сосредоточенные коэффициенты демпфирования, разные из-за разных амплитуд, а следовательно разных зазорах в системах возбуждения колебаний, ─ разные жёсткости из-за эллиптичности сечения трубки.

    Код программы для решения каждого дифференциального уравнения системы (6), с последующем сложением для получения траектории движения конца трубки.
    import numpy as np
    from sympy import *
    from IPython.display import *
    import matplotlib.pyplot as plt
    import matplotlib as mpl
    mpl.rcParams['font.family'] = 'fantasy'
    mpl.rcParams['font.fantasy'] = 'Comic Sans MS, Arial'
    def solution(w,v,i,n1,n2,B,f,N):
             t=Symbol('t')
             var('t C1 C2')
             u = Function("u")(t)   
             de = Eq(u.diff(t, t) +2*B*u.diff(t) +w**2* u, f*sin(w*t+v))        
             des = dsolve(de,u)
             eq1=des.rhs.subs(t,0)
             eq2=des.rhs.diff(t).subs(t,0)     
             seq=solve([eq1,eq2],C1,C2)
             rez=des.rhs.subs([(C1,seq[C1]),(C2,seq[C2])])
             g= lambdify(t, rez, "numpy")         
             t= np.linspace(n1,n2,N)
             plt.figure(1)         
             if i==1:
                      plt.subplot(221)
                      plt.plot(t,g(t),color='b', linewidth=3,label='x=%s*sin(%s*t+%s)' %(str(f),str(w),str(v)))
                      plt.legend(loc='best')
                      plt.grid(True)
             else:                 
                      plt.subplot(222)
                      plt.plot(t,g(t),color='g', linewidth=3,label='y=%s*sin(%s*t+%s)' %(str(f),str(w),str(v)))
                      plt.legend(loc='best')
                      plt.plot(t,g(t),color='r', linewidth=3)
                      plt.grid(True)                 
             return g(t)        
    N=1000#Число точек оцифровки временного интервала
    B=0.2#Установка демпфирования 
    f=1#Установка амплитуды
    n1=0#Нижняя граница временной развертки
    n2=20#Верхняя граница временной развёртки
    w1=5.0#Частота колебаний трубки вдоль оси ОХ
    w2=10.0#Частота колебаний трубки вдоль оси ОУ
    v1=0#Начальная фаза при колебании вдоль оси ОХ
    v2=0#Начальная фаза при колебании вдоль оси ОУ
    g1=solution(w1,v1,1,n1,n2,B,f,N)
    g2=solution(w2,v2,2,n1,n2,B,f,N)
    plt.subplot(223)
    plt.plot(g1,g2,color='b', linewidth=3,label='w1/w2=%s'%str(w1/w2))
    plt.legend(loc='best')        
    plt.grid(True)
    plt.subplot(224)
    x=[w for w in np.linspace(0,1,100)]
    k=1.875
    z=[sin(k*x)-sinh(k*x)+((cosh(k) -cos(k))/(sin(k)-sinh(k)))*(cos(k*x)-cosh(k *x) )for x in np.linspace(0,1,100)]
    plt.plot(z, x, 'g',label='Форма  -%s'%str(k))
    plt.legend(loc='best') 
    plt.grid(True)
    plt.show()
    

    Программа позволяет менять все параметры модели, например, для:
    N=1000, B=0.2, f=1, n1=0, n2=20, w1=5.0, w2=10.0, v1=0, v2=0


    Для отношения частот 0.5 переходной процесс множит фигуры. Поставим “ворота” времени n15=0, n2=20, получим.


    Снимем” ворота” и введём начальную фазу v2=-pi/2, получим:


    С учётом изложенного выше, графики комментарий не требую.

    Для интриги


    Если эта статья найдёт своих читателей или читатели её найдут, не устрашившись теней прошлого, то я опубликую трёхмерные анимационные графики сложных пространственных колебаний трубки при изменении в ней уровня заполняющей жидкости.

    Вместо выводов


    Изобретение Жюля Антуана Лиссажу продолжает свой путь во времени, но уже и на Python. Надеюсь, что представленная интерпретация, конечно далёкая от совершенства, позволит продолжить знакомство с работами гениального математика Лиссажу.

    Ссылки


    1. Биографии учёных физиков.
    2. Что такое фигуры Лиссажу?
    3. Фигуры Лиссажу.
    4. Вибрационный уровнемер.А.С.№777455
    5. Модель колебательного звена с применением символьного и численного решений дифференциального уравнения на SymPy и NumPy.
    6. Модель колебательного звена в режиме резонансных колебаний на Python.
    • +19
    • 6,1k
    • 2
    Поделиться публикацией

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

      0
      Отлично. Как раз ребенку показываю всякие несложные штуки на питоне. Спасибо!
        0
        Занятно, спасибо. Фигуры Лиссажу были моей первой серьёзной программой, ещё на Бейсике.

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

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