36424

Компьютерное моделирование физических процессов

Лабораторная работа

Информатика, кибернетика и программирование

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

Русский

2013-09-21

161.5 KB

35 чел.

Лабораторная работа 7

Компьютерное моделирование физических процессов

Материалы для заданий 1, 2

Потери пучка при прохождении через вещество

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

Мы будем изучать, естественно, очень простой вариант задачи — прохождение пучка тяжелых частиц (A) через слой газа, состоящего из легких частиц (O). К тому же будем полагать скорости частиц газа до столкновения с частицами пучка равными нулю. Не будем также учитывать изменений, происходящих в газе вследствие прохождения пучка. Это не ограничивает существенно общности задачи.

1. Эффективные сечения

Частицы можно представлять шариками радиусов Ra и Ro. Пусть частица A летит так, что центр ее должен пролететь на расстоянии ρ от центра частицы O. Столкновение произойдет, если ρ < R, где R = Ra + RO. Фактически частица О образует для частиц A «преграду», площадь которой σ = πR2 . Эта величина называется эффективным сечением столкновения . В зависимости от величины р (называемой прицельным параметром) определяются угол отклонения, передача энергии от частицы A частице O и т.п. Соответственно можно определить величины эффективных сечений, например, для потери частицей A определенной энергии от ε1 до ε2, отклонения на угол, больший данного и т.п. Определяются и дифференциальные эффективные сечения, например, дифференциальное эффективное сечение потери энергии частицей A

определяется так, что величина площадки  = f(ε) отвечает потере энергии в интервале от ε до ε +  (при достаточно малой величине ).

Характерные размеры пучка частиц — сантиметры или даже микроны — во много раз превосходят характерные прицельные параметры. Именно поэтому здесь уместен статистический подход. Пусть концентрация газа-мишени составляет nO частиц O на 1 см3 . Тогда в очень тонком слое (толщины dx ) на площадь S приходится nOS dx частиц O, а перекрываемая ими площадь равна dS = σnOSdx. Можно сказать также, что dW = dS/S = σnOdx представляет долю площади, перекрытую частицами O в слое dx. Иначе говоря, это вероятность столкновения частицы в данном слое. Имея в виду возможность выбрать достаточно тонкий слой, мы можем пренебречь возможностью, что какая-то из частиц O будет «затенена» другими. Фактически для этого достаточно условия dW <C 1.

(Для столкновений атомов характерные величины эффективных сечений имеют порядок 10 см2 , для нейтронов и атомных ядер — 10 см2 , на современных ускорителях изучаются сечения еще на 10 - 12 порядков меньшие.)

2. Потери частиц пучка при прохождении слоя

Начнем с простейшей задачи о выбывании частиц из пучка. Пусть любое столкновение ведет к выбыванию частицы — сечение выбывания равно σ . Задача об ослаблении пучка легко решается в таком случае аналитически. Пусть количество частиц A в исходном пучке равно Na0 . После прохождения слоя dx это число уменьшится на dNA = NAdW. Введя число частиц, оставшихся в пучке, Na(x), имеем

dNA(x) = −NA(x)nOσdx.                  (1)

Решение этого дифференциального уравнения, удовлетворяющее условию , т.е.

                               (2)

где b = nOσ. Величина 1/b определяет толщину слоя, в котором пучок потеряет значительное число частиц (ослабнет в e =2.7 раза). Приведенное решение задачи относится к средним значениям числа частиц.

Теперь покажем, как решается эта задача методом Монте-Карло. Будем «пропускать» частицу через тонкие слои толщиной dx, пока она не пройдет слой x (или поглотится). «Розыгрыш» для слоя dx состоит в том, что мы выбираем случайное число с равномерным распределением от 0 до 1 (функция rand в MATLAB) и сравниваем с величиной dW; если rand< dW, то частица выбывает, в противном случае перемещается в следующий слой. Пропустив Na частиц через слой толщины x, мы можем получить число оставшихся частиц Na(x).

Приведем отрезок программы, в котором через слой толщины Xmax «пропускается» Na частиц. При реализации этой части программы предполагается, что для всех частиц заведены два вектора - вектор координат, достигнутых каждой из частицей к текущему моменту времени, и вектор состоящий из 0 и 1. Каждая 1 соответствует еще не выбывшей частице, а 0 - уже выбывшей из потока.

Следующая программа описывает моделирование выбывания частиц из пучка (листинг А).

