Методичка по МС

Рисунок 8.2 — Блок-схема алгоритма моделирования явления радиоактивного распада

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

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

Среднее время жизни ядра можно определить как среднее арифметическое всех времен жизни

	N
	∑τi		(8.4)
τ =	i=1	.
	N

Для подсчета среднего времени жизни в процессе расчетов будем суммировать все времена жизни ядер и после каждого распада будем прибавлять к переменной τp время, ко-

торое	просуществовала
частица	до	распада
τp(i+1) =τpi +τi .		Среднее

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

(8.4)).

Активность вещества вычисляется согласно определению — количество распадов в единицу времени. Мгновенная активность вещества вычисляется по

(8.3).

Период полураспада определяется в процессе расчетов, как время, через

которое распадется половина первоначальных частиц.

Составим блок-схему алгоритма моделирования распада (рисунок 8.2).

8.2 Модель ядерного реактора

Будем считать, что в реакторе выполняются следующие условия:

1.Исходными параметрами является масса ядерного топлива — m , процентное со-

держание U 235 в U 238 — p0 . Масса топлива заведомо превышает критическую массу.

2.В процессе реакции p1 (доля частиц) всех нейтронов покидает реактор, уходя в окружающую среду.

31

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

3. Уровень управления определяется параметром p2 — доля поглощенных графи-

том (управляющие стрежни) нейтронов, выведенных из зоны реакции.

4. Каждый не поглощенный нейтрон с вероятностью p3 p0 вызовет деление ядра

U 235 , а иначе нейтрон поглощается U 238 с образованием Pu239 .

5. В каждом акте деления выделяется 168 МэВ тепловой энергии, 18 МэВ энергии β - и γ -излучения, 2,5 нейтрона (то есть в каждом акте деления обязательно образуется 2 нейтрона и еще один нейтрон с вероятностью 50%).

6.Тепловая энергия должна отводиться водой, при этом температура воды после охлаждения реактора должна быть в пределах от 800 C до 950 C . На реактор вода подается с температурой 10 C . Управление отводом тепла производиться при помощи изменения протока воды Q, м3с.

7.Оболочка реактора позволяет поглотить p4 энергии β - и γ -излучения. Энергия излучения за пределами реактора не должна превышать 0,1 Дж.

8.В процессе работы реактора температура должна находиться в пределах нормы и реакция не должна переходить в область не управляемой реакции. Иметь возможность регулировать температуру жидкости получаемой с реактора.

Алгоритм построения компьютерной модели следующий:

1.Задаем параметры: m — масса топлива, p0 — начальное содержание U 235 , p1 — число нейтронов покинувших реактор, p2 — число нейтронов поглощенных гра-

фитом, p3 — вероятность осуществления реакции захвата нейтрона, λ — посто-

	янная распада (вероятность распада ядра),			dt —		малый промежуток времени,
	Q — начальный расход воды.
	2. Вычисляем количество ядер всех веществ:			N =	m(1− p0 )		N	A	— количество ато-
				1		M
		mp0
	мов U 238 , N2 =		NA — количество ядер	U 235 ,		N3 = 0 — количество атомов
		M
	Pu239 .

3.Определяем количество распавшихся ядер U 235 по закону N = N2e−λdt . Тогда количество нейтронов Nn = 2,5N .

4.Количество нейтронов покинувших область реакции Nn1 = p1Nn .

5.	Количество	нейтронов	поглощенных	управляющими		стержнями
	Nn2 = (1− p1 ) p2 Nn . Итого в области реакции осталось (1− p1 )(1− p2 )Nn нейтронов.
6.	Количество	нейтронов	вступивших	в	ядерную	реакцию
	Nn3 = (1− p1 )(1− p2 ) p3 p0 Nn .
7.	Количество	нейтронов поглощенных U 238		и превративших		его в Pu239

Nn4 = (1− p1 )(1− p2 )(1− p3 p0 )Nn .

8. Выделяется тепловая энергия E =168 МэВ N , после отвода этой энергии при

	расходе воды	Q температура охлаждающей жидкости t	кон	=t	нач	+	E	, где
							Q dt c

