\newpage
\section[Metodika]{Metodika statické a dynamické analýzy}
\label{sec:metodika}

Statická a dynamická analýza je provedena numericky metodou konečných prvků. 
MKP model je vytvořen podle zadaných podkladů v~programu ANSYS.

\subsection{Analýza metodou konečných prvků}
	
Metoda konečných prvků (MKP) je numerická metoda pro řešení spojitých
prostředí (kontinua). Princip spočívá v~diskretizaci kontinua do konečného počtu
prvků, přičemž zjišťované parametry jsou určovány v~jednotlivých uzlech.

Diskretizací spojitého modelu se získá konečný počet prvků, které jsou
definovány uzly, materiálovými a~geometrickými charakteristikami.
Každý uzel může mít v~obecném případě až 6 stupňů volnosti
(3 posuny a 3 rotace). Zároveň lze pro každý prvek určit matici
tuhosti, která definuje závislost deformací na vnějším zatížení. Pomocí
okrajových podmínek je možné řešit soustavu rovnic a~získat tak pole
deformací, ze kterého se vypočítá pole napětí.

Sestavení matic tuhosti je založeno na Lagrangeově principu virtuálních
prací: \uv{Těleso je v~rovnováze, jestliže celková potenciální energie
deformace soustavy je minimální.}

Základní rovnice je:
\begin{equation}
	 \mathbf R =  \mathbf K\cdot \mathbf r,
\end{equation}
\begin{flalign*}
 \text{kde} \quad \mathbf R & \quad \text{\ldots vektor vnějšího zatížení} & \notag \\
 \mathbf K & \quad \text{\ldots matice tuhosti MKP modelu} \notag \\
 \mathbf r & \quad \text{\ldots vektor deformací.} \notag
\end{flalign*}

\nomenclature{{\bfseries R}}{Vektor vnějšího zatížení}
\nomenclature{{\bfseries K}}{Matice tuhosti MKP modelu}
\nomenclature{{\bfseries r}}{Vektor deformací}

Podrobnější informace lze získat v~publikaci {\em A~First Course in the Finite
Element Method}~\cite{FEM}.

\subsection{Modelování proudění vzduchu}
\label{sec:metodika:cfd}

Odvození Navier-Stokesových rovnic je podrobně uvedeno v \cite{cfdintro}. Zde
budou shrnuty pouze základní body.

\smallskip

Obecný model proudění tekutiny je založen na třech zákonech zachování:

\begin{itemize}
    \item {\bfseries Zákon zachování hmotnost} a z něho vyplívající vztah:
        Změna hmotnosti elementu je rovna celkovému přítoku/odtoku tekutiny do elementu.

    \item {\bfseries Zákon zachování hybnosti} a sním související druhý Newtonův pohybový 
        zákon: Časová změna hybnosti je rovna součtu sil působící na element tekutiny.

        \[ \dv{p}{t} = \va{F} \]

    \item {\bfseries První termodynamický zákon} -- celkové množství energie izolované
        soustavy je konstantní. Z pohledu změny energie je možné zákon přeformulovat:
        Změna vnitřní energie U je rovna součtu přidaného nebo odebraného tepla
        Q a vykonané práce W.

        \[ \dd{U} = \delta Q + \delta W \]
\end{itemize}

Zatímco první dva zákony platí pro každou částici, termodynamický zákon je
statistické povahy a platí pro \uv{dostatečný} počet molekul. Infinitezimální
element v mechanice tekutin je tedy nejmenší možný, pro který jeho makroskopické
vlastnosti nejsou ovlivněny jednotlivými molekulami. Rozměr takového elementu
v kartézském systému souřadnic značíme $\delta x$, $\delta y$ a $\delta z$.

\smallskip

