Моделирование маятника Фуко

Открыть модель

В этом примере показано, как моделировать маятник Фуко. Маятник Фуко был детищем французского физика Леона Фуко. Предполагалось доказать, что Земля вращается вокруг своей оси. Плоскость колебаний маятника Фуко вращается в течение суток в результате осевого вращения Земли. Плоскость колебаний завершает целый круг во временном интервале T, который зависит от географической широты.

Самый известный маятник Фуко был установлен внутри парижского Пантеона. Это была металлическая сфера 28 кг, прикрепленная к проволоке длиной 67 метров. В этом примере имитируется маятник длиной 67 метров на географической широте Парижа.

Модель Simulink ®

Самый простой способ решения проблемы маятника Фуко в Simulink ® - это создание модели, которая решает связанные дифференциальные уравнения для системы. Эта модель показана на рис. 1. Уравнения, описывающие маятник Фуко, приведены ниже. Подробные сведения о физике модели и выводе этих уравнений см. в разделе Анализ и физика.

$
\ddot{x} = 2\Omega \sin{\lambda} \dot{y} - \frac{g}{L} x
$

$
\ddot{y} = - 2\Omega\sin{\lambda} \dot{x} - \frac{g}{L} y
$

$x, y = \mbox{ pendulum bob coordinates as seen by an observer on Earth}$

