Lekce 40 - Fyzikální simulace lana

Přichází druhá část dvoudílné série o fyzikálních simulacích. Základy už známe, a proto se pustíme do komplikovanějšího úkolu - klávesnicí ovládat pohyby simulovaného lana. Zatáhneme-li za horní konec, prostřední část se rozhoupe a spodek se vláčí po zemi. Skvělý efekt.

Celé demo založíme na jednoduchém fyzikálním enginu z lekce 39. Nyní už byste měli umět aplikovat libovolné síly na jakýkoli hmotný objekt, přepočítat jeho novou pozici i rychlost a samozřejmě provádět operace s 3D vektory. Pokud něčemu z toho nerozumíte, vraťte se k minulé lekci, popř. zkuste jiné zdroje.

Předpokladem pro fyzikální simulace je implementace fyzikálních podmínek a závislostí, kdy se snažíme o to, aby vše vypadalo jako v reálném prostředí. Nejvíce na očích bývá vždy dynamika - pohybové reakce objektů na uživatelovy příkazy. Právě na nich hodnotí, zda se naše práce podařila, či ne. Na úplném začátku se vždy musí najít vhodný kompromis mezi rychlostí a kvalitou. Jsme schopni všimnout si atomů, elektronů nebo protonů? Ne. Určitě bude stačit aproximace pohybu skupiny částic.

Matematika pohybů

Klasická mechanika reprezentuje předměty jako částice v prostoru, které mají určitou hmotnost. Jejich zrychlení závisí na působících silách a pozice na uplynulém čase od minulých výpočtů. Můžeme ji použít k simulaci chování objektů, které jsou viditelné prostým okem (pro mikrosvět platí jiné fyzikální zákony). V lekci 39 jsme s její pomocí implementovali pohyb v gravitaci a objekt zavěšený na pružině. Nyní zkusíme simulovat složitější předmět - lano.

Výkon počítače, který používáme k simulaci

Rychlost počítače je hlavním omezením pro množství detailů simulace. Například při simulaci chodícího člověka eliminujeme na pomalém počítači pohyb prstů na nohou, které jsou sice důležité, ale výsledek vypadá přesvědčivě i bez nich. Musíme je vynechat z jednoduchého důvodu: počítač by nestíhal provádět potřebné výpočty, kterých je i bez nich až příliš.

Jako minimální požadavek pro simulaci určíme počítač s frekvencí procesoru kolem 500 MHz. Z toho plyne omezení počtu detailů. Při implementaci použijeme knihovnu Physics1.h z lekce 39. Tato knihovna obsahuje třídu Mass (hmota), která reprezentuje jeden hmotný bod. Spojením několika za sebe získáme fyzikální model, který reprezentuje lano. Můžeme usoudit, že se bude kývat a různě vlát, ale nebude moci kroužit, protože kroužení nejde pomocí hmotných bodů implementovat (nemohou rotovat okolo os). Určíme si, že body spojíme po vzdálenostech 10 cm. Tato hodnota vychází z toho, že kvůli rychlosti použijeme maximálně 50 až 100 částic na 3 až 4 metrové lano. Z toho plynou 3 až 8 cm velké mezery, což bude ještě větší přesnost, než jsme původně zamýšleli.

Odvození rovnice pohybu

Rovnice pohybu je v matematice vyjádřena diferenciální rovnicí druhého stupně. V modelu lana si můžeme každé dvě sousední částice, ze kterých je složeno, představit jako konce jedné pružiny. Znak o představuje částici a pomlčka pružinu.

o----o----o----o

První částice poutá druhou, ta zase třetí, ta čtvrtou atd. Vzniká jakýsi řetězec, složený ze čtyř částic a tří pružin. Pružiny představují zdroje síly mezi každými dvěma částicemi. Zapamatujte si, že síla pružiny se dá vyjádřit takto:

síla = -k * x

k: konstanta určující tuhost pružiny

x: vzdálenost mezi body pružiny