c — удельная теплоемкость воды.

32

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

9. Энергия β - и γ -излучения W =18 МэВ N p4 за пределами оболочки реактора. 10.Определяем новые значения количества частиц каждого вещества: N1 = N1 − Nn4 ,

N

2

= N

2

− N

n3

, N

3

= N

3

+ N

n4

и концентрации U 235 в U 238 p = N

N .

0

2 1

11.Переходим к пункту 3.

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

8.3 Скатывание шариков

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

	f (x) =		1	e−	( x−x )2	(8.5)
					2σ2
		σ	2π				Рисунок 8.3

где f (x) — функция распределения, площадь под которой определяет вероятность

того, что результат эксперимента попадает в интервал ∆x (рисунок 8.3), x — наиболее вероятное значение (среднее арифметическое результатов всех экспериментов),

σ2 — дисперсия случайной величины. При проведении эксперимента дисперсию

случайно величины можно определить по формуле

		n
	σ =	∑(xi − x )2	(8.6)
		i=1
		N

Примером событий, определяемых нормальным распределением является ска-

тывание множества шариков по наклонной плоскости, в которую вбито множество

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

	падания шарика в эту точку. Сама вероятность			px попадания шарика в точку с ко-
	ординатой x будет равна	Nx
	px =		,	(8.7)
		N

где Nx — количество шариков упавших в точку с координатой x , N — общее количество скатившихся шариков.

33

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

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

диться только один шарик. При ударе о первый гвоздь с вероятностью 0,5 он будет

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

34

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

8.4 Игры

Смоделируем процесс лотереи схожей по правилам с лотереей «Русское лото». Для этого сгенерируем шесть группы неповторяющихся цифр от одного до 90. В каждой группе должно быть по 5 цифр. Каждая группа соответствует строке лотерейного билета. Первые три строки являются первой карточкой, следующие три —

второй. Набор из этих тридцати цифр представляет собой один билет. Всего таких

билетов должно быть N (не менее 500000)

Согласно правилам игры выигрышным является билет, в котором раньше других выпадут все числа строки, карточки, билета. При этом количество цифр, не ис-

пользованных в лотерее равно шести.

Чтобы определить, какова вероятность выигрыша Вашего билета необходимо запустить процесс лотереи, то есть выбирать по одной цифре от 1 до 90 и сравнивать полученный ряд с рядами в билетах. При этом необходимо проводить подсчет общего числа ходов и выигрышных билетов.

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

35

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

Приложения

Приложение 1. Основные типы данных и часто используемые команды

Pascal 7.0.

	Тип	Диапазон							Размер в	Примечание
									байтах
								Целые типы	*
	ShortInt	-128…127							1	Знаковый тип
	Integer	-32768…32767							2	Знаковый тип
	LongInt	-2147483648…2147483647							4	Знаковый тип
	Byte	0…255							1	Беззнаковый тип
	Word	0…65535							2	Беззнаковый тип
							Вещественные		типы**
	Real	2,9 10−39 …1,7 1038							6	Знаковый тип, число значащих
										цифр: 11–12
	Single***	1,5 10−45 …3,4 1038							4	Знаковый тип, число значащих
										цифр: 7–8
	Double***	5,0 10−324 …1,7 10308							8	Знаковый тип, число значащих
										цифр: 15–16
	Extended***	3,4 10	−4932	…1,1 10	4932				10	Знаковый тип, число значащих
										цифр: 19–20
										Знаковый тип, число значащих
	Comp***	−9,2 1018 …9.2 1018				(	−263	+1…263 +1	8	цифр: 19–20. Фактически явля-
								)		ется типом целых чисел увели-
										ченного диапазона
								Логический	тип
	Boolean	False, True							1	False<True
								Символьный	тип
	Char	Один ASCII-символ							1
								Строковый тип
										Длина строковой переменной
	String	Строка длиной до 255 символов							1…255	(константы) может быть указа-
										на явно String[75].
								Указатель
	Pointer	Ссылка на ячейки памяти с адресами							4	Является адресом величины
		от 0000:0000 до FFFF:FFFF								базового типа

