КП: Молекула углекислого газа
Курсовой проект по Теоретической механике
Исполнитель: Смирнов Александр
Группа: 09 (23604)
Семестр: весна 2015
Содержание
Аннотация проекта
Моделирование молекул, даже самых простых - сложная задача. Для их моделирования необходимо использовать многочастичные потенциалы, но их поограммирование - тоже очень сложная задача. Встает вопрос о том, можло ли найти более простой путь моделирования простейших молекул.
Для моделирования хорошо подходят парные потенциалы, ибо они имеют простой вид и легко программируются. Но как их применить к моделирования молекул ? Моя работа и посвещена решению данной проблемы.
Формулировка задачи
Смоделировать с помощью многочастичного потенциала молекулу углекислого газа ( 2D модель ) и рассмотреть ее простейшую динамику частицы.
Общие сведения по теме
Для решения задачи будем использовать два потенциала: потенциал Морзе и потенциал Леннард-Джонса.
Потенциал Морзе
Парный силовой потенциал взаимодействия. Определяется формулой:
где
- — энергия связи,
- — длина связи,
- — параметр, характеризующий ширину потенциальной ямы.
Потенциал имеет один безразмерный параметр
. При взаимодействия Морзе и Леннард-Джонса близки. При увеличении ширина потенциальной ямы для взаимодействия Морзе уменьшается, взаимодействие становится более жестким и хрупким. Уменьшение приводит к противоположным изменениям — потенциальная яма расширяется, жесткость падает. Сила, соответствующая потенциалу Морзе, вычисляется по формуле:Или в векторной форме:
Потенциал Леннард-Джонса
Также парный силовой потенциал взаимодействия. Определяется формулой:
где
- — расстояние между частицами,
- — энергия связи,
- — длина связи.
Сила взаимодействия, соответствующая потенциалу Леннард-Джонса, вычисляется по формуле
Векторная сила взаимодействия определяется формулой
Углекислый газ
Углекислый газ ( диоксид углерода ) - газ без запаха и цвета. Молекула углекислого газа имеет линейное строение и ковалентные полярные связи, хотя сама молекула не является полярной. Дипольный момент = 0.
Решение
Взяв за основу программу Balls v4, было получено:
- программа, в которой можно рассмотреть динамику одной молекулы углекислого газа в другом газе ( например в азоте ).
- 2 версии программы с различными потенциалами.
- создание молекулы за конечное время.
Первая версия (потенциал Морзе)
- Синий шар - кислород .
- Коралловый шар - углерод .
- Желтый шар - азот.
Вторая версия (потенциал Морзе)
- Цвет шаров аналогичный.
- Левая кнопка мыши - добавление желтой частицы.
- Правая кнопка мыши - удаление любого шара.
- В системе теперь присутствует несколько молекул.
Файл CO2.js <syntaxhighlight lang="javascript" line start="1" enclose="div"> // m: Вязкоупругий шар // Версия 4.5 от 05.11.2014
window.addEventListener("load", MainBalls, true); function MainBalls() {
// Предварительные установки
var canvas = canvasBalls; var context = canvas.getContext("2d"); // на context происходит рисование canvas.oncontextmenu = function (e) {return false;}; // блокировка контекстного меню
var Pi = 3.1415926; // число "пи"
var m0 = 1; // масштаб массы var T0 = 1; // масштаб времени (период колебаний исходной системы) var a0 = 1; // масштаб расстояния (диаметр шара)
var g0 = a0 / T0 / T0; // масштаб ускорения (ускорение, при котором за T0 будет пройдено расстояние a0) var k0 = 2 * Pi / T0; // масштаб частоты var C0 = m0 * k0 * k0; // масштаб жесткости var B0 = 2 * m0 * k0; // масштаб вязкости
// *** Задание физических параметров ***
var Ny = 5; // число шаров, помещающихся по вертикали в окно (задает размер шара относительно размера окна) var m = 1 * m0; // масса var Cwall = 10 * C0; // жесткость стен var B = 0.003 * B0; // вязкость среды var Bwall = 0.03 * B0; // вязкость на стенках var Cball = 0.1 * Cwall; // жесткость между частицами
var mg = 0//0.25 * m * g0; // сила тяжести
var r = 0.1 * a0; // радиус частицы в расчетных координатах
var K = 0.85; // сила взаимодействия ограничивается значением, реализующимся при r/a = K
var a = 5 * r; // равновесное расстояние между частицами
var a_1 = 2.1 * r;
var a_2 = 2.1 * a ;
var aCut = 2 * r ; // радиус обрезания var alfa = 5 ;
// *** Задание вычислительных параметров ***
var fps = 50; // frames per second - число кадров в секунду (качечтво отображения) var spf = 100; // steps per frame - число шагов интегрирования между кадрами (скорость расчета) var dt = 0.045 * T0 / fps; // шаг интегрирования (качество расчета)
// Выполнение программы
var scale = canvas.height / Ny / a0 ; // масштабный коэффициент для перехода от расчетных к экранным координатам var r2 = r * r ; // ___в целях оптимизации___ var aCut2 = aCut * aCut ; // ___в целях оптимизации___ var a2 = a * a ; // ___в целях оптимизации___ var a22 = a_2 * a_2 ; var a11 = a_1 * a_1 ;
var D = a2 * Cball / 72 ; // энергия связи между частицами
var LJCoeff = 12 * D / a2 ; // коэффициент для расчета потенциала Л-Дж
var a1 = alfa / a ; var MorzCoeff = 2 * a1 * D ;
var Ka = K * r ; // ___в целях оптимизации___ var K2a2 = Ka*Ka ; // ___в целях оптимизации___
var w = canvas.width / scale ; // ширина окна в расчетных координатах var h = canvas.height / scale ; // высота окна в расчетных координатах
var dNd = null ; // ссылка на захваченный курсором шар (drag & drop)
// Работа с мышью
var mx_, my_; // буфер позиции мыши (для расчета скорости при отпускании шара)
canvas.onmousedown = function(e) { // функция при нажатии клавиши мыши var m = mouseCoords(e); // получаем расчетные координаты курсора мыши // цикл в обратную сторону, чтобы захватывать шар, нарисованный "сверху" // (т.к. цикл рисования идет в обычном порядке) for (var i = balls.length - 1; i >= 0; i--) { var b = balls[i]; var rx = b.x - m.x; var ry = b.y - m.y; var rLen2 = rx * rx + ry * ry; // квадрат расстояния между курсором и центром шара if (rLen2 <= r2) { // курсор нажал на шар if (e.which == 1) { // нажата левая клавиша мыши dNd = b; dNd.xPlus = dNd.x - m.x; // сдвиг курсора относительно центра шара по x dNd.yPlus = dNd.y - m.y; // сдвиг курсора относительно центра шара по y mx_ = m.x; my_ = m.y; canvas.onmousemove = mouseMove; // пока клавиша нажата - работает функция перемещения } else if (e.which == 3) { // нажата правая клавиша мыши balls.splice(i, 1); // удалить шар } return; } }
for (var i = C.length - 1; i >= 0; i--) {
var b = C[i]; var rx = b.x - m.x; var ry = b.y - m.y; var rLen2 = rx * rx + ry * ry; // квадрат расстояния между курсором и центром шара if (rLen2 <= r2) { // курсор нажал на шар if (e.which == 1) { // нажата левая клавиша мыши dNd = b; dNd.xPlus = dNd.x - m.x; // сдвиг курсора относительно центра шара по x dNd.yPlus = dNd.y - m.y; // сдвиг курсора относительно центра шара по y mx_ = m.x; my_ = m.y; canvas.onmousemove = mouseMove; // пока клавиша нажата - работает функция перемещения } else if (e.which == 3) { // нажата правая клавиша мыши C.splice(i, 1); // удалить шар } return; } }
for (var i = O1.length - 1; i >= 0; i--) {
var b = O1[i]; var rx = b.x - m.x; var ry = b.y - m.y; var rLen2 = rx * rx + ry * ry; // квадрат расстояния между курсором и центром шара if (rLen2 <= r2) { // курсор нажал на шар if (e.which == 1) { // нажата левая клавиша мыши dNd = b; dNd.xPlus = dNd.x - m.x; // сдвиг курсора относительно центра шара по x dNd.yPlus = dNd.y - m.y; // сдвиг курсора относительно центра шара по y mx_ = m.x; my_ = m.y; canvas.onmousemove = mouseMove; // пока клавиша нажата - работает функция перемещения } else if (e.which == 3) { // нажата правая клавиша мыши O1.splice(i, 1); // удалить шар } return; } }
for (var i = O2.length - 1; i >= 0; i--) {
var b = O2[i]; var rx = b.x - m.x; var ry = b.y - m.y; var rLen2 = rx * rx + ry * ry; // квадрат расстояния между курсором и центром шара if (rLen2 <= r2) { // курсор нажал на шар if (e.which == 1) { // нажата левая клавиша мыши dNd = b; dNd.xPlus = dNd.x - m.x; // сдвиг курсора относительно центра шара по x dNd.yPlus = dNd.y - m.y; // сдвиг курсора относительно центра шара по y mx_ = m.x; my_ = m.y; canvas.onmousemove = mouseMove; // пока клавиша нажата - работает функция перемещения } else if (e.which == 3) { // нажата правая клавиша мыши O2.splice(i, 1); // удалить шар } return; } }
// если не вышли по return из цикла - нажатие было вне шара, добавляем новый if (e.which == 1) { dNd = addNewBall(m.x, m.y); // добавляем шар и сразу захватываем его курсором if (dNd == null) return; // если шар не добавился (из за стен или других шаров) - возвращаемся dNd.xPlus = 0; dNd.yPlus = 0; // держим шар по центру mx_ = m.x; my_ = m.y; canvas.onmousemove = mouseMove; // пока клавиша нажата - работает функция перемещения }
};
document.onmouseup = function(e) { // функция при отпускании клавиши мыши canvas.onmousemove = null; // когда клавиша отпущена - функции перемещения нету dNd = null; // когда клавиша отпущена - захваченного курсором шара нету };
function mouseMove(e) { // функция при перемещении мыши, работает только с зажатой ЛКМ var m = mouseCoords(e); // получаем расчетные координаты курсора мыши dNd.x = m.x + dNd.xPlus; dNd.y = m.y + dNd.yPlus; dNd.vx = 0.6 * (m.x - mx_) / dt / fps; dNd.vy = 0.6 * (m.y - my_) / dt / fps; mx_ = m.x; my_ = m.y; }
function mouseCoords(e) { // функция возвращает расчетные координаты курсора мыши var m = []; var rect = canvas.getBoundingClientRect(); m.x = (e.clientX - rect.left) / scale; m.y = (e.clientY - rect.top) / scale; return m; }
// Работа с массивом
var balls = []; // массив шаров
var C =[] ; var O1 = [] ; var O2 = [] ;
var addNewBall = function(x, y) {
// проверка - не пересекается ли новый шар со стенами или уже существующими шарами if (x - r < 0 || x + r > w || y - r < 0 || y + r > h) return null; for (var i = 0; i < balls.length; i++) { var rx = balls[i].x - x; var ry = balls[i].y - y; var rLen2 = rx * rx + ry * ry; if (rLen2 < 4 * r2) return null; } var b = [];
b.x = x; b.y = y; // расчетные координаты шара b.fx = 0; b.fy = mg; // сила, действующая на шар b.vx = 0; b.vy = 0; // скорость
balls[balls.length] = b; // добавить элемент в конец массива
return b;
};
var addNewC = function(x, y) { // проверка - не пересекается ли новый шар со стенами или уже существующими шарами if (x - r < 0 || x + r > w || y - r < 0 || y + r > h) return null; for (var i = 0; i < balls.length; i++) { var rx = balls[i].x - x; var ry = balls[i].y - y; var rLen2 = rx * rx + ry * ry; if (rLen2 < 4 * r2) return null; } var b = [];
b.x = x; b.y = y; // расчетные координаты шара b.fx = 0; b.fy = mg; // сила, действующая на шар b.vx = 0; b.vy = 0; // скорость
C[C.length] = b; // добавить элемент в конец массива
return b;
};
var addNewO1 = function(x, y) { // проверка - не пересекается ли новый шар со стенами или уже существующими шарами if (x - r < 0 || x + r > w || y - r < 0 || y + r > h) return null; for (var i = 0; i < balls.length; i++) { var rx = balls[i].x - x; var ry = balls[i].y - y; var rLen2 = rx * rx + ry * ry; if (rLen2 < 4 * r2) return null; } var b = [];
b.x = x; b.y = y; // расчетные координаты шара b.fx = 0; b.fy = mg; // сила, действующая на шар b.vx = 0; b.vy = 0; // скорость
O1[O1.length] = b; // добавить элемент в конец массива
return b;
};
var addNewO2 = function(x, y) {
// проверка - не пересекается ли новый шар со стенами или уже существующими шарами if (x - r < 0 || x + r > w || y - r < 0 || y + r > h) return null; for (var i = 0; i < balls.length; i++) { var rx = balls[i].x - x; var ry = balls[i].y - y; var rLen2 = rx * rx + ry * ry; if (rLen2 < 4 * r2) return null; } var b = [];
b.x = x; b.y = y; // расчетные координаты шара b.fx = 0; b.fy = mg; // сила, действующая на шар b.vx = 0; b.vy = 0; // скорость
O2[O2.length] = b; // добавить элемент в конец массива
return b;
};
// Основной цикл программы
function control() { physics(); draw(); }
// Расчетная часть программы
function powers(b ,b2 , hkk , jk ) { var rx = b.x - b2.x; var ry = b.y - b2.y; // вектор смотрит на первый шар (b)
var r2 = rx * rx + ry * ry; // квадрат расстояния между шарами
var rLen = (Math.sqrt(r2));
if (r2 < K2a2) {
if (rLen > 0.00001) { // проверка, чтобы избежать деления на 0 rx = rx / rLen * Ka; ry = ry / rLen * Ka; } r2 = K2a2; rLen = Ka; // корень K2a2 }
// сила взаимодействия
var u = Math.exp( -a1 * ( rLen - hkk )) ; var F = jk * MorzCoeff * u* ( u - 1 ) /rLen ;
var Fx = F * rx; var Fy = F * ry; b.fx += Fx; b.fy += Fy; b2.fx -= Fx; b2.fy -= Fy;
}
function physics() { // то, что происходит каждый шаг времени for (var s = 1; s <= spf; s++) {
// пересчет сил идет отдельным массивом, т.к. далее будут добавляться силы взаимодействия между шарами
for (var i0 = 0; i0 < balls.length; i0++) {
balls[i0].fx = - B * balls[i0].vx; balls[i0].fy = mg - B * balls[i0].vy; }
for (var i0 = 0; i0 < C.length; i0++) {
C[i0].fx = - B * C[i0].vx; C[i0].fy = mg - B * C[i0].vy; }
for (var i0 = 0; i0 < O1.length; i0++) {
O1[i0].fx = - B * O1[i0].vx; O1[i0].fy = mg - B * O1[i0].vy; }
for (var i0 = 0; i0 < O2.length; i0++) {
O2[i0].fx = - B * O2[i0].vx; O2[i0].fy = mg - B * O2[i0].vy; }
for (var i = 0; i < balls.length; i++) { // dlya Balls // расчет взаимодействия производится со всеми следующими шарами в массиве, // чтобы не считать каждое взаимодействие дважды
var b = balls[i];
for (var j = i + 1; j < balls.length; j++) {
var b2 = balls[j]; var rx = b.x - b2.x; var ry = b.y - b2.y; // вектор смотрит на первый шар (b) var r2 = rx * rx + ry * ry; // квадрат расстояния между шарами
if ( r2 > aCut2 ) continue ;
powers(b,b2,a_1,0.3) ;
} for (var j = 0; j < C.length; j++) {
var b2 = C[j]; var rx = b.x - b2.x; var ry = b.y - b2.y; // вектор смотрит на первый шар (b) var r2 = rx * rx + ry * ry; // квадрат расстояния между шарами
if ( r2 > aCut2 ) continue ;
powers(b,b2,a_1,0.3) ;
} for (var j = 0; j < O1.length; j++) {
var b2 = O1[j]; var rx = b.x - b2.x; var ry = b.y - b2.y; // вектор смотрит на первый шар (b) var r2 = rx * rx + ry * ry; // квадрат расстояния между шарами
if ( r2 > aCut2 ) continue ;
powers(b,b2,a_1,0.3) ;
}
for (var j = 0; j < O2.length; j++) {
var b2 = O2[j]; var rx = b.x - b2.x; var ry = b.y - b2.y; // вектор смотрит на первый шар (b) var r2 = rx * rx + ry * ry; // квадрат расстояния между шарами
if ( r2 > aCut2 ) continue ;
powers( b, b2, a_1 ,0.3) ;
}
if (b == dNd) continue; // если шар схвачен курсором - его вз. со стенами и перемещение не считаем
if (b.y + r > h) { b.fy += -Cwall * (b.y + r - h) - Bwall * b.vy; } if (b.y - r < 0) { b.fy += -Cwall * (b.y - r) - Bwall * b.vy;} if (b.x + r > w) { b.fx += -Cwall * (b.x + r - w) - Bwall * b.vx; } if (b.x - r < 0) { b.fx += -Cwall * (b.x - r) - Bwall * b.vx; }
b.vx += b.fx / m * dt; b.vy += b.fy / m * dt; b.x += b.vx * dt; b.y += b.vy * dt;
}
for (var i = 0; i < C.length; i++) { // dlya C
// расчет взаимодействия производится со всеми следующими шарами в массиве, // чтобы не считать каждое взаимодействие дважды
var b = C[i];
for (var j = i + 1; j < C.length; j++) {
var b2 = C[j]; var rx = b.x - b2.x; var ry = b.y - b2.y; // вектор смотрит на первый шар (b) var r2 = rx * rx + ry * ry; // квадрат расстояния между шарами
if ( r2 > aCut2 ) continue ;
powers(b,b2,a_1, 0.3) ;
} for (var j = 0; j < O1.length; j++) {
var b2 = O1[j]; var rx = b.x - b2.x; var ry = b.y - b2.y; // вектор смотрит на первый шар (b) var r2 = rx * rx + ry * ry; // квадрат расстояния между шарами
if ( j == i ) { powers(b,b2,a, 6) ;} else { if ( r2 > aCut2 ) continue ; powers(b,b2,a_1, 0.3) ; }
}
for (var j = 0; j < O2.length; j++) {
var b2 = O2[j]; var rx = b.x - b2.x; var ry = b.y - b2.y; // вектор смотрит на первый шар (b) var r2 = rx * rx + ry * ry; // квадрат расстояния между шарами
if ( j == i ) { powers(b,b2,a,6) ;} else { if ( r2 > aCut2 ) continue ; powers(b,b2,a_1, 0.3) ; }
}
if (b == dNd) continue; // если шар схвачен курсором - его вз. со стенами и перемещение не считаем
if (b.y + r > h) { b.fy += -Cwall * (b.y + r - h) - Bwall * b.vy; } if (b.y - r < 0) { b.fy += -Cwall * (b.y - r) - Bwall * b.vy;} if (b.x + r > w) { b.fx += -Cwall * (b.x + r - w) - Bwall * b.vx; } if (b.x - r < 0) { b.fx += -Cwall * (b.x - r) - Bwall * b.vx; }
b.vx += b.fx / m * dt; b.vy += b.fy / m * dt; b.x += b.vx * dt; b.y += b.vy * dt;
}
for (var i = 0; i < C.length; i++) { // dlya O1
// расчет взаимодействия производится со всеми следующими шарами в массиве, // чтобы не считать каждое взаимодействие дважды
var b = O1[i];
for (var j = i + 1; j < O1.length; j++) {
var b2 = O1[j];
var rx = b.x - b2.x; var ry = b.y - b2.y; // вектор смотрит на первый шар (b) var r2 = rx * rx + ry * ry; // квадрат расстояния между шарами
if ( r2 > aCut2 ) continue ;
powers(b,b2,a_1,0.3) ;
}
for (var j = 0; j < O2.length; j++) {
var b2 = O2[j]; var rx = b.x - b2.x; var ry = b.y - b2.y; // вектор смотрит на первый шар (b) var r2 = rx * rx + ry * ry; // квадрат расстояния между шарами
if ( i == j ) { powers(b,b2,a_2, 3) ;} else { if ( r2 > aCut2 ) continue ; else powers(b,b2,a_1 , 0.3) ; }
}
if (b == dNd) continue; // если шар схвачен курсором - его вз. со стенами и перемещение не считаем
if (b.y + r > h) { b.fy += -Cwall * (b.y + r - h) - Bwall * b.vy; } if (b.y - r < 0) { b.fy += -Cwall * (b.y - r) - Bwall * b.vy;} if (b.x + r > w) { b.fx += -Cwall * (b.x + r - w) - Bwall * b.vx; } if (b.x - r < 0) { b.fx += -Cwall * (b.x - r) - Bwall * b.vx; }
b.vx += b.fx / m * dt; b.vy += b.fy / m * dt; b.x += b.vx * dt; b.y += b.vy * dt;
}
for (var i = 0; i < C.length; i++) { // dlya O2
// расчет взаимодействия производится со всеми следующими шарами в массиве, // чтобы не считать каждое взаимодействие дважды
var b = O2[i];
for (var j = i + 1; j < O2.length; j++) {
var b2 = O2[j]; var rx = b.x - b2.x; var ry = b.y - b2.y; // вектор смотрит на первый шар (b) var r2 = rx * rx + ry * ry; // квадрат расстояния между шарами
if ( r2 > aCut2) continue ;
powers(b,b2,a_1, 0.3) ;
}
if (b == dNd) continue; // если шар схвачен курсором - его вз. со стенами и перемещение не считаем
if (b.y + r > h) { b.fy += -Cwall * (b.y + r - h) - Bwall * b.vy; } if (b.y - r < 0) { b.fy += -Cwall * (b.y - r) - Bwall * b.vy;} if (b.x + r > w) { b.fx += -Cwall * (b.x + r - w) - Bwall * b.vx; } if (b.x - r < 0) { b.fx += -Cwall * (b.x - r) - Bwall * b.vx; }
b.vx += b.fx / m * dt; b.vy += b.fy / m * dt; b.x += b.vx * dt; b.y += b.vy * dt;
} } }
// Рисование
var rScale13 = r * scale * 1.3; // ___в целях оптимизации___ var rScaleShift = r * scale / 5; // ___в целях оптимизации___ function draw() { context.clearRect(0, 0, w * scale, h * scale); // очистить экран for (var i = 0; i < balls.length; i++){ var xS = balls[i].x * scale; var yS = balls[i].y * scale; context.fillStyle = "#FFD700";
context.beginPath(); context.arc(xS, yS, r * scale, 0, 2 * Math.PI, false); context.closePath(); context.fill(); }
for (var i = 0; i < C.length; i++){
var xS = C[i].x * scale; var yS = C[i].y * scale; context.fillStyle = "#f08080";
context.beginPath(); context.arc(xS, yS, r * scale, 0, 2 * Math.PI, false); context.closePath(); context.fill(); }
for (var i = 0; i < O1.length; i++){
var xS = O1[i].x * scale; var yS = O1[i].y * scale; context.fillStyle = "#3070d0";
context.beginPath(); context.arc(xS, yS, r * scale, 0, 2 * Math.PI, false); context.closePath(); context.fill(); }
for (var i = 0; i < O2.length; i++){
var xS = O2[i].x * scale; var yS = O2[i].y * scale; context.fillStyle = "#3070d0";
context.beginPath(); context.arc(xS, yS, r * scale, 0, 2 * Math.PI, false); context.closePath(); context.fill(); }
}
// Запуск системы for (var i = 0; i < 8; i++){ // добавляем 20 частиц, сдвинув их от стен //addNewBall(Math.random() * (w - 2 * r) + r, Math.random() * (h - 2 * r) + r , 10 );
addNewC(Math.random() * (w - 2 * r) + r, Math.random() * (h - 2 * r) + r , 10 ); addNewO1(Math.random() * (w - 2 * r) + r, Math.random() * (h - 2 * r) + r , 10 ); addNewO2(Math.random() * (w - 2 * r) + r, Math.random() * (h - 2 * r) + r , 10 );
}
setInterval(control, 1000 / fps);
} </syntaxhighlighting>
Обсуждение результатов и выводы
Можно смоделировать молекулу с помощью парных потенциалов. За конечное время при небольшой вязкости среды частицы собираются в молекулы. При моделировании частиц только с помощью парных потенциалов молекула получается очень неустойчивой. Но ее поведение в системе при небольших скоростях получается достаточно правдаподобной. Это означает что можно изучать динамику системы при малых скоростях частиц. При больших скоростях частицы разрушаются. После разрушения они собираются заного, образуя те же самые частицы.
Скачать презентацию: Молекула углекислого газа в среде инертного газа