Kdybychom použili tuto rovnici, lano by se za chvíli smrštilo, protože z rovnice vyplývá, že dokud není vzdálenost částic nulová, působí na ně síla. Všechny by tíhly k ostatním a to nechceme. Představte si lano položené na stole. Chceme, aby to naše mělo stejnou pevnost a tudíž musíme explicitně udržovat jeho délku konstantní. Abychom toho dosáhli, musí být při určité kladné vzdálenosti síla pružiny nulová. Nic těžkého:

síla = -k * (x - d)

k: konstanta určující tuhost pružiny

x: vzdálenost mezi body pružiny

d: konstanta označující kladnou vzdálenost částic, při které pružina zůstane ve stálé poloze

Z rovnice vyplývá, že pokud se bude vzdálenost mezi částicemi rovnat konstantě d, nebudou aplikovány žádné síly. Definovali jsme si lano složené ze sta částic. Zvolíme-li d = 5 cm (0,05 metrů), získáme pevné pětimetrové lano. Pokud bude x větší než d, pružina se začne natahovat a při menším x naopak smršťovat. Každopádně se bude neustále nacházet v blízkosti bodu rovnováhy.

Máme zajištěn celkem slušný pohyb, ale něco mu schází. A to něco jsou ztráty - napětí vláken, jejich tření a podobně. Beze ztrát si fyzikální systém uchovává veškerou energii, kterou mu dodáme - lano se nikdy nepřestane houpat. Než se začneme ztrátám věnovat, pojďme se nejprve podívat na kód.

Třída pružiny

Třída pružiny (anglicky spring) popisuje dvě částice a silové působení pružiny na každou z nich.

class Spring// Třída pružiny

{

public:

Mass* mass1;// Částice na prvním konci pružiny

Mass* mass2;// Částice na druhém konci pružiny

float springConstant;// Konstanta tuhosti pružiny

float springLength;// Délka, při které nepůsobí žádné síly

float frictionConstant;// Konstanta vnitřního tření

V konstruktoru nastavíme vnitřní datové členy na hodnoty, které byly předány v parametrech.

Spring(Mass* mass1, Mass* mass2, float springConstant, float springLength, float frictionConstant)// Konstruktor

{

// Nastavení členských proměnných

this->springConstant = springConstant;

this->springLength = springLength;

this->frictionConstant = frictionConstant;

this->mass1 = mass1;

this->mass2 = mass2;

}

Nejdůležitější částí třídy je metoda solve(), ve které se aplikují síly. Za číslo x se má dosadit vzdálenost mezi okrajovými body, v našem případě se jedná o délku 3D vektoru, kterou vypočteme odečtením pozic bodů.

void solve()// Aplikování sil na částice

{

Vector3D springVector = mass1->pos - mass2->pos;// Vektor mezi částicemi

float r = springVector.length();// Vzdálenost částic

Vytvoříme další vektor, který bude označovat výslednou sílu. Konstruktor ji automaticky vynuluje. Protože v další části dělíme, musíme ošetřit, jestli se číslo r nerovná nule. Pokud je vše v pořádku, můžeme přistoupit k výpočtu.

Vector3D force;// Pomocný vektor síly

if (r != 0)// Proti dělení nulou

{

Chceme-li dosáhnout rovnice uvedené výše, potřebujeme získat jednotkový vektor, který reprezentuje směr působení síly. Defakto ho už máme uložen v objektu springVector, ale je v něm navíc započítána i jeho délka a to nechceme. Dá se však velice jednoduše odstranit dělením (springVector / r). Dále se pokusíme implementovat část (x - d). Máme jak vzdálenost bodů, tak délku pružiny, nic nám proto nebrání, abychom je odečetli (r - springLength). Konečný výsledek ještě vynásobíme tuhostí pružiny. Záporná hodnota označuje, že se lano bude spíše vléci než odrážet.

force += (springVector / r) * (r - springLength) * (-springConstant);// Výpočet síly

}

Vyřešili jsme silové působení pružiny, ale ještě nám chybí ztráty energie v materiálu. Pokud se na hmotu aplikuje síla v opačném směru, než se pohybuje, zpomalí. Kam se ztratila pohybová energie? Mohla se například přeměnit ve tření a následně v tepelnou energii.

třecí síla = -k * rychlost

k: konstanta představující velikost ztrát (např. v závislosti na drsnosti povrchu)

rychlost: rychlost hmoty, na kterou působí třecí síla

Tato rovnice by se dala napsat i jinak (více složitě), ale nám bude tato verze bohatě stačit. Všimněte si, že lze dosadit pouze rychlost jednoho bodu, ale pružina se skládá ze dvou. Co dělat? Vypočteme rozdíl rychlostí a předáme ho jako relativní rychlost. Ztráty tedy budou představovat vnitřní tření v materiálu.

force += -(mass1->vel - mass2->vel) * frictionConstant;// Zmenšení síly o tření

Podle Newtonova zákona akce a reakce aplikujeme na jeden bod pružiny kladnou sílu a na druhý zápornou. (Působí-li jedno těleso na druhé silou, působí i druhé na první. Obě síly jsou stejně velké, ale opačně orientované.) Představte si dvě loďky na jezeře. Po odstrčení se nezačne pohybovat jenom jedna, ale obě. Pokud by bylo jedno těleso mnohonásobně těžší než to druhé, mohlo by se silové působení na něj zanedbat, jeho zrychlení by se blížilo k nule. Představte si například raketu v gravitačním poli planety, která je přitahována dolů do jejího středu. Raketa se zároveň snaží přitáhnout planetu nahoru, nicméně nemá nejmenší šanci :-) a tak se druhé působení jednoduše zanedbává. To ale není náš případ, protože oba naše objekty mají stejnou hmotnost.

mass1->applyForce(force);// Aplikování síly na částici 1

mass2->applyForce(-force);// Aplikování opačné síly na částici 2

}

};