Všechny vlastnosti tekutiny jsou funkcí času a polohy: hustota $\rho(t, x, y, z)$,
tlak $p(t, x, y, z)$, teplota $T(t, x, y, z)$ a rychlost $\va{u}(t, x, y, z)$.
Pro přehlednost se tato závislost implicitně předpokládá a zápis se zkracuje
na $\rho$, $p$, $T$ a $\va{u}$. Složky vektoru rychlosti se označují $u$, $v$,
$w$ pro směr $x$, $y$ a $z$.

\subsubsection{Rovnice kontinuity}

Zákona zachování hmotnosti v třírozměrném proudění můžeme vyjádřit pomocí rovnici kontinuity:

\begin{equation}\label{eq:metodika:cfd:kontinuum}
    \underbrace{\pdv{\rho}{t}}_{\text{Změna hmotnosti}} 
  + \underbrace{\nabla\cdot(\rho\va{u})\vphantom{\pdv{\rho}{t}}}_{\text{Přítok tekutiny}}
  = 0
\end{equation}

Rovnice \ref{eq:metodika:cfd:kontinuum} platí pro stlačitelné tekutiny. Pro
nestlačitelné se zjednoduší na:

\begin{equation}\label{eq:metodika:cfd:kontinuum}
    \nabla\cdot\va{u} = 0
\end{equation}

\subsubsection{Materiálová derivace}

Změnu určité fyzikální veličiny $\Phi$ spjaté s pohybem hmotné částice kontinua
můžeme v čase vyjádřit pomocí materiálové derivace:

\begin{equation}
    \label{eq:metodika:cfd:zmena_veliciny}
    \frac{D\Phi}{Dt} = \pdv{\Phi}{t} + \nabla\cdot(\rho\Phi\va{u})
\end{equation}

kde $\Phi$ je odpovídající veličina -- $u$, $v$, $w$ (momenty ve směru x, y a z)
a energie $E$.

\subsubsection{Momentové rovnice}

Podle druhého Newtonova pohybového zákona je změna momentu rovna součtu sil
působící na částici. Mezi tyto síly patří povrchové síly (tlakové a viskózní) a 
objemové síly (gravitační, Coriolisovi, atp.).

\smallskip

Vyjádříme-li stav napjatosti elementu tekutiny pomocí tlaku a devíti viskózních
komponent, můžeme celkovou sílu, působící na daný element, vyjádřit:

\begin{equation}
    \label{eq:metodika:cfd:zmena_moment}
    \begin{aligned}
    \pdv{(-p + \tau_{xx})}{x} &+ \pdv{\tau_{yx}}{y} + \pdv{\tau_{zx}}{z} \qq{ve směru x,}\\    
    \pdv{(-p + \tau_{xy})}{x} &+ \pdv{\tau_{yy}}{y} + \pdv{\tau_{zy}}{z} \qq{ve směru y,} \\
    \underbrace{\pdv{(-p + \tau_{xz)}}{x}}_{\text{normálové napětí}} &+ \underbrace{\pdv{\tau_{yz}}{y}
    + \pdv{\tau_{zz}}{z}}_{\text{smyková napětí}}  \qq{ve směru z.}
    \end{aligned}
\end{equation}

Změnu momentu pak můžeme zapsat:

\begin{equation}
    \label{eq:metodika:cfd:zmena_moment}
    \begin{aligned}
        \rho\frac{Du}{Dt} &= \pdv{(-p + \tau_{xx})}{x} &+& \pdv{\tau_{yx}}{y} &+& \pdv{\tau_{zx}}{z} &+ S_{Mx} \\
        \rho\frac{Dv}{Dt} &= \pdv{\tau_{xy}}{x} &+& \pdv{(-p + \tau_{yy})}{y} &+& \pdv{\tau_{zy}}{z} &+ S_{My} \\
        \rho\frac{Dw}{Dt} &= \pdv{\tau_{xz}}{x} &+& \pdv{\tau_{yz}}{y} &+& \pdv{(-p + \tau_{zz})}{z} &+ S_{Mz}    
    \end{aligned}