% Программа расчета потерь пучка %

%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear;

Na=250; k=(1:Na)'; no=100; SecA=0.05; dx=0.01;

x=zeros(size(k)); % Начальный вектор координат

sc=ones(size(k)); % Начальный вектор-счетчик

dW=no*SecA*dx;          % Вероятность поглощения на шаге dx
hl=line(x,k); % Подготовка Na точек для рисования

set(hl, 'Marker', 'o', 'MarkerSize',3);
axis([0 1 0 Na+1]);          % Масштабирование осей
pause; % Пауза перед запуском основного цикла

% Далее выполняется цикл, на каждом шаге которого

% координата x возрастает на dx до тех пор, пока

% частица не пройдет путь Xmax или не поглотится

while (any(sc)> 0)

 ra=rand(size(k)); % Расчет вероятностей с помощью

% датчика случайных чисел для

% всех частиц

 k1=find(ra-dW< 0);       % Определение номеров частиц

% подлежащих отбраковке

 sc(k1)=0; % Зануление счетчика отбракованной частицы

 x=x+sc*dx; % Продвижение на dx остальных частиц

 set(hl,'XData',x); % Рисование

end; % Конец цикла while

%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

Число частиц Nf, прошедших через достаточно толстый слой, оказывается малым и сильно флуктуирует (изменяется от запуска к запуску). Распределение вероятностей для доли прошедших частиц p = Nf/Naэто распределение Бернулли.

3. Превращение частиц пучка при прохождении слоя

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

Приведем участок программы, моделирующей «судьбу» частиц в этом случае (листинг B).

% Моделирующая часть программы с поглощением A   %
% и превращением A -> B %

% SecA - сечение поглощения A

% SecAB - сечение превращения A в B

% scA - счетчик числа частиц типа A

% scB - счетчик числа частиц типа B

clear; % Очистка рабочей области

Na=25; % Начальное число частиц A

L=5; k=(1:Na)'; no=100;dx=0.01; SecA=0.1; SecAB=0.15;

x=zeros(size(k)); % Начальный вектор координат

scA=ones(size(k)); % Вектор-счетчик частиц A

scB=zeros(size(k)); % Вектор-счетчик частиц B

dWa=no*dx*SecA; % Вероятность поглощения A на шаге dX

dWab=no*dx*SecAB;        % Вероятность A -> B на шаге dX

% Имитация поглощения и превращения частиц

% Пока есть частицы A и не пройден путь L

while (any(scA)> 0 & all(x)< L) ra=rand(size(k));

 ka=find(ra-dWa< 0);       % Номера поглощенных частиц A
 
scA(ka)=0; % Выбывание частиц типа A

 kb=find(dWa< ra & ra< dWa+dWab)

 scB(kb)=scA(kb);          % Превращение A в B
 
scA(kb)=0; % Исчезновение этих A

 x=x+scA*dx+scB*dx;

 pause(1);

end;

Для подобной модели несложно также составить и решить уравнения, описывающие изменения среднего числа частиц A и B, подобно уравнениям (1), (2).

Стоит также отметить близость рассматривавшейся задачи к задаче о радиоактивном распаде атомных ядер (для перехода к ней нужно понимать под x время, прошедшее с начала наблюдения).

Материалы для задания 3

Метод Монте-Карло при нахождении значения определенного интеграла

Во многих исследованиях для построения модели некоторого физического процесса необходимо решать задачи нахождения значений определенных интегралов. Рассмотрим метод статистических испытаний (метод Монте-Карло) для решения указанной задачи.

Рассмотрим метод Монте-Карло применительно к задаче нахождения значения определенного интеграла

.

Используем известный факт, что для нахождения значения определенного интеграла необходимо найти площадь криволинейной трапеции, ограниченной отрезком [a,b] оси OX и графиком функции f(x) на этом отрезке.

Для простоты изложения пусть известна аналитическая запись подынтегральной функции f(x) и задан интервал интегрирования [a,b].

Решение задачи состоит в следующем:

на основании равномерного распределения  случайных величин (например, функция rand в MatLab) случайным образом выбираем координаты точки (x,y),

где

Если

,

то данная точка (x,y) вносит вклад в значение определенного интеграла (такие точки образуют множество M1).

Если

,

то точка (x,y) не засчитывается при нахождении значения определенного интеграла (такие точки образуют множество M2).

