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

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

В этом примере показано, как смоделировать маятник Фуко. Маятник Фуко был детищем французского физика Леона Фуко. Она была призвана доказать, что Земля вращается вокруг своей оси. Плоскость колебаний маятника Фуко вращается в течение дня в результате осевого вращения Земли. Плоскость колебаний завершает целый круг за временной интервал 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)

Выполнение симуляции

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

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

Результаты

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

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

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

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

Примечание: Фрагмент «Animate Results» примера требует Signal Processing Toolbox™. Двойной клик по блоку анимации приведет к ошибке, если он не установлен. Все другие части примера будут функционировать правильно без Signal Processing Toolbox.

The 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 есть только две силы, которые действуют на bob - натяжение провода T и силу тяготения Fg.

Векторная r задает положение маятника bob, B (см. фигуру 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 относительно Omega.

$
\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-double-точка). Компоненты натяжения строки могут быть выражены с помощью малых угловых приближений, которые также значительно упрощают задачу, делая ее двумерной (см. ниже).

$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 определяют положение маятника bob, как видно наблюдателем на Земле.

$
 \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
$

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

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