\end{equation}

kde $S_{Mx}$, $S_{My}$ a $S_{Mz}$ jsou účinky objemových sil.

\subsubsection{Termodynamické rovnice}

Obdobně je možné vyjádřit první termodynamický zákon pro element tekutiny
a získat tak energetickou rovnici.

\subsubsection{Stavová rovnice}

Pohyb tekutiny v prostoru je zatím popsán pěti parciálními diferenciálními
rovnicemi -- rovnicí kontinuity, třemi momentovými rovnicemi a energetickou
rovnicí, celkem o šesti neznámých ($p$, $\rho$, $E$, $u$, $v$ a $w$).

\smallskip

Celý systém doplňuje stavová rovnice, která určuje vztah mezi termodynamickými
veličinami $p$, $T$ a $\rho$:

\begin{equation}
    p = p(\rho, T)
\end{equation}

Pro ideální plyn má tvar:

\begin{equation}
    p = \rho RT
\end{equation}

Kapaliny a plyny se při nízkých rychlostech chovají jako nestlačitelné tekutiny
a proudění může být vyřešeno pouze pomocí rovnice kontinuity a momentových rovnic.

\subsubsection{Navier-Stokesovy rovnice}

Uvažujeme-li Newtonské tekutiny, pro které platí, že napětí je přímo úměrné
rychlosti deformace:

\begin{equation}
    \tau = \nu\dd{y}{x}
\end{equation}

kde $\nu$ je dynamická viskozita, můžeme určit složky tenzoru viskózního napětí
v rovnicích \ref{eq:metodika:cfd:zmena_moment}.

\smallskip

Výslednou rovnici, popisující pohyb tekutiny, poprvé navrhl v roce 1822
francouzský inženýr Claude Louise Marie Henry Navier a nezávisle v roce 1845
irský matematik George Gabriel Stokes.

\begin{equation}
 \overbrace{\rho\Big(
            \underbrace{\pdv{\mathbf{v}}{t}}_{\parbox{2cm}{\scriptsize\centering místní zrychlení}} +
            \underbrace{\mathbf{v}\cdot\nabla\mathbf{v}\vphantom{\pdv{\mathbf{v}}{t}}}_{\parbox{2cm}{\scriptsize\centering konvektivní zrychlení}}\Big)}
^{\text{Setrvačnost}} =
 \overbrace{\underbrace{-\nabla p}_{\text{Tlakový gradient}} +
            \underbrace{\mu\nabla^2\boldsymbol{\mathsf{v}}}_{\text{Viskozita}}}^{\text{Změna napětí}} + 
            \underbrace{\mathbf{f}}_{\parbox{2cm}{\scriptsize\centering Ostatní objemové síly}}
\end{equation}

Rovnice popisuje {\em rychlost} tekutiny v libovolném bodě na rozdíl od klasických
pohybových rovnic mechaniky, které obvykle popisují {\em polohu} částice.

\subsubsection{Metoda konečných objemů}

Jedna z numerických metod pro řešení parciálních diferenciálních rovnic je
metoda konečných objemů. Kontinuum je stejně jako v metodě konečných prvků
rozděleno do konečných oblastí zvané kontrolní objemy. Základní rovnice
popisující spojité prostředí (Navier-Stokesova rovnice, rovnice kontinuity a turbulence)
jsou aproximovány na těchto objemech a řešeny v jejich středech. Hodnoty
veličin na hranicích objemů se získávají interpolací. Diskretizace zaručuje
zachování proudění skrz kontrolní objem.

\smallskip

Řídící rovnice má tvar:

\begin{equation}
    \pdv{t}\iiint \vb{Q} \dd{V} + \iint \vb{F} \dd{A} = 0
\end{equation}

kde $\vb{Q}$ je vektor proměnných, $\vb{F}$ je vektor proudění,
$V$ je objem kontrolního elementu a $A$ je povrch elementu.