Следовательно, множество M1 образуют точки, находящиеся ниже графика функции f(x), множество M2 образуют точки, находящиеся выше графика функции f(x).

Обозначим, S1 – площадь криволинейной трапеции, т.е. площадь фигуры, расположенной ниже графика f(x), S2 – площадь фигуры, расположенной над графиком функции f(x).

Рис. 1. Фигуры S1 и S2, ограниченные графиком функции f(x)

Пусть после проведения серии N испытаний в множестве M1 имеем N1 точек, а в множестве M2 имеем N2 точек.

Тогда, отношение площадей фигур S1 и S2 и отношение количества точек N1 и N2 связаны следующим равенством

.

Сумма площадей S1 и S2 равняется площади прямоугольника S0 со сторонами (b-a) и  ymax.

Следовательно, так как

и

то

,

,

.     (3)

Формула (3) позволяет находить значение определенного интеграла методом Монте-Карло.

Материалы для задания 4

Маятник

Задачи о колебаниях встречаются во всех областях физики. Во многом колебания совершенно различных физических объектов сходны друг с другом. Простейшие примеры — малые колебания маятника и электрические колебания в цепи, составленной из конденсатора и катушки.

Такое движение маятника хорошо известно — это гармонические колебания. Закон движения можно записать в виде x = acos(ω0t +0), где ω0 – частота колебаний, a — амплитуда, 0  — начальная фаза. (Угол отклонения маятника мы обозначили здесь x). Малые колебания описываются уравнением

линейным относительно функции x, поэтому их обычно называют линейными.

В этом задании будем исследовать движение математического маятника при больших углах отклонения.

Свободные колебания

Результаты исследования движения маятника удобно представить в виде набора кривых на плоскости (x,p), где x - угол отклонения маятника, p = скорость изменения угла.

Плоскость (x,p) называется фазовой плоскостью, переменная p — импульсом, а кривые, определяемые параметрически законом движения как x = x(t), p = p(t), — фазовыми траекториями.

Фазовая траектория определяется, например, начальными значениями координаты x(0) и импульса p(0).

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

Угол отклонения маятника достаточно задавать в некоторых конечных пределах, например, принимая − π ≤ x < π.

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

При некоторой энергии колебания сменяются вращением.

Свободное движение математического маятника без трения описывается дифференциальным уравнением

       (4)

Имеется аналитическое решение этого уравнения (довольно сложное), однако мы будем исследовать движение маятника численно.

Запишем это дифференциальное уравнение в виде системы уравнений первого порядка

,

.    (5)

Основная идея численного расчета чрезвычайно проста (метод Эйлера): зная значения координаты и скорости в момент времени t, можно приближенно найти их значения через малый промежуток ∆t

,

,

взяв значения величин х и р из уравнений (5). Многократно повторяя такие вычисления, можно найти зависимости x(t) и p(t). Большей точности можно достигнуть, практически не усложняя расчеты, если использовать так называемую вычислительную схему с перешагиванием. В этой схеме вычисляются значения координат в моменты времени

t−1/2∆t,     t + 1/2∆t,     t + 3/2∆t,     t + 5/2∆t,...,

а значения скоростей в моменты времени

t, t + ∆t,      t + 2∆t,      t + 3∆t,...

Данная вычислительная схема и использована в предлагаемой ниже программе.

Чтобы отклонение от точного решения уравнения (1) было небольшим, должен быть достаточно малым шаг по времени ∆t. Точность вычислений можно проверить, применяя метод повторного счета с уменьшенным шагом ∆t. Если при этом решение остается прежним, значит, шаг был выбран достаточно малым и можно быть уверенным в правильности результата.

Далее приводится текст простейшей программы на языке MATLAB.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Фазовая траектория математического маятника     %

clear; % Очистка рабочей области

% Задание начальных значений

x1=2.2; % Начальная координата

p1=0.0; % Начальный импульс

dt=0.025; % Шаг по времени

axis([-pi pi -pi pi]); % Задание диапазона осей

hl=line(x1,p1); % Задание дескриптора линии

% Задание параметров выводимой линии

% 'EraseMode',’none' - это режим вывода без стирания

% 'LineStyle',':'  - вывод линии в виде пунктира

% по умолчанию выводится сплошная линия

% 'Color',’r’      - задание цвета линии

set(hl,'EraseMode', 'none', 'LineStyle',':','Color', 'r');

grid on; % Задание вывода координатной сетки

pause; % Пауза, обеспечивает немедленный вывод

