Размещено на http://www.allbest.ru/

Размещено на http://www.allbest.ru/

КУРСОВА РОБОТА

з дисципліни «Обчислювальна техніка та програмування»

на тему

Моделювання руху тіла у в'язкому середовищі



Зміст

Вступ

Розділ 1

1.1 Рух тіла у в'язкому середовищі

1.2 В'язкість (внутрішнє тертя) і в'язкопружність

1.3 Актуальність данної роботи

1.4 Поставлена задача

Розділ 2

2.1 Метод Рунге -Кутти

2.2 Інструмент. Мова програмування Сі

2.3 Алгоритм вирішення

2.4 Структура програми

Розділ 3. Результати роботи програми

Висновки

Список використаних джерел

Додаток А. Текст програми

Додаток Б. Блок схеми програми



Вступ

Моделювання є важливим елементом більшості сучасних досліджень і процесів розробки. Однак, якщо раніше моделювання позначало тільки створення зменшеного або збільшеного макета досліджуваної системи або її частини, то до справжнього моменту набули поширення також математичні методи моделювання.

Суть їх полягає в побудові фізичної моделі системи, створення для неї різних параметрів і початкових умов і знаходженні стану цієї системи в заданих умовах. Однак, найчастіше модель виходить досить складною, і рішення рівнянь, що описують її поведінку, виявляється неможливо отримати в точному вигляді. Саме в цих випадках використовується чисельне рішення рівнянь, тобто відбувається чисельне моделювання системи.

Методи чисельного розв'язання рівнянь (зокрема, диференціальних та інтегральних) відомі давно, але тільки з розвитком комп'ютерів стало можливо використовувати їх стосовно до складних систем рівнянь і отримання їх рішень з високою точністю, необхідної в реальних завданнях.

Зміст і рішення поставленої задачі

У цій роботі досліджується рух твердого тіла, що знаходиться всередині в'язкої рідини. Виробляється чисельне рішення рівняння руху і порівнюється з точним рішенням, отриманим математичними методами. Для чисельного рішення диференціального рівняння використовується метод Рунге-Кутти 4-ого порядку. Програма написана на мові програмування Сі.



Розділ 1

1.1 Рух тіла у в'язкому середовищі

Розглянемо занурене у в'язку рідину тверде тіло (для простоти - кулькообразної форми). На це тіло діють наступні сили:

сила земного тяжіння

сила Архімеда з боку рідини

і сила внутрішнього тертя, яка при досить малих швидкостях руху виявляється прямо пропорційна швидкості і спрямована протилежно напрямку руху. При цьому коефіцієнт опору середовища k прямо пропорційний в'язкості середовища і площі дотику тіла і рідини. Для тіл сферичної форми радіуса R він буде дорівнює . Таким чином, сила внутрішнього тертя дорівнює

Рівнодіюча цих сил буде дорівнювати

Рух буде рівномірним, тому, якщо направити вісь OY вертикально вниз, то в проекції на цю вісь рівняння матиме вигляд:

Можна привести його до вигляду задачі Коші

з початковими умовами , . Вважаючи праву частину рівняння рівним f (y, v), можна розв'язати його чисельними методами, що й було зроблено в даній роботі. Крім того, можна знайти точне рішення рівняння, щоб перевірити точність методу Рунге-Кутти.

Якщо ввести позначення

рівняння можна записати у вигляді

або

Проінтегрувавши, отримуємо

Або

Враховуючи, що , рівняння записується у вигляді

Оскільки в початковий момент часу швидкість тіла дорівнює нулю, то остаточно для швидкості отримуємо рівняння

Проінтегрувавши даний вираз за часом і враховуючи, що в початковий момент тіло перебувало на початку координат , знаходимо вираз для залежності координати тіла від часу

1.2 В'язкість (внутрішнє тертя) і в'язкопружність

У фізичної моделі ідеальної рідини передбачається, що в рідкому середовищі відсутнє внутрішнє тертя, тобто взаємодія між елементами середовища, спрямована уздовж поверхні їхнього зіткнення. Однак, у реальній рідині така взаємодія присутня, і така властивість називається в'язкістю. В'язкість обумовлена ??тим, що дотичні шари рідини, що рухаються з різними швидкостями, за рахунок дифузії обмінюються молекулами, а як наслідок - імпульсом. Шар, що рухається з більшою швидкістю, втрачає частину своїх швидких молекул, і натомість отримує повільні молекули, а шар, що рухається з меншою швидкістю, навпаки - втрачає повільні і отримує швидкі молекули. За рахунок цього швидкість руху першого шару знижується, а швидкість руху другого - зростає.