\smallskip

Tato metoda se využívá hlavně pro numerická řešení mechaniky tekutin a je
použita i v softwaru ANSYS-CFX.

\subsubsection{Mezní vrstva}

Experimentálně bylo zjištěno, že rychlost tekutiny na povrchu tělesa je
rovna nule. Přechod mezi laminárním a turbulentním prouděním v blízkosti
stěny je řešen v tzv.\,mezní vrstvě.

\begin{figure}[H]
    \begin{center}
        \includegraphics[width=12cm]{images/boundary_layer.png}
    \end{center}
    \caption{Proudění v mezní vrstvě}
    \label{fig:cfd:boundary_layer}
\end{figure}

Při vytváření sítě konečných objemů je nutné brát ohled na použitou funkci v mezní vrstvě (Wall function). Manuál ANSYS-CFX\cite{Ansys} doporučuje, aby vrstva obsahovala alespoň 10 uzlů. Vzdálenost prvního uzlu od stěny je dána vztahem:

\begin{equation}
    \Delta y = L \cdot \Delta y^+ \sqrt{74} \cdot Re_{L}^{-13/14}
    \label{eq:boundary:first}
\end{equation}

kde $Re_{L}$ je Reynoldsovo číslo, $L$ je délka domény (obě hodnoty lze získat z řešiče) a $\Delta y^+$ je cílová hodnota bezrozměrné veličiny $y^+$. V případě použití Automatic Wall-Function je možné použít $\Delta y^+ = 1000$.

\nomenclature{Re}{Reynoldsovo číslo}
\nomenclature{$y^+$}{Bezrozměrná vzdálenost od stěny}

Tloušťka mezní vrstvy je pak určena vztahem:

\begin{equation}
    \delta = 0.035 \cdot L \cdot Re_{L}^{-1/7}
    \label{eq:boundary:thickness}
\end{equation}

a musí platit:

\begin{equation}
    n(10) - n(1) \le \delta
    \label{eg:boundary:thickness}
\end{equation}

kde $n(x)$ je vzdálenost x-tého uzlu od stěny.



\subsubsection{Kontrola kvality sítě konečných objemů}

Numerické řešení vnáší do výsledku diskretizační chybu, kterou je možné výrazně snížit za použití sítě odpovídající kvality. Požadavky na prvky sítě jsou uvedeny v manuálu ANSYS-CFX\cite{Ansys}.

Kvalitu sítě je možné ohodnotit na základě vlastností jednotlivých prvků. Kontrolují se především tři parametry: mesh orthogonality, mesh expansion a mesh aspect ratio.

\begin{description}
    \item [Mesh orthogonality] -- určuje úhlovou odchylku dvou sousedních prvků od ideálního úhlu ($90^\circ$ pro prvky s obdélníkovou stěnou, $60^\circ$ pro prvky s trojúhelníkovou stěnou). Vysoká odchylka od ideálního úhlu, tj. příliš zkosené prvky, snižuje míru konvergence a zvyšuje diskretizační chybu.
    \begin{figure}[H]\centering
    \includegraphics[width=7cm]{images/mesh_orthogonality_measurements.jpg}
    \caption{Mesh orthogonality}
\end{figure}

    \item [Mesh expansion] -- vyjadřuje poměr mezi největším a nejmenším objevem prvku sousedících s uzlem. Prvky by se měly zvětšovat, resp. změnšovat, plynule tak, aby tento poměr byl maximálně 20.
    
    \begin{figure}[H]\centering
    \includegraphics[width=7cm]{images/mesh_expansion_measurements.jpg}
    \caption{Mesh expansion}