Nyní už máme vyřešenu rovnici pohybu, kterou si můžeme představit jako silové působení pružin. Abychom simulaci dokompletovali, přidáme ještě gravitační sílu, tření lana se vzduchem a plochý povrch, po kterém lanem posunujeme. První dvě jsou jednoduché, nejprve odpor vzduchu:

odpor vzduchu = -k * rychlost

k: konstanta, která představuje velikost odporu

rychlost: rychlost pohybu

A gravitace...

gravitační síla = gravitační zrychlení * hmotnost

Gravitace i odpor vzduchu působí na každou částici lana zvlášť. Co se zemí? Můžeme si vytvořit pomyslnou rovinu a testovat, kolize s částicemi lana. Pokud se ocitne pod úrovní roviny, vyvýšíme ji a rozšíříme působící sílu o tření s podlahou.

Nastavení počátečních hodnot simulace

V tuto chvíli je prostředí připraveno pro simulaci, ale potřebujeme definovat jednotky používaných fyzikálních veličin. Vzdálenost bude specifikována v metrech, čas v sekundách a hmotnost v kilogramech. Gravitaci nastavíme tak, aby působila ve směru záporné části osy y se zrychlením 9,81 m*s-2 (odpovídá gravitaci na zemi). Před spuštěním umístíme lano rovnoběžně se zemí ve vzdálenosti čtyř metrů od ní. Aby se jí mohlo dotknout bude mít délku (v klidovém stavu) také čtyři metry, Z toho vyplývá 5 cm vzdálenost mezi jednotlivými částicemi (4 m / 80 částic = 0,05 m = 5 cm). Normální délku pružiny mezi částicemi (délka bez působení žádných sil) tedy nastavíme na těchto 5 cm, aby lano na začátku simulace nebylo ani napnuté ani prohnuté. Celkovou hmotnost lana určíme na 4 kg a to dává 0,05 kg (= 50 gramů) na každou částici. Pro přehlednost si vše shrneme:

Nyní zkusíme vypočítat konstantu určující tuhost pružiny. Pokud budeme držet lano za jeden konec, tak se natáhne. Představte si elastické lano se závažím. Je to úplně stejné, akorát my nemáme pouze jednu pružinu ale hned několik. Nejvíce se natáhne horní pružina, protože drží hmotnost všech ostatních částic - prakticky celé lano. Naopak nejméně se prodlouží ta spodní, protože je na ní zavěšena jenom jedna jediná částice. Nechceme, aby se ta horní natáhla více než o 1 cm.

f = hmotnost lana * gravitační zrychlení = (4 * 9,81) N ~= 40 N

Síla pružiny odpovídá přibližně 40 N. Dá se ale vyjádřit i jinak:

síla pružiny = -k * x = -k * 0,1 metrů

Suma těchto sil by se měla rovnat nule.

40 N + (-k * 0,01 metrů) = 0

Vypočteme a získáme

k = 4000 N/m

Pro snadnější zapamatovatelnost budeme předpokládat 10000 N/m, což nám dává větší tuhost lana, které se v horní části natáhne jen o 4 mm.

Aby se nalezla konstanta vnitřního tření v laně, museli bychom podniknout mnohem více komplikovanější výpočet než je ten výše. Proto jsem použil metodu pokusů a omylů a získal...

konstanta vnitřního tření = 0,2 N/(m/s)

... což vypadá celkem realisticky.

Třída simulace lana

Předtím než začneme zkoumat tření se vzduchem a se zemí, pojďme se podívat na třídu simulace lana, která je odvozená od obecné třídy simulace z lekce 39. Tato třída obsahuje čtyři metody potřebné pro běh simulace.

V potomku základní třídy přepíšeme funkci solve() a funkci simulate(float dt), protože kvůli lanu potřebujeme jejich speciální implementaci. Solve() slouží k aplikování sil a simulate(float dt) k ovládání lana za jeden konec pověšený v prostoru. Jak už bylo řečeno, třída RopeSimulation je potomkem třídy Simulation. Obecnou simulaci rozšiřuje o lano složené z hmotných bodů (částic) pospojovaných pružinami. Tyto pružiny mají vnitřní tření a klidovou délku. Jeden konec lana je udržován v prostoru na souřadnicích ropeConnectionPos a je jím možno pohybovat pomocí metody setRopeConnectionVel(Vector3D ropeConnectionVel). Třída dále zajišťuje tření se vzduchem a s rovinným povrchem (nebo-li se zemí), jehož normála z něj vychází ve směru kladné části osy y.

class RopeSimulation : public Simulation// Třída simulace lana