Властивість в'язкопружності деяких рідин означає, що вони можуть вести себе і як пружні тіла, і як в'язкі рідини. При цьому якщо з боку зовнішнього середовища на них виробляється вплив, протягом досить короткого часу, то рідина веде себе, як пружне тверде тіло, і тільки через деякий час після початку дії починає рухатися, як в'язка рідина.

1.3 Актуальність данної роботи

Існуючі рішення

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

Виникає рівняння виду

з початковою умовоюназивається задачею Коші. При чисельному рішенні необхідно знайти значення функції в заданій точці. При цьому послідовно шукається значення в точках множини xn, (n = 1, ... N), де xn як правило і є потрібна нам точка.

Найпростіший з таких методів це метод Ейлера. У ньому весь відрізок від x0 до xn розбивається на рівні відрізки довжини h, і наближено покладається, що



звідки по рекурентного співвідношення

знаходиться значення yn функції u у всіх точках аж до xn.

Однак, метод Ейлера є досить приблизним і має малу точність. Існують модифікації цього методу, названі виправленим методом Ейлера і модифікованим методом Ейлера. Вони мають більш високу точність, проте існує сімейство ще більш точних методів, так званих методами Рунге-Кутти. При цьому згадані три методи виявляються частими випадками методу Рунге-Кутти.

Так, метод Ейлера являє собою метод Рунге-Кутти першого порядку, а його модифікації - методами Рунге-Кутти другого порядку. У даній роботі використовується найпоширеніший з представників цього сімейства, який є методом Рунге-Кутти четвертого порядку (одним з багатьох) і називається просто методом Рунге-Кутти. Він є, з одного боку, досить простим в реалізації, в порівнянні з аналогічними більш високого порядку, з іншого боку - має хорошу точність, достатню для більшості випадків.

1.4 Поставлена задача

У цій роботі вирішується наступна задача. Нехай є сферичне тверде тіло радіуса R і щільності ?т. Необхідно розглянути рух цього тіла у в'язкій рідині щільності ? рід. з коефіцієнтом в'язкості ?, вважаючи, що в початковий момент тіло знаходиться на початку координат і має нульову швидкість. При цьому враховується дія на тіло наступних сил: сила тяжіння, сила Архімеда і сила внутрішнього тертя рідини.

Необхідно скласти рівняння руху тіла в цій системі, написати на мові програмування Сі програму для вирішення цього диференціального рівняння методом Рунге-Кутти. Також необхідно знайти точне рішення рівняння і порівняти обидва рішення.

Програма повинна приймати наступні вхідні параметри: параметри фізичної моделі - радіус тіла, густина тіла, щільність рідини, коефіцієнт в'язкості рідини; параметри методу Рунге-Кутти - кінцевої точки, в якій і обчислюється значення, і величина кроку h.



Розділ 2

2.1 Метод Рунге-Кутти

У загальному вигляді сімейство методів Рунге-Кутти для одного диференціального рівняння першого порядку полягає в наступному. Нехай є задача Коші

для х > x0, u (x0 )= u0

і нехай чисельне рішення задачі Коші відомо в точці xn, а саме

y (x n )= y n. Тоді m-етапний метод Рунге-Кутти (метод Рунге-Кутти m-го порядку) буде полягати в наступному. Задаються деякі числові коефіцієнти a i, b ij, i = 2, 3, ... m; j = 1, 2, ... (m-1), ? l, l = 1, 2, ... m і послідовно обчислюються наступні функції:

k 1 = f (x n, y n )

k 2 = f (x n + a 2 h, y n + b 21 hk 1)

k 3 = f (x n + a 3 h, y n + b 31 hk 1 + b 32 hk 1)

...

k m = f (x n + a m h, y n + b m1 hk 1 + ... + b m, m-1 hk m-1)

Потім з формули

знаходиться рішення задачі в наступній точці (x n+1): y (x n+1) = y n+1.

При цьому потрібно, щоб виконувалась умова



Методи порядку 5-го і вище використовуються рідко, найчастіше застосовуються методи четвертого порядку. При цьому найбільш популярним з них, є так званий метод Рунге-Кутти:

k 1 = f (x n, y n );

k 2 = f (x n + h / 2, y n + k 1 h / 2);