% рисунка на экран для продолжения программы
% необходимо нажать любую клавишу
while 1 % Бесконечный цикл

% Основной алгоритм расчета

x2=x1+p1*dt;

p2=p1-sin(x2)*dt;

% Склейка граничных условий

if x2> pi

x2=x2-2*pi; end; if x2< -pi

x2=x2+2*pi; end;

% Вывод очередного участка фазовой траектории

set(hl,'XData',x2,'YData',p2);

% Переприсвоение начальных значений

x1=x2;p1=p2;

end;

(Здесь шаг по времени ∆t обозначен dt, pi=π.)

Следует сделать несколько замечаний. После запуска программа откроет графическое окно и нарисует в нем одно окно с системой координат, осями и координатной сеткой, после чего перейдет в режим ожидания (команда pause). Для продолжения расчета необходимо нажать любую клавишу, и программа начнет вывод фазовой траектории. Поскольку используется бесконечный цикл, выполнение программы можно прервать клавишей Ctr+C. Для повторного прогона с другими параметрами необходимо ввести соответствующие изменения в текст программы и перезапустить ее.

Задания лабораторной работы

Написать файл-функции для решения поставленных далее задач.

Задачи:

1. Построить модель прохождения частицами пучка некоторого слоя заданной ширины и имеющим заданное эффективное сечение поглощения, используя метод Монте-Карло. Условия выполнения приведены далее.

2. Построить модель (используя метод Монте-Карло) прохождения частицами пучка некоторого слоя заданной ширины и имеющим заданное эффективное сечение поглощения с учетом превращения частиц А пучка в частицы В (эффективное сечение превращения задано). Условия выполнения приведены далее.

В задачах 1 и 2 для каждого варианта определены следующие значения:

- ширина слоя Xmax,

- количество частиц в исходном пучке Na,

- ширина элементарного слоя dx,

- концентрация частиц слоя (преграды) no,

- эффективное сечение поглощения (A) SecA,

- эффективное сечение превращения (AB) SecB.

Используя текст листингов А и В, необходимо выполнить следующие задания:

- В задачах 1 и 2 необходимо визуализировать процесс прохождения частицами слоя и превращения частиц А пучка в частицы В.

- В задаче 2 необходимо визуализировать процесс превращения частиц А пучка в частицы В.

- Ширина слоя в задачах 1 и 2 не должна превышать величины Xmax.

3. Найти значение определенного интеграла методом Монте-Карло.

Построить график функции f(x) и отметить на графике «случайные» точки, образованные при решении.

4. Построить фазовую траекторию математического маятника. Визуализируйте процесс движения самого математического маятника.

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

Варианты:

Значения переменных для задания 1, 2.

Вариант

Xmax

Na

dx

no

A

AB

1

4

30

0.01

100

0.05

0.1

2

5

40

0.01

100

0.07

0.08

3

6

45

0.01

100

0.08

0.17

4

7

35

0.01

100

0.09

0.15

5

8

30

0.01

100

0.1

0.2

6

7

35

0.01

100

0.05

0.15

7

6

40

0.01

100

0.1

0.25

8

5

30

0.01

100

0.07

0.3

9

4

45

0.01

100

0.01

0.15

10

7

35

0.01

100

0.15

0.2

Значения переменных для задания 3.

Для всех вариантов – количество испытаний N=2000.

Вариант

f(x)

a

b

ymax

1

sin(x)+2

0

3

max f(x) на [a,b]

2

cos(x)+3

1

5

max f(x) на [a,b]

3

x2-3x+10 

2

6

max f(x) на [a,b]

4

sin(2x)+1

3

7

max f(x) на [a,b]

5

cos(x/2)+5

0

4

max f(x) на [a,b]

6

5x2-8x+14 

1

3

max f(x) на [a,b]

7

sin(4x)+2

2

5

max f(x) на [a,b]

8

cos(x/3)+15

3

7

max f(x) на [a,b]

9

7x2-x+2 

1

5

max f(x) на [a,b]

10

sin(3x)+5

2

8

max f(x) на [a,b]

Значения переменных для задания 4.

Вариант

x1

Импульс p1

dt

1

3

5

0.02

2

2

1

0.02

3

2.5

1

0.02

4

3.1

0.1

0.02

5

3

0

0.02

6

2.5

7

0.02

7

0

3

0.02

8

1

4

0.02

9

0.5

5

0.02

10

3

0

0.02

