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

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

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

Рисунок 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. Это симулирует маятник Фуко в течение 3 600 секунд (можно изменить время симуляции). Модель использует относительный допуск по умолчанию RelTol = 1e-6.

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

Результаты

Результаты симуляции показывают в рисунке 2 выше. Симуляция вычисляет координаты X и Y маятника и скоростные компоненты X и Y маятника.

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

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

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

Примечание: "Анимационные Результаты" фрагмент примера требуют Signal Processing Toolbox™. Двойной клик на блоке анимации вызовет ошибку, если это не будет установлено. Все другие части примера будут функционировать правильно без Signal Processing 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 дает положение боба маятника, 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 относительно Омеги.

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

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

Аналитическое (аппроксимированное) решение

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

$\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 и добавьте необходимое выражение. Это изменение введет очень маленькую полную коррекцию числовому результату.

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