k 3 = f (x n + h / 2, y n + k 2 h / 2);

k 4 = f (x n + h, y n + k 13 h);

2.2 Інструмент. Мова програмування Сі

Мова програмування Сі є одним з найпопулярніших мов високого рівня. Її особливостями є висока ефективність, універсальність, відносна простота, а також широка підтримка з боку розробників програмного забезпечення. Існує величезна кількість компіляторів мови для більшості платформ, а також безліч бібліотек, що розширюють функціональність програм, написаних на ньому.

Завдяки цим особливостям Сі є найпопулярнішою мовою системного програмування, а також широко використовується при створенні прикладного програмного забезпечення.

Мова Сі справила великий вплив на розвиток мов програмування. Так, вона послужила основою для багатьох мов, у тому числі для популярних мов C + + і нащадків останньої мови (наприклад, C #); також багато відомих мови запозичили синтаксис мови Сі.

На даний момент Сі є стандартом де-факто швидкодії та ефективності серед мов програмування високого рівня.



2.3 Алгоритм вирішення

Алгоритм вирішення поставленого завдання по моделюванню руху твердого тіла у в'язкому рідини заснований на рішенні рівняння руху методом Рунге-Кутти. Рівняння записується у вигляді (явної залежності від часу немає)

Безпосередньо застосовуючи метод Рунге-Кутти, можна знайти значення швидкості тіла в будь-який момент часу t, якщо відомо початкове значення (в даній задачі воно покладається рівним нулю).

Однак, для повного опису руху тіла необхідно знати не його швидкість, а координати. При точному вирішенні диференціального рівняння руху знайти формулу для координат дуже просто - потрібно тільки проінтегрувати формулу для швидкості за часом. При чисельному рішення це зробити неможливо.

З іншого боку, швидкість є похідна за часом від координати

Це рівняння теж можна вирішити методом Рунге-Кутти, оскільки ми вже знаємо значення v у потрібні моменти часу. Однак, треба враховувати, що в методі використовується значення функції в моменти, кратні не тільки довжині інтервалу, але і половині інтервалу. Тому, задаючи інтервал часу при рішенні рівняння для координат, необхідно перш за все задати вдвічі менший інтервал при вирішенні рівняння для швидкості, так як для знаходження значення координат у N моментів часу буде потрібно значення швидкості в 2N моментів часу.

2.4 Структура програми

Програма складається з чотирьох функцій:

1. Функція main - основна функція програми, в який відбувається введення параметрів алгоритму (максимального моменту часу, для якого потрібно знайти координати тіла і інтервалу між двома сусідніми моментами часу) і виклик всіх інших функцій;

2. Функція ReadData, в якій відбувається зчитування фізичних параметрів задачі - густин тіла і рідини, радіуса тіла і в'язкості рідини. Ознайомитися зі структурою функції можна у додатку Б (рис. 4);

3. Функція RungeCuttaV, в якій за допомогою методу Рунге-Кутта чисельно вирішується рівняння для швидкості тіла; при цьому повертається масив зі значеннями швидкості в потрібні нам моменти часу; (рис. 5)

4. Функція RungeCuttaY, в якій аналогічно вирішується рівняння для координати тіла - при цьому необхідні значення швидкості у відповідні моменти часу виходять, як результат роботи функції RungeCuttaV, яку ми викликаємо. (рис. 5)

Крім того, у функціях RungeCuttaV і RungeCuttaY відбувається виведення на екран одержуваних значень швидкості і координати в кожен з моментів часу. Одночасно з рішенням, отриманим чисельно, виводиться рішення для того ж моменту часу, отримане з точних формул. Це дозволяє оцінити, наскільки точним виходить рішення, отримане чисельними методами.

Варто відзначити, що оскільки в методі Рунге-Кутти для отримання значення функції в точці xn необхідні значення похідної (функції f (x, u) в задачі Коші) в точках xn і xnh / 2, тобто вдвічі більшому числі точок, то функція RungeCuttaV, чиї результати використовуються у функції RungeCuttaY, виконується з удвічі меншим кроком h між моментами часу.

У програмі підключені слідуючи бібліотеки:

#include <windows.h> - файл бібліотеки для коректної роботи в операційному середовищі Windows.

#include <stdio.h> - файл вводу та виводу стандартної бібліотеки мови Сі, у якому визначені функції та макроси для вводу та виводу даних у стилі мови Сі.

#include <stdlib.h> - файл вводу, який містить у собі різні стандартні функції для виділення пам'яті та контролю процесів. У даній програмі його треба підключати для використання функції malloc, яка виділяє пам'ять.

#include <math.h> - файл стандартної бібліотеки, підключаємо для використання математичних операцій та маніпуляцій. Містить у собі математичні функції та константи. У даній програмі потрібен для використання функції exp

#include <conio.h> - файл для консольного вводу та виводу, потрібен для роботи с текстовою інформацією у консолі. Потрібен для коректної роботи функції getch()



Розділ 3. Результат роботи програми

Щоб оцінити роботу програми, вирішимо чисельно, а також за допомогою точної формули конкретну задачу.

Розглянемо падіння сталевої кульки (щільності 7800 кг / м 3 і радіуса 1 см) в гліцерині (щільність 1261 кг / м 3, в'язкість при температурі 20 о С - 1.49 Па * с), в момент часу t 0 = 0 (тіло знаходиться у спокої на початку координат). Знайдемо його координати в момент часу 1секунда після початку руху.

Задамо для методу Рунге-Кутти інтервал часу, рівний однієї десятитисячної секунди (0,0001 с).

Значення, отримане методом Рунге-Кутти, складає 0,845197 м, а точне значення - 0,845149 м. Різниця між цими величинами складає 0,000048 - всього 0,057%.

Для наочності розглянемо графіки значень координат при значно меншому числі кроків (і відповідно, меншою точності). Нехай інтервал часу буде дорівнювати 0,05 с.

Рис. 1 Графік координат при меншому числі кроків



Значення координат, отриманих для моменту часу 1сек при цьому інтервалі (і, отже, 20 кроках), складає 0,818847 м. Різниця з точним значенням становить 0,021508, тобто 2,54%.

Результат роботи програми з такими параметрами

Рис. 2 Результат роботи програми

Цікаво відзначити, що навіть при такому грубому наближенні рішення для швидкості виявляється на кілька порядків точніше. Так, з удвічі меншим інтервалом (0,025 с) в момент часу 0,5 с (тобто також після 20 кроків) значення, отримане чисельними методами, виходить рівним 0,943383 м / c , а точне значення - 0,953384 м / c і відрізняється всього на 0,000001 м / с або близько 0,0001%.

Причина цієї розбіжності в точності для координати і швидкості полягає в тому, що задача Коші

у разі швидкості зводилася до рівняння



тобто похідна швидкості не залежала явно від часу, а залежала тільки від значення самої швидкості. У той же час для координат ми мали диференціальне рівняння

тобто в правій частині ми, чисельно розв'язуючи задачу, могли отримати лише залежність від часу. Явної залежності від координати тут не було.

Таким чином, різниця в точності рішення обох рівнянь пояснюється відмінністю їх правих частин, а саме - від яких параметрів ці праві частини залежать.

в'язкість фізичний рунге чисельний



Висновки

У цій роботі було вивчено рух твердого тіла у в'язкій рідини, складено рівняння руху і знайдені точні вирази для швидкості і координат тіла. Був також вивчений і застосований для чисельного рішення диференціального рівняння метод Рунге-Кутти.

Чисельний метод застосовувався для вирішення завдання в два етапи - спочатку було вирішено диференціальне рівняння для швидкості (при цьому в рівнянні похідна явно залежала тільки від швидкості, і не залежала від часу), потім було записано диференціальне рівняння, що зв'язує поняття координати і швидкості , знаючи значення швидкості в потрібні моменти часу, знову вирішено чисельними методами (у цьому випадку похідна виходила функцією тільки від часу, і явно не залежала від координати).

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

Рішення ж для координати при малому числі кроків відрізнялися від точних рішень на значення порядку декількох відсотків, однак зі збільшенням числа кроків похибка помітно знизилася - аж до сотих часток відсотка при 10.000 кроків.

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



Список використаних джерел

1. О.Г. Ревинская, Н.С. Кравченко; Томский политехнический университет. Движение тела в вязкой среде -- Томск: Изд-во Томского политехнического университета, 2011 г.

2. Савельев И.В. Курс общей физики - М.:Наука - 2008 г.

3. Мышенков В.И., Мышенков Е.В. Численные методы. - 2009 г.

4. Керниган Б., Ритчи Д. Язык программирования Си -- 2-е изд. -- М.: Вильямс, 2007. -- С. 304.

5. Мак-Кракен Д.Численные методы и программирование. - 1990 г.



Додаток А

Текст програми

#include <windows.h>

#include <stdio.h>

#include <stdlib.h>

#include <math.h>

#define g 9.80665

double ro_fluid = 10, ro_body = 100;

double teta = 100;

double R = 2;

double f(double x, double V)

{

return g * (1.0 - ro_fluid / ro_body) - 4.5 * teta * V / (ro_body * R * R);

}

double* RungeCuttaV(double t0, double tend, double h, double Y0, double V0) {

double tn, Vn = V0;

// Функція f з диф. рівняння

// Знаходження швидкості методом Рунге-Кутти

double tau = ro_body * R * R / (4.5 * teta

double U = (1 - ro_fluid / ro_body) * g * tau;

int N = (tend - t0) / h;

double* VArr = (double*)malloc(N * sizeof(double));

int i = 0;

// Коефіціенти з диф. рівняння

// Створюємо масив

// Лічільник елементів

for (tn = t0; tn < tend; tn += h)

{

double k1, k2, k3, k4

double Vn1;

k1 = f(tn, Vn);

k2 = f(tn + h / 2, Vn + k1 * h / 2);

k3 = f(tn + h / 2, Vn + k2 * h / 2);

k4 = f(tn + h, Vn + k3 * h);

Vn1 = Vn + h * (k1 + 2 * k2 + 2 * k3 + k4) / 6.0;

// Коефіцієнти метода Рунге-Кутти

printf("V(runge)(%lf) = %lf, V(real)(%lf) = %lf\n", tn, Vn, tn, U * (1 - exp (- tn / tau)));

Виводимо значення швидкості, знайдені методом Рунге та із точної формули.

VArr[i++] = Vn1;

Vn = Vn1;

}

return VArr;

}

void RungeCuttaY(double t0, double tend, double h, double Y0, double V0)

{ // Записуємо значення у масив.

// Знаходимо координату методом Рунге із рівняння dY/dt = V, де V - це значення знайдені функцією RungeCuttaV

double* VArr = RungeCuttaV(t0, tend, 0.5 * h, Y0, V0);

int i = 0;

double tn, Yn = Y0;

double tau = ro_body * R * R / (4.5 * teta);

double U = (1 - ro_fluid / ro_body) * g * tau;

printf("\n time | Y (runge) | Y (pres.form)\n"

"-----------+-------------+--------------\n");

for (tn = t0; tn < tend; tn += h, i++)

{

double k1, k2, k3, k4;

double Yn1;

// У методі потрібно використати значення у момент часу кратні не лише інтервалу, та і його половині, тому треба знайти значення швидкості у два рази більшим числі точок (крок робимо у два рази меншим)

k1 = VArr[i * 2]

k2 = k3 = VArr[i * 2 + 1]

k4 = VArr[i * 2 + 2];

Yn1 = Yn + h * (k1 + 2 * k2 + 2 * k3 + k4) / 6.0; // VArr - значення швидкості у потрібні моменти часу.

Парні елементи - у моменти, кратні інтервалу.

printf("%9.5lf | %10.6lf | %10.6lf \n", tn, Yn, U * (tn - tau * (1 - exp (-tn / tau))));

Yn = Yn1;

}

}

// Виводимо значення координат, знайдених методом Рунге-Кутти та значення із точної формули.

void ReadData()

{

printf("Plotnost jitkosti (kg/m^3): ");

scanf("%lf", &ro_fluid);

printf("Plotnost tela (kg/m^3): ");

scanf("%lf", &ro_body);

printf("Koeficient vyazkosti (Pa*c): ");

scanf("%lf", &teta);

printf("Radius shara (m): ");

scanf("%lf", &R);

} // Фізичні параметри.

// Щільність рідини

// Щільність тіла

// Коефіцієнт в'язкості

// Радіус кулі

int main()

{

double t, dt;

ReadData();

printf("Enter vremennoi interval: ");

scanf("%lf", &t);

printf("Enter shag: ");

scanf("%lf", &dt);

RungeCuttaY(0, t, dt, 0, 0);

getch();

return 0;

}

// Кінцевий момент часу

//Часовий інтервал(крок) між двома послідовними моментами



Додаток Б

Блок схеми програми

Рис. 4 Блок схема функції Read Data()



Рис. 5 Блок схеми функцій RungeCuttaV та RungeCuttaY

Размещено на Allbest.ur