$\Omega = \mbox{ Earth's angular speed of revolution about its axis } (rad/sec)$

$g = \mbox{ acceleration of gravity } (m/sec^2)$

$L = \mbox{ length of the pendulum string } (m)$

$\lambda = \mbox{ geographic latitude } (rad)$

Открытие модели

Напечатать sldemo_foucault в окне команд MATLAB ® для открытия этой модели. Эта модель регистрирует данные моделирования в переменнойsldemo_foucault_output. Записанные сигналы имеют синий индикатор. Дополнительные сведения см. в разделе Настройка сигнала для регистрации.

Рис. 1: Маятниковая модель Фуко

Начальные условия

Эта модель загружает константы и начальные условия из sldemo_foucault_data.m файл. Содержимое этого файла показано в таблице 1 ниже. Параметры моделирования можно изменять непосредственно в рабочем пространстве MATLAB. Начальная амплитуда маятника должна быть небольшой по сравнению с длиной маятника, поскольку дифференциальные уравнения действительны только для малых колебаний.

Таблица 1: Исходные условия

g = 9.83;          % acceleration of gravity (m/sec^2)
L = 67;            % pendulum length (m)
initial_x = L/100; % initial x coordinate (m)
initial_y = 0;     % initial y coordinate (m)
initial_xdot = 0;  % initial x velocity (m/sec)
initial_ydot = 0;  % initial y velocity (m/sec)
Omega=2*pi/86400;  % Earth's angular velocity of rotation about its axis (rad/sec)
lambda=49/180*pi;  % latitude in (rad)

Выполнение моделирования

Нажмите кнопку «Воспроизведение» на панели инструментов в окне модели для запуска моделирования. При моделировании будет использоваться жесткий решатель с переменным шагом ode23t. Он будет моделировать маятник Фуко в течение 3600 секунд (можно изменить время моделирования). Модель использует относительный допуск по умолчанию RelTol = 1e-6.

Рисунок 2: Результаты моделирования маятника Фуко (время моделирования 3600 с)

Результаты

Результаты моделирования показаны на фигуре 2 выше. При моделировании вычисляются координаты x и y маятника, а также компоненты скорости x и y маятника.

Плоскость колебаний маятника завершает 360 градусную развертку более чем за 24 часа. Период развертки является функцией географической широты lambda (см. деривацию в анализе и физике).

Рисунок 3: Блок анимации показывает, насколько плоскость колебаний маятника вращается за час

После выполнения моделирования дважды щелкните блок анимации, чтобы анимировать результаты.

Примечание.Часть примера «Анимация результатов» требует Toolbox™ обработки сигналов. Двойной щелчок по блоку анимации приведет к ошибке, если он не установлен. Все остальные части примера будут работать правильно без панели инструментов обработки сигналов.

sldemo_foucault_animate.m файл отображает положение маятника в разные моменты времени. Хорошо видно, как вращается плоскость колебаний маятника.

Примечание.Если моделирование выполняется с большим относительным допуском, результат будет численно нестабильным в течение длительного периода времени. Убедитесь, что используется решатель жестких переменных шагов. Подробнее о численной нестабильности жестких проблем и производительности решателя читайте в примере «Исследование решателей с переменным шагом с использованием жесткой модели».

Закрытие модели

Закройте модель. Очистить сгенерированные данные.

Анализ и физика

В этом разделе анализируется маятник Фуко и описывается физика за ним. Маятник может быть смоделирован как точечная масса, подвешенная на проволоке длиной L. Маятник расположен на географической широте lambda. Удобно использовать опорные кадры, показанные на рис. 4: инерционный кадр I (относительно центра Земли), и неинерционный кадр N (относительно наблюдателя на поверхности Земли). Неинерциальный каркас ускоряется в результате вращения.

Рис. 4: Инерционные и неинерциальные кадры для решения проблемы

Точка O является началом неинерциального каркаса N. Это точка на поверхности земли под точкой подвешивания маятника. Неинерциальный каркас выбран таким образом, что ось Z указывает в сторону от центра Земли и перпендикулярна поверхности Земли. Ось X указывает на юг, а ось Y - на запад.

Как упоминалось во введении, плоскость колебаний маятника Фуко вращается. Плоскость колебаний завершает полный поворот во времени Trot по следующей формуле, где Tday - продолжительность одного дня (т.е. время, в течение которого Земля вращается вокруг своей оси один раз).

$T_{rot}=T_{day} \cdot \sin \lambda$

Синусоидальный фактор требует дальнейшего обсуждения. Часто неверно предполагается, что плоскость колебаний маятника закреплена в инерциальном каркасе относительно центра Земли. Это справедливо только на северном и южном полюсах. Чтобы устранить эту путаницу, подумайте о точке S (см. рис. 4), где маятник подвешен. В инерционном кадре I точка S перемещается по окружности. Маятниковый боб подвешен на проволоке постоянной длины. Для простоты игнорируйте трение воздуха. В инерционном кадре I есть только две силы, которые действуют на боб - натяжение проволоки T и гравитационной силы Fg.

Вектор r В показано положение маятника В (см. рис. 4). Второй закон Ньютона гласит, что сумма всех сил, действующих на тело, равна массе, умноженной на ускорение тела.

$m \ddot{\overrightarrow{r}} = \overrightarrow{T} + \overrightarrow{F_g}$

На протяжении всего этого доказательства точки обозначают производные времени, стрелки - векторы, шапки - унитарные векторы (i, j и k по осям x, y и z). Точка над векторной стрелкой указывает производную вектора по времени. Стрелка над точкой указывает вектор производной времени. См. разницу между общим ускорением и радиальным ускорением ниже.

Общее ускорение:

$
\ddot{\overrightarrow{r}} = \frac{d^2 \overrightarrow{r}}{d t^2}=
\frac{d^2}{d t^2} \left( x\mathbf{\hat{i}} + y\mathbf{\hat{j}} + z\mathbf{\hat{k}} \right)
$

Радиальное ускорение:

$
\overrightarrow{ \ddot{r} } = \overrightarrow{ \left( \frac{d^2 r}{dt^2}\right)}=
\ddot{x} \mathbf{\hat{i}}+
\ddot{y} \mathbf{\hat{j}}+
\ddot{z} \mathbf{\hat{k}}
$

Ускорение силы тяжести указывает к центру земли (отрицательное z-направление).

$$g=9.83 m/sec^2$$

$\overrightarrow{g} = -g\mathbf{\hat{k}}$

$
m \ddot{\overrightarrow{r}} =
\overrightarrow{T} - mg\mathbf{\hat{k}}
$

Разложить член ускорения:

$
\overrightarrow{r}=
r_x \mathbf{\hat{i}}+
r_y \mathbf{\hat{j}}+
r_z \mathbf{\hat{k}}
$

$
\dot{\overrightarrow{r}}=
\left(
\dot{r}_x \mathbf{\hat{i}}+
\dot{r}_y \mathbf{\hat{j}}+
\dot{r}_z \mathbf{\hat{k}}
\right) +
r_x \dot{ \mathbf{\hat{i}} } +
r_y \dot{ \mathbf{\hat{j}} } +
r_z \dot{ \mathbf{\hat{k}} }
$

Производные времени единичных векторов появляются, потому что неинерциальный опорный кадр N вращается в пространстве. Это означает, что унитарные векторы i, j и k вращаются в пространстве. Их производные времени приведены ниже. Омега - угловая скорость вращения Земли вокруг своей оси. Скалярный Омега - это значение угловой скорости. Векториал Омега - векторная угловая скорость. Его направление определяется правилом правой руки.

$
\dot{\mathbf{\hat{i}}}=\overrightarrow{\Omega} \times \mathbf{\hat{i}}
$

$
\dot{\mathbf{\hat{j}}}=\overrightarrow{\Omega} \times \mathbf{\hat{j}}
$

$
\dot{\mathbf{\hat{k}}}=\overrightarrow{\Omega} \times \mathbf{\hat{k}}
$

Перепишите производную по времени вектора r относительно Омеги.

$
\dot{\overrightarrow{r}}=
\left(
\dot{r}_x \mathbf{\hat{i}}+
\dot{r}_y \mathbf{\hat{j}}+
\dot{r}_z \mathbf{\hat{k}}
\right) +
\overrightarrow{\Omega} \times \overrightarrow{r}=
\overrightarrow{\dot{r}}+
\overrightarrow{\Omega} \times \overrightarrow{r}
$

Аналогично, выразить вторую производную по времени вектора r.

$
\ddot{\overrightarrow{r}}=
\overrightarrow{\ddot{r}} +
2 \overrightarrow{\Omega} \times \overrightarrow{\dot{r}} +
\overrightarrow{\Omega} \times
\left( \overrightarrow{\Omega} \times \overrightarrow{r} \right)
$

$\overrightarrow{\dot{r}} = \mbox{ acceleration in the non-inertial frame N (x, y, z components)}$

$2 \overrightarrow{\Omega} \times \overrightarrow{\dot{r}} = \mbox{ Coriolis acceleration}$

$\overrightarrow{\Omega} \times
\left( \overrightarrow{\Omega} \times \overrightarrow{r} \right)
= \mbox{ additional term due to rotation of non-inertial frame N}
$

Чтобы упростить это уравнение, предположим, что Омега для Земли очень мала. Это позволяет игнорировать третий член в уравнении выше. На самом деле второй член (который уже намного меньше первого члена) на четыре порядка больше третьего члена. Это сводит уравнение к следующей форме:

$
\ddot{\overrightarrow{r}} \simeq
\overrightarrow{\ddot{r}} +
2 \overrightarrow{\Omega} \times \overrightarrow{\dot{r}}
$

Второй закон Ньютона можно записать и разложить на компоненты x, y и z следующим образом:

$
m \overrightarrow{\ddot{r}} =
\overrightarrow{T} - mg\mathbf{\hat{k}} -
2 m \overrightarrow{\Omega} \times \overrightarrow{\dot{r}}
$

$
m \ddot{x} = T_x + 2m\Omega \dot{y} \sin{\lambda}
$

$
m \ddot{y} = T_y - 2m\Omega \left(\dot{x} \sin{\lambda}+\dot{z}\cos{\lambda}\right)
$

$
m \ddot{z} = T_z - mg + 2m \Omega \dot{y} \cos{\lambda}
$

Угловая амплитуда колебаний мала. Поэтому мы можем игнорировать вертикальную скорость и вертикальное ускорение (z-точка и z-двойная точка). Компоненты натяжения струны могут быть выражены с помощью малых угловых аппроксимаций, которые также значительно упрощают задачу, делая её двумерной (см. ниже).

$T_z = mg - 2m\Omega \dot{y} \cos{\lambda} \simeq mg$

$T_x= -T\frac{x}{L}$

$T_y= -T\frac{y}{L}$

$T_z= T\frac{L-z}{L}\simeq T$

Характеристические дифференциальные уравнения

Наконец, физика задачи может быть описана системой связанных уравнений, приведенных ниже. Координаты x и y определяют положение маятникового боба, видимое наблюдателем на Земле.

$
 \ddot{x} - 2\Omega \sin{\lambda} \dot{y} + \frac{g}{L} x =0
$

$
 \ddot{y} + 2\Omega\sin{\lambda} \dot{x} + \frac{g}{L} y =0
$

Аналитическое решение (приблизительное)

Ниже приведено аналитическое решение проблемы маятника Фуко. К сожалению, она не точна. Если попытаться подставить аналитическое решение в дифференциальные уравнения, останутся несансированные члены порядка Omega в квадрате. Однако, поскольку Омега очень мала, мы можем игнорировать неизбранные термины в практических целях.

$\eta = x+ i\cdot y \mbox{ (complex number)}$

$\ddot{\eta}+(2i\Omega \sin{\lambda})\dot{\eta} + \frac{g}{L} \eta = 0$

$\eta = \left( C_1 e^{i\sqrt{g/L}t} + C_2 e^{-i\sqrt{g/L}t}\right)
e^{-i\Omega t \sin{\lambda}}$

$C_1, C_2 \mbox{ are complex integration
constants}$

Фактическая система дифференциальных уравнений асимметрична

Во время деривации термины с участием Омеги в квадрате игнорировались. Это привело к симметрии xy в дифференциальных уравнениях. Если принять во внимание квадратные члены Омеги, система дифференциальных уравнений становится асимметричной (см. ниже).

$
\ddot{x} - 2\Omega \sin{\lambda} \dot{y} + (\frac{g}{L}-\Omega^2 \sin^2{\lambda}) x =0
$

$
\ddot{y} + 2\Omega\sin{\lambda} \dot{x} + (\frac{g}{L} - \Omega^2) y =0
$

Можно легко изменить текущую модель маятника Фуко, чтобы учесть асимметричные дифференциальные уравнения. Просто отредактируйте соответствующие блоки усиления, которые содержат g/L и добавьте необходимое выражение. Это изменение внесет очень небольшую общую коррекцию в численный результат.

Документация