*При описании команд под типом Integer будем понимать все целые типы.

**При описании команд под типом Real будем понимать все вещественные типы.

***Для использования этих вещественных типов необходимо использовать сопроцессор. Подключение производиться при помощи глобального ключа компилятора {$N+} или меню

Options/Compiler/Numeric Processing/[x] 8087/80287. Операции над числами этих типов производятся значительно бастре чем над числами типа Real.

Модуль System.

Abs(x:Real):Real {Abs(x:Integer):Integer} — функция, возвращает абсолютное значение аргумента.

ArcTan(x:Real):Real — возвращает значение арктангенса аргумента в радианах, в пределах от −π / 2 до π / 2 .

Cos(x:Real):Real — возвращает значение косинуса аргумента, x— угол в радианах.

36

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

Exp(x:Real):Real — возвращает значение ex . Frac(x:Real):Real — возвращает дробную часть аргумента. Int(x:Real):Real — возвращает целую часть аргумента. Ln(x:Real):Real — возвращает натуральный логарифм аргумента. Pi:Real — возвращает число π .

Sqr(x:Real):Real {Sqr(x:Integer):Integer} — возвращает квадрат аргумента.

Sqrt(x:Real) — возвращает значение квадратного корня.

Dec(i[,N]:Integer) — уменьшает значение переменной на N, если N не указано, то на 1. Inc(i[,N]:Integer) — увеличивает значение переменной на N, если N не указано, то на 1. Delete(s:String; i,c:Integer) — удаляет из строки s начиная с i-ого c символов.

Insetr(s1:String; s2:String; i:Integer) — помещает подстроку s1 в строку s2,

начиная с i-ой позиции строки s2.

Str(i[:M[:N]]:Integer{Real},s:String) — преобразует число i в последовательность символов в строке s, M, N — форматы вывода числа.

Val(s:String; i:Integer{Real}; c:Integer) — преобразует символьное представле-

ние числа s в двоичное i. В переменной c номер неправильного символа при ошибке, иначе нуль.

Copy(s:String; i,c:Integer):String — возвращает подстроку строки s, выделяя с i-

ого символа с символов.

Length(s:String):Integer — возвращает длину строки s.

Pos(s1,s2:String):Integer — возвращает номер первого символа с которого строка s2 входит в строку s1.

Chr(b:Byte):Char — возвращает символ с указанным кодом.

Round(x:Real):Integer — возвращает округленное до целого значение вещественного числа. Dispose(p:Pointer) — освобождает память, занятую динамической переменной. FreeMem(p:Pointer; w:Word) — удаляет динамическую переменную p размером w, выделенную процедурой GetMem.

GetMem(p:Pointer; w:Word) — создает новую динамическую переменную размером w и помещает ее адрес в p.

New(p:Pointer) — создает новую динамическую переменную p.

ParamCount:Word — возвращает количество параметров командной строки программы. ParamStr(i:Word):String — возвращает указанный параметр командной строки. Randomize — инициализирует генератор случайных чисел.

Random:Real — возвращает случайное число в диапазоне от 0 до 1. SizeOf(x):Word — возвращает размер аргумента в байтах.

Модуль CRT.

TextBackGround(c:Byte) — задает цвет фона, цвета устанавливаются в зависимости от выбранного видеорежима.

TextColor(c:Byte) — задает цвет символов.

GotoXY(x,y:Byte) — устанавливает текстовый курсор в заданную позицию. KeyPressed:Boolean — анализирует нажатие клавиш на клавиатуре, True, если клавиша нажата, иначе False.

ReadKey:Char — считывает символ с клавиатуры.

Delay(i:Word) — задержка на i миллисекунд (зачастую задержка выполняется пропорционально i, но не в миллисекундах).

Модуль Graph.

InitGraph(gd,gm:Integer;s:String) — инициализирует графическую систему, устанавливая графический режим, gd — графический драйвер (часто gd:=Detect — автоопределение

37

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

графического драйвера), gm — номер графического режима, s — путь к драйверу (файл с расширением .bgi).