PAGE  1


 

А также другие работы, которые могут Вас заинтересовать

24074. Строение гемоглобина 82 KB
  Строение гемоглобина. В молекуле гемоглобина белковый компонент представлен белком глобином небелковый компонент гем. За счет еще одной координационной связи к атому железа может присоединяться молекула кислорода с образованием оксигемоглобина.
24075. Синтез гема 41 KB
  При восстановлении биливердина НАДФ Н2 образуется билирубин. Билирубин плохо растворимое соединение и в крови связывается с альбумином. В виде комплекса альбуминбилирубин идет транспорт билирубина кровью в клетки печени. В печени билирубин соединяется с глюкуроновой кислотой с образованием моно 20 и диклюкуронидов 80 они хорошо растворимы в воде.
24076. Распад гема 27 KB
  остающихся в сыворотке крови после осаждения белков. в сыворотке крови является ценным диагностическим показателем при многих заболеваниях. определяют в надосадочной жидкости после удаления осажденных белков сыворотки крови центрифугированием с помощью азотометрического метода Кьельдаля в его многочисленных модификациях колориметрических и гипобромитных методов. в сыворотке крови равна 2040 мг 100 мл или 143286 ммоль л.
24077. Витамины и Антиметаболиты 54.26 KB
  Согласно современным предтавлениям все клеточные и внутриклеточные мембраны устроены сходным образом: основу мембраны составляет двойной молекулярный слой липидов липидный бислой на котором и в толще которого находятся белки см. В состав липидов мембран входят в основном фосфолипиды сфингомиелины и холестерин. Например в мембранах эритроцитов человека их содержание составляет соответственно 36 30 и 22 по весу; еще 12 приходится на гликолипиды Примером амфифильной молекулы может служить молекула фосфатидилэтаноламина структура...
24078. Биохимия печени 32.5 KB
  Биохимия печени Печень самый крупный из паренхиматозных органов. Роль печени в метаболизме углеводов Печень играет ведущую роль в поддержании физиологической концентрации глюкозы в крови. При физиологической гипогликемии в печени активируется распад гликогена. В печени активно протекает глюконеогенез при котором предшественниками глюкозы являются пируват и аланин поступающий из мышц глицерол из жировой ткани и с пищей ряд глюкогенных АК.
24079. Метаболизм белков 35 KB
  Детоксицирующая функция печени Детоксикация ядовитых метаболитов и чужеродных соединений ксенобиотиков протекает в гепатоцитах в две стадии. Реакции первой стадии катализируются монооксигеназной системой компоненты которой встроены в мембраны эндоплазматического ретикулума. На первой стадии биотрансформации происходит образование или высвобождение гидрокси карбоксильных тиоловых и аминогрупп которые являются гидрофильными и молекула может подвергаться дальнейшему превращению и выведению из организма. Кроме цх Р450 в первой...
24080. Биологическая ценность белков 30 KB
  Для оценки состояния обмена белков используется понятие азотистый баланс. Азот остается в организме и расходуется на синтез белков. Встречается при голодании белковой недостаточности тяжелых заболеваниях когда происходит интенсивный распад белков тела. Биологическая ценность белков.
24081. Переваривание белков. Пути превращения аминокислот в печени 105 KB
  Расщепление белков происходит при участии нескольких групп ферментов: Экзопептидазы катализирует разрыв концевой пептидной связи с образованием одной какойлибо аминокислоты. В результате расщепления образуются свободные аминокислоты которые затем подвергаются всасыванию. Аминокислоты всасываются свободно с ионами натрия. Некоторые аминокислоты обладают способностью конкурентно тормозить всасывание других аминокислот: Лизин тормозит всасывание аргинина.
24082. Токсическое действие аммиака-инактивация альфа-кетоглутарата в цикле кребса,энергетическое голодание,к которому чувствителна очень нервная ткань 57.5 KB
  Возможны 4 типа дезаминирования: Восстановительное RCHCOOH RCH2COOH NH3 NH2 Гидролитическое RCHCOOH RCHCOOH NH3 NH2 OH Внутримолекулярное RCH2CHCOOH RCH=CHCOOH NH3 NH2 Окислительное RCHCOOH RCCOOH NH3 NH2 O Окислительное дезаминирование бывает 2 видов: прямое и непрямое трансдезаминирование. R R1 R R1 HCNH2 C=O C=O HCNH2 COOH COOH COOH COOH Реакция трансаминирования...