\end{figure}

    \item [Mesh aspect ratio] -- podobně jako mesh orthogonality popisuje tvar prvků. Je definován jako poměr největší a nejmenší plochy stran prvků sousedídích s uzlem. Nežádoucí je poměr vyšší než 100, tj. příliš ploché prvky.
    \begin{figure}[H]\centering
    \includegraphics[width=7cm]{images/mesh_aspect_ratio_measurements.jpg}
    \caption{Mesh aspect ratio}
\end{figure}

\end{description}

Preprocesory často nabízejí souhrnou analýzu sítě pod veličinou \textbf{Quality}. V~případě nevyhovující sítě může v lepším případě trvat déle, než výpočet zkonverguje. V~horším případě řešič skončí chybou a~není možné síť použít pro výpočet. Proto je nutné klást důraz na kvalitu sítě a~řídit se pokyny solveru.


\subsubsection{Rychlostní profil}

Funkce $v_m(z)$ střední rychlosti větru pro vstupní hraniční podmínku je převzata 
z~platné normy ČSN EN 1991-1-4 Eurokód 1: Zatížení konstrukcí 
-- Část 1-4: Obecná zatížení -- Zatížení větrem~\cite{EN114}.

Základní rychlost větru je určena vztahem:

\begin{equation}
v_b = c_{dir} \cdot c_{season} \cdot v_{b,0}
\label{eq:vb}
\end{equation}

\nomenclature{$v_{b,0}$}{Základní rychlost větru}

kde $c_{dir}$ součinitel směru větru a $c_{season}$ součinitel ročního období.
Hodnoty součinitelů jsou určeny v národní příloze a pro ČR jsou oba rovny 1,0.

Střední rychlost větru je pak stanovena:

\begin{equation}
v_m(z) = c_r(z) \cdot c_0(z) \cdot v_b
\label{eq:vmz}
\end{equation}

kde $c_r(z)$ je součinitel drsnosti terénu a $c_0(z)$ je součinitel orografie,
který je roven $1,0$ pokud není uvedeno jinak. Součinitel drsnosti terénu je
logaritmickou funkcí výšky nad terénem:

\begin{equation}
c(z) = \left\{ 
  \begin{array}{l l}
    k_r \cdot \ln(\displaystyle \frac{z}{z_0}) & \quad \text{pro $z_{min} \le z \le z_{max}$}\\
    c_r(z_{min}) & \quad \text{pro $z \le z_{min}$}
  \end{array} \right.
  \label{eq:cz}
\end{equation}

kde $k_r$ je součinitel drsnosti terénu:

\begin{equation}
k_r = 0,19 \left( \frac{z_0}{z_{0,II}} \right)^{0,07} 
\label{eq:kr}
\end{equation}

Parametry $z_0$ a $z_{min}$ závisí na kategorii terénu definované
v~tab.~\ref{tab:metodika:teren} a~$z_{0,II} = 0,05$ ($z_0$~pro kategorii terénu II).
Maximální výška, pro kterou platí odvozený vztah, je $z_{max} = \SI{200}{\metre}$.

\begin{table}[H]\small\centering
    \begin{tabular}{|>{\centering\arraybackslash}m{0.5cm}| m{11.5cm}| >{\centering\arraybackslash}m{1cm}| >{\centering\arraybackslash}m{1.2cm}|}
    \hline
     & Kategorie terénu & $z_0$ [m] & $z_{min}$ [m] \\ \hline
    0 & Moře nebo pobřežní oblasti vystavené otevřenému moři & 0,003 & 1 \\ \hline
    I & Jezera nebo vodorovne oblasti se zanedbatelnou vegetací a bez překážek & 0,01 & 1 \\ \hline
    II & Oblasti s nízkou vegetací jako je tráva a s izolovanými překážkami (stromy)
        budovy), jejichž vzdálenost je větší než 20násobek výšky překážek & 0,05 & 2 \\ \hline
    III & Oblasti rovnoměrně pokryté vegetací nebo budovami nebo s izolovanými překážkami,
        jejichž vzdálenost je maximálně 20násobek výšky překážek (jako jsou vesnice,
        předměstský terén, souvislý les) & 0,3 & 5 \\ \hline
    IV & Oblasti, ve kterých je nejméně 15 \% povrchu pokryto pozemními stavbami,
        jejichž průměrná výška je větší než 15 m. & 1 & 10 \\ \hline
    \end{tabular}
    \caption[Kategorie terénu a jejich parametry]{Kategorie terénu a jejich parametry\cite{EN114}}
    \label{tab:metodika:teren}