CloseGraph — завершает работу в графическом режиме и осуществляет переход в текстовый. ClearDevice — очистка графического экрана.

SetActivePage(w:Word) — устанавливает номер активной (для записи информации) страницы. SetVisualPage(w:Word) — устанавливает номер отображаемой на экране страницы. SetBkColor(w:Word) — устанавливает цвет фона экрана.

SetColor(w:Word) — устанавливает текущий цвет для отображения контуров фигур и текста. SetRGBPalette(c,r,g,b:Integer) — устанавливает цвет c в соответствии с установками интенсивности красного r, зленного g и синего b цветов.

GetMaxX:Word — возвращает максимальную координату x для текущего драйвера в выбранном режиме.

GetMaxY:Word — возвращает максимальную координату y для текущего драйвера в выбранном режиме.

GetMaxColor:Word — возвращает максимальный номер цвета для данного драйвера. GetPixel(x,y:Integer) — возвращает цвет пикселя с указанными координатами PutPixel(x,y:Integer;c:Word) — устанавливает пиксель с заданными координатами x,y и цветом c.

Line(x1,y1,x2,y2:Integer) — проводит линию текущего цвета между точками (x1,y1) и (x2,y2). Цвет задается — SetColor, параметры линии — SetLineStyle, режим построе-

ния — SetWriteMode.

LineRel(dx,dy:Integer) — проводит линию с относительными координатами (dx,dy) из текущей точки.

LineTo(x,y:Integer) — проводит линию в точку с абсолютными координатами (x,y) из текущей точки.

MoveRel(dx,dy:Integer) — перемещает указатель в точку с относительными координатами

(dx,dy).

MoveTo(x,y:Integer) — перемещает указатель в точку с абсолютными координатами (x,y). SetLineStyle(ls,p,t:Word) — задает текущие параметры линии, ls — стиль линии, p — шаблон линии, t — толщина линии.

SetWriteMode(i:Integer) —устанавливает режим построения линии. Bar(x1,y1,x2,y2:Integer) — строит прямоугольник, закрашенный текущим орнаментом и цветом заполнения (SetFillStyle, SetFillPattern), (x1,y1) — координаты верхнего левого угла, (x2,y2) — координаты правого нижнего угла прямоугольника. DrawPoly(np:Word; pp) — стрит контур многоугольника с заданными параметрами и цветом (см. Line), np — количество вершин многоугольника, pp — переменная без типа, содержащая np+1 пар вершин многоугольника (массив).

FillPoly(np:Word; pp) — строит закрашенный многоугольник (см. DrawPoly). Rectangle(x1,y1,x2,y2:Integer) — строит контур многоугольника с текущими параметрами (см. Line).

SetFillStyle(p,c:Word) — задает орнамент p и цвет c заполнения фигур. Arc(x,y:Integer;s,e,r:Word) — строит дугу с координатами центра (x,y), s — начальный угол, e —конечный угол, r — радиус.

Circle(x,y:Integer;r:Word) — строит окружность с координатами центра (x,y) и радиусом r. параметры вывода задаются процедурами SetColor и SetLineStyle. Ellipse(x,y:Integer;s,e,xr,yr:Word) — строит дугу эллипса с координатами центра (x,y), начальным углом s, конечным углом e и полуосями xr,yr.

OutText(s:String) — выводит строку s на экран, в текущую позицию курсора.

38

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

OutTextXY(x,y,:Integer;s:String) — выводит строку s на экран начиная с заданных координат x,y.

SetTextStyle(f,d,cs:Word) — устанавливает текущий тип шрифта f, направление вывода текста d и размер символа cs.

GetImege(x1,y1,x2,y2:Word;Map) — сохраняет изображение заданного участка (x1,y1, x2,y2) в переменную Map без типа.

PutImege(x,y:Integer;Map;s:Word) — выдает на экран образ изображения, записанный

вбуфере Map, x,y — координаты верхнего левого угла для вывода, s — способ вывода на экран. ImageSize(x1,y1,x2,y2:Integer):Word — возвращает количество байт для размещения

впамяти изображения заданного размера.

