Simulace - 8. díl

Osmý díl seriálu se zabývá nejjednodušší simulací z mechaniky – simulací pružiny.

V minulých dílech jsme se zabývali výlučně elektrotechnickými záležitostmi. Není na škodu si udělat malou exkurzi do mechaniky a nasimulovat např. chování pružiny.

Opět teorie

Jak jsem psal v jednom z minulých dílů, základem simulace je vždy teorie. Vycházíme z matematických opisů skutečných jevů, jejichž chování v simulaci odpovídá skutečnosti. Naší soustavou skutečných jevů bude břemeno zavěšené na pružině, k čemuž budeme hledat matematický opis.

Na naše břemeno, kromě síly pružiny, působí i gravitační síla vlivem vlastní hmotnosti umístěné v gravitačním poli Země (na Měsíci by se soustava chovala jinak). Jelikož se jedná pouze o exkurzi, bude nám to prozatím stačit a zapomeneme na jevy v materiálu pružiny, jako i na její vlastní hmotnost, a na další drobné jevy. Bude nám stačit tvrzení známé z mechaniky, že síla pružiny je úměrná jejímu stlačení nebo roztažení (deformaci):

kde „k“ je parametrem pružiny, nazývaný též tuhost pružiny. Proměnná „x“ je délková deformace pružiny buď ve směru tažném nebo tlačném, přičemž výchozím bodem (x = 0) je klidový stav pružiny, kdy nepůsobí žádné síly (ani vlastní hmotnost pružiny). Abychom se alespoň trochu přiblížili reálnému světu, budeme uvažovat tlumení pružiny, které zabrání tomu, abychom získali harmonický oscilátor. Tlumení je způsobeno chováním materiálu při pružné deformaci a je závislé na její rychlosti, tedy:

Gravitační síla je vyvolána polem Země, jehož projevem je gravitační zrychlení, často označovaném „g“, v dané zeměpisné šířce konstantní (9,81 m/sˆ2). Síla vyvolaná gravitačním zrychlením je podle Newtonova zákona:

kde „m“ je hmotnost břemene v kilogramech. A protože Newton byl chytrý pán, použijeme další jeho zákon, o rovnováze sil:

Získali jsme tedy výslednici zúčastněných sil. Z této výslednice bychom měli získat polohu břemene a tedy i deformaci pružiny, ze které budeme odvozovat sílu pružiny. Nezbývá nám, než se opět obrátit na Newtonův zákon:

Poslední rovnice je úpravou do tvaru lineární diferenciální rovnice druhého řádu s konstantní pravou stranou. Řečeno srozumitelnější formou: jsme opět u podobného problému jako v předchozích dvou dílech. Danou rovnici si tedy upravíme na dvě diferenciální rovnice prvního řádu a poté i do tvaru pro numerické řešení, na diferenční rovnice. Níže uvedená substituce dx/dt = v není náhodná, protože popisuje definici rychlosti. Můžeme si tedy vykreslit i průběh rychlosti pohybu tělěsa v čase.

Rutina…

Poslední dvě rovnice stačí jen naprogramovat, uvést počáteční podmínky a parametry systému a můžeme „pružit“. Ještě je nutné poznamenat, že tak jak jsme zavedli konvenci směru působení sil, tak musíme počítat s určováním parametrů systému. V našem konkrétním případě z fyziky víme, že síla pružiny působí proti gravitační síle břemene. Tento fakt musíme zahrnout i do rovnic a koeficient „k“ pružiny (tuhost) definovat se záporným znaménkem, protože jsme v rovnici definovali směr všech sil shodný. Stejným způsobem definujeme i koeficient tlumení B, který působí ve stejném směru jako síla pružiny, tedy záporný.

Výpis programu považuji za zbytečný, protože se od toho v minulém díle prakticky neliší, ale přesto ho zde uvedu. Parametry jsem volil odhadem.

// Simulace tělesa zavěšeného na pružině
// www.programujte.com 2006 - Matěj Pácha
// Visual C++ 2005 Express Edition
// Výsledek ukládá do souboru "sim.txt"

#include "stdafx.h"
#include <fstream>
#include <iostream>

using namespace std;

int _tmain(int argc, _TCHAR* argv[])
{
	ofstream fout;
	fout.open("sim.txt");

	double t,h,x,v,g,m,k,b,konec;
	
	t=0;			//proměnná času
	x=0;			//nastavíme počáteční podmínky
	v=0;			//
	g=9.81;			//gravitační zrychlení
	h=0.01;			//krok, se kterým počítáme
	konec=10;		//čas konce simulace
	m=1;			//Hmotnost tělesa m = 1 kg
	k=-100;			//Tuhost pružiny k = -100
	b=-10;			//Tlumení pružiny B = -10
	
	for(t = 0; t < konec ; t += h)
	{
					//simulační rovnice
		v += h*(g + k/m*x + b*v/m);
		x += h*v;
					//výstup do souboru
		fout << t << "\t" << x << "\t" << v << "\n";
	}
	fout.close();
	return 0;
}

Interpretace polohy x v prvním grafu je kladná (v kladném směru), tj. ve směru gravitační síly (dolů). Druhý graf znázorňuje průběh rychlosti v čase.

Pokud by se někomu nelíbila představa záporných koeficientů k a B, není nic jednoduššího, než změnit znaménka příslušných sil v rovnici. Vše je pouze ve zvolené konvenci, ale je nutné tuto konvenci dodržovat (podobně jako jazyk a význam jeho slov).

Závěr

Při pohledu na tento díl je opět patrné, že většinu času zabírá teoretický rozbor dané problematiky. Se základním matematickým aparátem a znalostí programovacího jazyka je sestavení simulačního programu otázkou několika minut. Každopádně je ale nutné si uvědomit všechny souvislosti a výsledek simulace správně vyhodnotit. Nemluvím zde jen o různých případech divergence výsledku, zaokrouhlovacích chybách a chybách z volby kroku. V simulaci i tak dost jednoduché soustavy jsme zanedbali spoustu vlivů, které se v reálném světě vyskytují. Přesný simulační model, i kdyby existoval, by byl značně složitý a na numerické řešení značně nestabilní, což je i důvod, proč zjednodušení zavádíme. Musíme ale nést tíhu důsledků zjednodušení a brát výsledky simulace s přiměřenou rezervou.