{

public:

Spring** springs;// Pružiny spojující částice

Vector3D gravitation;// Gravitační zrychlení

Vector3D ropeConnectionPos;// Bod v prostoru; pozice první částice pro ovládání lanem

Vector3D ropeConnectionVel;// Rychlost a směr požadovaného pohybu

float groundRepulsionConstant;// Velikost odrážení částic od země

float groundFrictionConstant;// Velikost tření částic se zemí

float groundAbsorptionConstant;// Velikost absorpce sil částic zemí (vertikální kolize)

float groundHeight;// Pozice roviny země na ose y

float airFrictionConstant;// Konstanta odporu vzduchu na částice

Konstruktorem inicializujeme všechny členské proměnné třídy, alokujeme paměť pro všechny potřebné pružiny a umístíme je v řadě rovnoběžně se zemí.

RopeSimulation(// Konstruktor třídy

int numOfMasses,// Počet částic

float m,// Hmotnost každé částice

float springConstant,// Tuhost pružiny

float springLength,// Délka pružiny v klidovém stavu

float springFrictionConstant,// Konstanta vnitřního tření pružiny

Vector3D gravitation,// Gravitační zrychlení

float airFrictionConstant,// Odpor vzduchu

float groundRepulsionConstant,// Odrážení částic zemí

float groundFrictionConstant,// Tření částic se zemí

float groundAbsorptionConstant,// Absorpce sil zemí

float groundHeight// Pozice země na ose y

) : Simulation(numOfMasses, m)// Inicializace předka třídy

{

this->gravitation = gravitation;

this->airFrictionConstant = airFrictionConstant;

this->groundFrictionConstant = groundFrictionConstant;

this->groundRepulsionConstant = groundRepulsionConstant;

this->groundAbsorptionConstant = groundAbsorptionConstant;

this->groundHeight = groundHeight;

for (int a = 0; a < numOfMasses; ++a)// Nastavení počáteční pozice částic

{

masses[a]->pos.x = a * springLength;// Offsety jednotlivých částic

masses[a]->pos.y = 0;// Rovnoběžně se zemí

masses[a]->pos.z = 0;// Rovnoběžně s obrazovkou

}

springs = new Spring*[numOfMasses - 1];// Alokace paměti pro ukazatele na pružiny

for (a = 0; a < numOfMasses - 1; ++a)// Vytvoření jednotlivých pružin

{

springs[a] = new Spring(masses[a], masses[a + 1], springConstant, springLength, springFrictionConstant);// Dvě částice na pružinu

}

}

Jak jsme si už řekli výše, funkce solve() slouží k aplikování sil na všechny objekty, ze kterých je lano složeno.

void solve()// Aplikování sil

{

Nejdříve ošetříme všechny pružiny; na jejich pořadí nezáleží. Třída obsahuje svoji vlastní funkci.

for (int a = 0; a < numOfMasses - 1; ++a)// Prochází pružiny

{

springs[a]->solve();// Aplikování sil na pružinu

}

V cyklu přes všechny částice zajistíme působení gravitace a odpor vzduchu.

for (a = 0; a < numOfMasses; ++a)// Prochází částice

{

masses[a]->applyForce(gravitation * masses[a]->m);// Gravitace

masses[a]->applyForce(-masses[a]->vel * airFrictionConstant);// Odpor vzduchu

Síly ze země vypadají trochu komplikovaněji, ale jsou právě tak jednoduché, jako všechny ostatní. Země na částici může působit pouze tehdy, pokud se navzájem dotknou, tj. částice se na ose y nachází pod její úrovní.

if (masses[a]->pos.y < groundHeight)// Kolize se zemí,

{

Smýkání lana na zemi je zajištěno třecí silou, která ze své podstaty zanedbává rychlost na ose y. Y je směr, kterým země (její normálový vektor) směřuje vzhůru; smýkání nemůže působit v tomto směru.

Vector3D v;// Pomocný vektor

v = masses[a]->vel;// Grabování rychlosti

v.y = 0;// Vynechání rychlosti na ose y

masses[a]->applyForce(-v * groundFrictionConstant);// Třecí síla země

Opakem smýkání je absorpční efekt, kdy se síla aplikuje pouze ve směru, kterým země směřuje vzhůru. Proto obě ostatní složky vynulujeme. Absorpce nemůže na částici působit tehdy, když se vzdaluje od země. Pokud bychom nepřidali podmínku v.y < 0, lano by tíhlo k zemi, i když by se jí už nedotýkalo.

v = masses[a]->vel;// Grabování rychlosti

v.x = 0;// Zanedbání rychlosti na osách x a z

v.z = 0;

if (v.y < 0)// Pouze při kolizi směrem k zemi

{

masses[a]->applyForce(-v * groundAbsorptionConstant);//Absorpční síla

}

Síla odrazu je poslední ze sil, kterou vyvolává kolize se zemí. Země odráží částice právě tak, jako kdyby se mezi nimi nacházela pružina. Její síla je přímo úměrná rychlosti částice při nárazu.

Vector3D force = Vector3D(0, groundRepulsionConstant, 0) * (groundHeight - masses[a]->pos.y);// Síla odrazu

masses[a]->applyForce(force);// Aplikování síly odrazu

}

}

}

Abychom vyvolali dojem média, které za jeden konec drží lano a pohybuje s ním, musíme přepsat metodu simulate(float dt).

void simulate(float dt)// Simulace lana

{

Nejdříve ze všeho zavoláme metodu předka, potom k aktuální pozici přičteme rychlost pohybu a nakonec ošetříme náraz do země.

Simulation::simulate(dt);// Metoda předka

ropeConnectionPos += ropeConnectionVel * dt;// Zvětšení pozice o rychlost

if (ropeConnectionPos.y < groundHeight)// Dostala se částice pod zem?

{

ropeConnectionPos.y = groundHeight;// Přesunutí na úroveň země

ropeConnectionVel.y = 0;// Nulování rychlosti na ose y

}

Pomocí právě získaných parametrů nastavíme vlastnosti částice na indexu nula.

masses[0]->pos = ropeConnectionPos;// Pozice první částice

masses[0]->vel = ropeConnectionVel;// Rychlost první částice

}

Potřebujeme funkci, pomocí které budeme moci nastavit rychlost první částice.

void setRopeConnectionVel(Vector3D ropeConnectionVel)// Nastavení rychlosti první částice

{

this->ropeConnectionVel = ropeConnectionVel;// Přiřazení rychlostí

}

Tím končí vysvětlování vnitřních závislostí třídy RopeSimulation. Její objekt v aplikaci vytváříme dynamicky pomocí operátoru new. Změnou hodnot předávaných konstruktoru můžete docílit téměř libovolného chování lana. Všimněte si, že zemi umisťujeme -1.5f jednotek pod osou y. Lano inicializujeme na nule rovnoběžně se zemí. To nám tedy dává možnost vidět hned na začátku efektní pád a kolizi se zemí.

RopeSimulation* ropeSimulation = new RopeSimulation(// Vytvoření objektu simulace lana

80,// 80 částic

0.05f,// Každá částice váží 50 gramů

10000.0f,// Tuhost pružin

0.05f,// Délka pružin, při nepůsobení žádné sily

0.2f,// Konstanta vnitřního tření pružiny

Vector3D(0, -9.81f, 0), // Gravitační zrychlení

0.02f,// Odpor vzduchu

100.0f,// Síla odrazu od země

0.2f,// Třecí síla země

2.0f,// Absorpční síla země

-1.5f);// Poloha země na ose y

Stejně jako v lekci 39 existuje maximální možná hodnota dt simulace. S výše uvedenými parametry konstruktoru činí přibližně 0,002 sekund. Pokud vaše změna tuto hodnotu sníží, může simulace vypadat poněkud nestabilně a lano nemusí pracovat správně. Abyste poměry stabilizovali, musíte najít nové maximální možné dt. Velké síly a/nebo malé hmotnosti znamenají větší nestabilitu, protože zrychlení bude vyšší (zrychlení = síla / hmotnost).

// Funkce Update(DWORD milliseconds)

float dt = milliseconds / 1000.0f;// Převod milisekund na sekundy

float maxPossible_dt = 0.002f;// Maximální možné dt

int numOfIterations = (int)(dt / maxPossible_dt) + 1;// Výpočet počtu opakování při této aktualizaci

if (numOfIterations != 0)// Proti dělení nulou

{

dt = dt / numOfIterations;// Aktualizace dt podle numOfIterations

}

for (int a = 0; a < numOfIterations; ++a)// Opakování simulace

{

ropeSimulation->operate(dt);

}

}

Lanem můžete pohybovat pomocí šipek a kláves HOME a END. Pohrajte si, stojí to za to. Simulační procedura hodně závisí na výkonu procesoru, tudíž jsou také důležité optimalizace kompilátoru. Při standardním Visual C++ Release nastavení běží program více než 10 krát rychleji než v Debug módu, pro který činí minimální frekvence procesoru cca. 500 MHz.

V tomto tutoriálu je představeno kompletní fyzikální nastavení, teoretická funkce, design a implementace. Více pokročilejší simulace uvnitř vypadají úplně stejně jako tato.

napsal: Erkin Tunca <erkintunca (zavináč) icqmail.com>
přeložil: Michal Turek - Woq <WOQ (zavináč) seznam.cz>

Zdrojové kódy

Lekce 40

<<< Lekce 39 | Lekce 41 >>>