\end{table}

Kombinací vztahů \ref{eq:vb}, \ref{eq:vmz}, \ref{eq:cz} a \ref{eq:kr} získáme
výsledný vztah pro rychlostní profil:

\begin{equation}
v_m(z) = \left\{ 
    \begin{array}{l l}
        0,19 \cdot \left(\displaystyle \frac{z_0}{z_{0,II}}\right)^{0,07} \cdot \ln \left(\displaystyle \frac{z}{z_{0}}\right) \cdot v_{b,0}
            & \quad \text{pro $z_{0} < z \le z_{max}$}\\
        0 & \quad \text{pro $z \le z_{0}$}\\
    \end{array} \right.
    \label{eq:metodika:cfd:rychlost}
\end{equation}

\nomenclature{$v_m$}{Výchozí rychlost vzduchu}

Střední rychlost je použita jako logaritmická funkce na otevřeném intervalu $(z_0, z_{max})$.
Oproti vztahu v normě tedy není dole omezena hodnotou $z_{min}$. Tato úprava
je provedena z důvodu, že proudění vzduchu bude modelováno a rychlostí
profil se tedy přizpůsobí modelované tekutině.

\subsection{Přenos zatížení z ANSYS-CFX do ANSYSu}

Existují dva způsoby, jak převést tlakové zatížení z ANSYS-CFX do ANSYSu:

\begin{enumerate}
    \item Jednosměrný přenos -- V MKP modelu jsou pokryty SURF prvkem všechny plochy, které se použijí pro zatížení větrem. V CFX se provede výpočet a výsledné tlakové zatížení se importuje do ANSYSu pro konečnou statickou analýzu.

        Toto řešení má výhodu, že je lineární. Všechny zatěžovací stavy tak
        je možné libovolně kombinovat a vytvořit návrhové kombinace odpovídající
        normě.

        CFX výpočet může být v tomto případě statický (steady state) nebo
        časově závislý (transient). V obou případech je výpočet zatížení
        větrem proveden na původní, nezdeformované konstrukce a síti.
        Tímto způsobem je tedy možné modelovat širokopásmovou odezvu na
        odtrhávání vírů.

    \item Obousměrný přenos -- Cílem tohoto procesu je, aby se deformace
        konstrukce, způsobene vypočteným zatížením, projevili zpět v~deformaci
        sítě v CFX. V průběhu výpočtu tedy dochází k synchronizaci tlakovému
        zatížení v ANSYSu a deformací v~CFX.

        Nevýhoda tohoto řešení je, že není lineární. V případě, že usilujeme
        o zpřesnění výsledku pomocí této metody, měl by MKP model konstrukce
        být zatížen všemi zatěžovacími stavy.

        I v tomto případě může v CFX být proveden statický (steady state)
        výpočet, kdy se hledá rovnovážný stav deformace a tlaku, nebo
        časově závislý (transient). Tento postup je vhodný pro modelování
        úzkopásmové odezvy na odtrhávání vírů a vliv ovalování.
\end{enumerate}

Dle doporučení normy ČSN EN 1991-1-4 ~\cite{EN114} je širokopásmová odezva 
důležitější pro těžké konstrukce
(vyztužené betonové a těžké ocelové), zatímto úzkopásmová je důležitá pro 
lehké konstrukce. Z tohoto důvodu je proudění modelováno na nezdeformované
konstrukci jednosměrným přenosem zatížení.