Модуль DOS.

GetIntVec(n:Byte;v:Pointer) — в переменной v возвращает адрес вектора прерывания с номером n.

SetIntVec(n:Byte;v:Pointer) — устанавливает адрес программы обработки прерывания, n — номер прерывания, v — адрес программы-обработчика.

GetTime(h,m,s,s100:word) — возвращает текущее время, установленное в DOS, h, m, s, s100 соответственно часы, минуты, секунды и сотые доли секунд. GetDate(y,m,d,dw:Word) — возвращает текущую дату, установленную в DOS, y, m, d, dw — год, месяц, день месяца и день недели соответственно.

SetTime(h,m,s,s100:Word) — устанавливает системное время (GetTime). SetDate(y,m,d:Word) — устанавливает системную дату (GetDate).

39

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

Приложение 2. Задание движения при помощи матриц

program	Line_And_Sqr_Matrix; {Название программы}
uses	CRT,Graph,Objects;	{Подключаемые модули}
const	N=3;	{Размерность матриц}
	FN=9;	{Количество точек в фигуре}
	dl=20;	{Масштабный коэффициент}
	Figure : array[1..FN] of TPoint=((x:0;y:4), (x:1;y:1), (x:4;y:0),
	(x:1;y:-1),(x:0;y:-4),(x:-1;y:-1),(x:-4;y:0),(x:-1;y:1),(x:0;y:4));
type	{ачальные координаты точек фигуры}		{Определяем тип квадратной матрицы}
	TMatrix=array[1..N,1..N] of Extended;
	Matrix=^TMatrix;		{Указатель на квадратную матрицу}
	TLineMatrix=array[1..N] of Extended;		{Определяем тип линейной матрицы}
var	LineMatrix=^TLineMatrix;		{Указатель на линейную матрицу}
	MatrixA,MatrixB,MatrixC : Matrix;		{Квадратные матрицы для проведения
арифметических операций}		: Matrix;	{Матрица поворота}
	MatrixR
	MatrixM,MatrixN	: LineMatrix;	{Матрицы точек}
	i,j	: Integer;	{Переменные циклов}
	gd,gm	: Integer;	{Переменные для определения графики}
	phi	: Extended;	{Угол поворота фигуры}
	color	: Byte;	{Цвет}
	x,y	: Integer;	{Координаты точек фигуры}
{********************************************************************}
procedure	ShowSqrMatrix(M:Matrix; n1,n2:Byte); {Процедура вывода матрицы на экран}
var i,j : Integer;
begin		{От первой строки до последней}
for j:=1 to N do
begin		{От первого столбца до последнего}
for i:=1 to N do
Write(M^[i,j]:n1:n2,' ');		{Выводим элемент матрицы переданной в виде указателя}
WriteLn;		{Перевод строки}
end;
end; {ShowAqrMatrix}
{********************************************************************}
procedure	AddMatrix(A,B,C:Matrix); {Сложение квадратных матриц}
var i,j : Integer;
begin
for j:=1 to N do
for i:=1 to N do C^[i,j]:=A^[i,j]+B^[i,j];

end; {AddMatrix} {********************************************************************} procedure MulSqrMatrix(A,B,C:Matrix); {Умножение квадратных матриц} begin

{Вписать текст самтостоятельно} end; {MulSqrMatrix}

{********************************************************************}

procedure MulMatrix(A:Matrix;MM,MN:LineMatrix); {Умножение матрицы A на матрицу MM,

результат в MN}

var i,j : Integer; begin

for j:=1 to N do begin MN^[j]:=0;

for i:=1 to N do MN^[j]:=MN^[j]+A^[i,j]*MM^[i]; end;

end; {MulMatrix} {********************************************************************} procedure SetMatrixR(p:Extended); {Заполнение матрицы поворота на угол p} begin

MatrixR^[1,1]:=cos(p); MatrixR^[2,1]:=-sin(p); MatrixR^[3,1]:=0; MatrixR^[1,2]:=sin(p); MatrixR^[2,2]:=cos(p); MatrixR^[3,2]:=0;

40