\chapter{معادلات حاکم و روش عددی}


\section{معادلات حاکم}

پدیده‌های کلاسیک الکترومغناطیس به وسیله معادلات ماکسول\LTRfootnote{Maxwell's Equations} قابل تشریح می‌باشند که عبارتند از:
\begin{equation}
\nabla \times \overrightarrow{E}=-\frac{\partial \overrightarrow{B}}{\partial t}\label{faraday}
\end{equation}
\begin{equation}
\nabla \times \overrightarrow{H}=\overrightarrow{j}+\frac{\partial \overrightarrow{D}}{\partial t}\label{Ampere}
\end{equation}
\begin{equation}
\nabla \cdot \overrightarrow{D}=q\label{gaussel}
\end{equation}
\begin{equation}
\nabla \cdot \overrightarrow{B}=0\label{gaussmag}
\end{equation}
معادلات \eqref{faraday} و \eqref{Ampere} به ترتیب قانون القای فارادی\LTRfootnote{Faraday Law} و قانون آمپر\LTRfootnote{Ampere Law} اصلاح شده می‌باشند و معادلات \eqref{gaussel} و \eqref{gaussmag} قانون گوس\LTRfootnote{Gauss Law} را به ترتیب در میدان‌های الکتریکی و مغناطیسی بیان می‌کنند.
جایی که $\overrightarrow{B}(Tesla)$ و $\overrightarrow{E}(V/m^3)$ به ترتیب میدان‌های مغناطیسی و الکتریکی و $\overrightarrow{H}$ و $\overrightarrow{D}$ به ترتیب میدان القایی برای میدان‌های الکتریکی و الکترومغناطیسی می‌باشند. همچنین $q (C/m^3)$ چگالی بار الکتریکی و $\overrightarrow{j}(A/m^2)$ بردار چگالی جریان الکتریکی می‌باشند.

میدان‌های القایی $\overrightarrow{H}$ و $\overrightarrow{D}$ به صورت زیر تعریف می‌شوند:
\begin{equation}
\overrightarrow{H}=\frac{1}{\mu}\overrightarrow{B}\label{H_eq}
\end{equation}
\begin{equation}
\overrightarrow{D}=\varepsilon\overrightarrow{E}\label{D_eq}
\end{equation}
که در این معادلات $\mu$ و $\varepsilon$ به ترتیب نفوذپذیری مغناطیسی و نفوذپذیری الکتریکی می‌باشند.

معادلات پیوستگی و مومنتوم برای یک جریان تراکم‌ناپذیر و رسانا به صورت زیر می‌باشد.
\begin{equation}
\nabla\cdot u=0\label{con}\\
\end{equation}

\begin{equation}
\frac{\partial u}{\partial t}+(u\cdot\nabla)u=-\frac{\nabla p}{\rho}+\nu\nabla^2u+\frac{F}{\rho}\label{mom}
\end{equation}
که در آن \lr{u} بردار سرعت، $\rho$ چگالی سیال و $\nu$ لزجت سینماتیکی جریان می‌باشد. همچنین \lr{F} نیرویی است  که در حضور میدان‌های الکترومغناطیسی به ذرات یک سیال باردار که با سرعت \lr{u} در حال حرکت می‌باشد، وارد می‌شود و نیروی لورنتز\LTRfootnote{Lorentz Force} نامیده می‌شود. 

برای مطالعه تقابل میدان جریان و میدان‌ الکترومغناطیسی لازم است که چگالی جریان را به دست آوریم. برای این منظور، به‌طور کلی دو روش وجود دارد، حل معادله القایی مغناطیسی و حل معادله پتانسیل الکتریکی.

در حالت کلی این نیرو به صورت ضرب بردار چگالی جریان الکتریکی $\overrightarrow{j}$ و بردار میدان مغناطیسی $\overrightarrow{B}$ نشان داده می‌شود:
\begin{equation}
F=\overrightarrow{j}\times\overrightarrow{B}\label{lorentz}
\end{equation}
که در آن چگالی جریان از قانون اُهم\LTRfootnote{Ohm's Law} به دست می‌آید.
\begin{equation}
J=\sigma(\overrightarrow{E}+V\times\overrightarrow{B})\label{ohm}
\end{equation}

با توجه به معادلات \ref{mom} و \ref{ohm} توزیع میدان الکترومغناطیس و توزیع کمیت‌های جریان در میدان جریان از همدیگر تاثیر می‌پذیرند. بنابراین جهت بررسی تاثیر متقابل آن‌ها، می‌بایست معادلات ماکسول و معادلات ناویر-استوکس به‌طور همزمان حل گردند. اگر سیال یک فلز مایع یا فلز مذاب با ضریب رسانایی بالا باشد، یک کنترل کامل و بدون تماس فلز مایع با استفاده از میدان‌های مغناطیسی ممکن می‌باشد. برخلاف فلزات مایع، الکترولیت‌ها مثل آب دریا، رسانایی پایینی در حدود $S/m$10 دارند. 

 در مباحث دینامیک سیالات مغناطیسی، تقریب‌های زیر برای مطالعه اثرات متقابل جریان سیال و میدان الکترومغناطیس در نظر گرفته می‌شوند\cite{crammer1973magnetofluid}.
\begin{equation}
\rho_{el}=0\label{assum1}
\end{equation}
\begin{equation}
\frac{\partial D}{\partial t}=0\label{assum2}
\end{equation}
\begin{equation}
\frac{\partial \rho_{el}}{\partial t}=0\label{assum3}
\end{equation}

از آنجا که سیال مورد نظر ما، ضریب هدایت الکتریکی پایینی دارد و با توجه به اینکه حداکثر شدت میدان مغناطیسی قابل کاربرد در عمل کم‌تر از 1 تسلا می‌باشد، از اینرو در سمت راست معادله \ref{ohm}، ترم دوم نسبت به ترم اول ناچیز بوده و قابل صرف نظر می‌باشد. با اعمال این فرض توزیع میدان الکترومغناطیسی در میدان جریان مستقل از کمیت‌های جریان می‌باشد.
\begin{equation}
J=\sigma E
\end{equation}

برای اینکه چگالی جریان برای اهداف کنترل جریانی به اندازه کافی بزرگ باشد، لازم است که میدان الکتریکی به اندازه $E_0$ را طوری اعمال کنیم که\cite{Weier2004control}:
\begin{equation}
\frac{E_0}{U_\infty B_0}\gg 1
\end{equation}
همچنین در تحقیق حاضر فرض شده است میدان مغناطیسی از به کارگیری آهنربای دائمی ایجاد شده باشد، که در این صورت تغییرات زمانی آن صفر خواهد بود و از معادله \ref{faraday} خواهیم داشت:
\begin{equation}
\nabla \times E=0 \label{af1}
\end{equation}

از معادله \ref{gaussel} و با توجه به \ref{assum1} داریم:
\begin{equation}
 \nabla \cdot E=0 \label{af2}
 \end{equation}

با توجه به معادلات \ref{af1} و \ref{af2} می‌توان نتیجه گرفت که میدان الکتریکی به صورت گرادیان یک تابع پتانسیل قابل بیان است:
\begin{equation}
E=-\nabla \phi \label{potential}
\end{equation}
همچنین با ادغام روابط \ref{af2} و \ref{potential} رابطه زیر نتیجه می‌شود:
\begin{equation}
\nabla^2 \phi=0 \label{potentia2}
\end{equation}
با اعمال اپراتور کرل در طرفین معادله \ref{Ampere} پس از بسط دادن طرفین داریم:
\begin{equation}
\nabla (\nabla \cdot H)-\nabla^2 H=-\sigma \nabla \times \nabla \phi \label{curl}
\end{equation}
با توجه به اینکه کرل گرادیان یک کمیت اسکالر صفر است و با توجه به روابط \ref{gaussmag} و \ref{H_eq} رابطه زیر به دست می‌آید:
\begin{equation}
\nabla^2 B=0 \label{divb}
\end{equation}

بنابراین معادلات حاکم بر یک جریان تراکم‌ناپذیر یک سیال رسانا با ضریب هدایت الکتریکی پایین و در حضور میدان‌های الکترومغناطیسی به معادلات (\ref{mom}،\ref{potentia2} و \ref{divb})
 به علاوه معادله پیوستگی تقلیل می‌یابد.
 
 اگر یک میدان جریان با عمق \lr{h} بر روی این عملگر در نظر گرفته شود، آن وقت شرایط مرزی الکترومغناطیسی را می‌توان به صورت زیر در نظر گرفت:
 \begin{equation}
(J_y)_{wall}=-\sigma (\frac{\partial \phi}{\partial y})=J_0 \sin{(\frac{\pi}{2a}x)}
\end{equation}
 \begin{equation}
(B_y)_{wall}=B_0 \cos{(\frac{\pi}{2a}x)}
\end{equation}
 \begin{equation}
(J_y)_h=(B_y)_h=0
\end{equation}
با اعمال شرایط مرزی فوق، توزیع میدان‌های الکترومغناطیسی از حل تحلیلی معادلات (\ref{potentia2} و \ref{divb})
به دست آمده و با اعمال فرض $h/a \gg 1$ نهایتاً نیروی لورنتز به رابطه زیر تقلیل می‌یابد\cite{berger2000turbulent}:
\begin{equation}
F=J\times B=J_0B_0exp(\frac{-\pi}{a}y)
\end{equation}
مطابق رابطه فوق، نیروی لورنتز به صورت نمایی با افزایش فاصله از دیوار کاهش می‌یابد.

یکی از مزیت‌های این روش این است که در عمل، تقریباً هر نوع توزیعی از نیروهای لورنتز را می‌توان با آرایش مناسب الکترودها و مگنت‌ها به دست آورد. علاوه‌بر این، اگر یک توزیع نیروی بهینه را برای کنترل جریان بتوان معین کرد، آن را می‌توان از حل معکوس مساله پیکربندی الکترودها و مگنت‌ها به دست آورد. این روش دارای محدودیت‌هایی، مانند انتگرال‌گیری الکترودها و مگنت‌ها برای هر شکل دلخواه می‌باشد\cite{Weier2004control}.

همان‌طور که در بخش \ref{refact} اشاره شد، برای کاهش پسا تا کنون سه نوع پیکربندی مورد تحقیق قرار گرفته است: نیروی لورنتز عمود بر دیواره، نیروی لورنتز موازی با دیواره و در جهت عرض جریان و نیروی لورنتز موازی با دیواره و در جهت جریان. در این پایان‌نامه از آخرین پیکربندی برای ایجاد نیروی لورنتز استفاده شده است. شکل \ref{thesisactuator} نحوه چیدمان این الکترودها و مگنت‌ها را برای تولید نیروی لورنتز موازی با دیوار و در جهت جریان نشان می‌دهد.
\begin{figure}
\centering
\subfigure[به کار گرفته شده بر روی صفحه تخت]{
\includegraphics[scale=0.3]{actuator}
}
\hspace*{0.5mm}
\subfigure[به کار گرفته شده بر روی ایرفویل]{
\includegraphics[scale=0.25]{airfoilactuator}
}
\caption{نحوه چیدمان الکترودها و مگنت‌ها برای ایجاد نیروی لورنتز موازی با دیواره در این پایان‌نامه}
\label{thesisactuator}
\end{figure}

از این نوع پیکربندی برای تاخیر در گذرش جریان\cite{gailitis1961on} و از بین بردن جدایش جریان و افزایش برآ\cite{weier2003control,Weier2004control} استفاده شده است. هر دو میدان الکتریکی و مغناطیسی مولفه‌هایی در جهت عمود بر دیواره (\lr{y}) و در جهت عرض جریان (\lr{z}) دارند در نتیجه نیروی لورنتز در جهت جریان (\lr{x}) اعمال می‌شود. شکل \ref{Lorentzforce} توزیع نیروهای لورنتز را نشان می‌دهد. در نزدیکی دیواره، تغییرات عرضی شدیدی در توزیع نیروی دیده می‌شود\cite{weier2003control}. 

\begin{figure}
\centering
\includegraphics[scale=0.45]{lorentzdis}
\caption{توزیع نیروی لورنتز در صفحه \lr{y-z}}
\label{Lorentzforce}
\end{figure}

در هیدرودینامیک مغناطیسی معمولاً از دو پارامتر برای مشخص کردن قدرت نیروی لورنتز استفاده می‌کنند\cite{davidson2001an}. 
\begin{enumerate}
\item
پارامتر تداخلی\LTRfootnote{Interaction Parameter}
\begin{equation}
N=\frac{J_0 B_0 c}{\rho {U_\infty}^2}
\end{equation}
که نسبت نیروهای الکترومغناطیسی به نیروهای اینرسی را نشان می‌دهد.
\item
عدد هارتمن\LTRfootnote{Hartmann Number}
\begin{equation}
Z=\frac{1}{8\pi}\frac{j_0 M_0 a^2}{\rho U_\infty \nu}
\end{equation}
\end{enumerate}

پارامتر تداخلی فقط بیشینه نیرو را بر روی دیوار محاسبه می‌کند، در حالی که نحوه کاهش آن را در جهت عمود بر دیوار در نظر نمی‌گیرد. به همین دلیل پارامتر هندسه\LTRfootnote{Geometry Parameter}\cite{Posdziech2001electromagnetic} به صورت زیر تعریف شد تا تاثیرات پیکربندی الکترودها و مگنت‌ها را در توزیع پروفیل نیروی لورنتز در نظر بگیرد.
\begin{equation}
G=\frac{a}{c}
\end{equation}

بسته به قدرت نیرو حتی می‌توان یک جت دیواره مشخص را تخمین زد. می‌توان نتیجه گرفت که نیروی لورنتز پایا همان اثراتی را بر روی کنترل جدایش جریان دارد که دمش جریان بر روی سطح مکشی ایرفویل دارد\cite{weier2003control,Weier2004control}. علاوه‌بر این در \cite{weier2003control} نشان داده شده است که این دو روش از این طریق قابل مقایسه با هم می‌شوند. به همین دلیل ضریب مومنتوم مغناطیسی طبق رابطه زیر تعریف می‌شود\cite{weier2003control}:
\begin{equation}
c_\mu=\frac{1}{4}\frac{aj_0M_0}{\rho {U_\infty}^2}\frac{x_e-x_s}{c}
\end{equation}
که در رابطه بالا $x_e-x_s$ نشان‌دهنده طول پوشیده شده توسط الکترودها و مگنت‌ها می‌باشد. $c_\mu$ مومنتوم کل تزریق شده به جریان را توسط نیروی لورنتز به فشار دینامیکی ربط می‌دهد. اگر الکترودها و مگنت‌ها کل سطح مکشی ایرفویل را بپوشانند، پارامتر تداخلی را می‌توان توسط رابطه زیر به ضریب مومنتوم مغناطیسی ربط داد:
\begin{equation}
c_\mu=\frac{a}{2c}N
\end{equation}

همچنین پارامتر تداخلی و عدد هارتمن اصلاح شده نیز مستقل از هم نیستند:
\begin{equation}
\frac{Z}{N}\propto Re
\end{equation}
 
\section{روش عددی}
برای حل عددی معادلات حاکم بر جریان سیال روش‌های مختلفی همچون اختلاف محدود، حجم محدود و المان محدود وجود دارد. هر کدام از این روش‌ مزایا و معایبی دارند. اما با توجه به توانایی کاربرد دو روش حجم محدود و به‌ویژه المان محدود در حل میدان‌های با هندسه‌های پیچیده، در دو دهه اخیر آن دو بیش از روش تفاضل محدود مورد توجه قرار گرفته‌اند.

هر چند روش المان محدود مدتی است که کارایی خود را در انطباق شبکه بر شکل‌های نامنظم نشان داده است، اما مشکلاتی مانع از پیشرفت آن شده است  \cite{patankar1387com}.

\subsection{روش عددی حجم المان محدود} 
حل معادلات جریان در بخش اول فصل نتایج که مربوط به ایرفویل شیاردار است، با برنامه توسعه داده شده به روش حجم المان محدود انجام شده است. روش حجم المان محدود اولین‌بار توسط بالیگا و پاتانکار برای برطرف کردن مشکلات فوق ارائه شد. این روش از مزیت‌های هر دو روش حجم محدود و المان محدود برخوردار است. اخیراً نیز نادری و همکاران \cite{darbandi2006multiblock,naderi2009developing,naderi2009phd} با بهره‌گیری از مزیت‌های این روش به توسعه روش حجم المان محدود برای حل میدان جریان آرام و آشفته بر روی شبکه‌ با مرز متحرک پرداخته‌اند. در ادامه به نکات مهم این روش اشاره خواهیم کرد و می‌توان برای به دست آوردن جزئیات بیش‌تر به مراجع \cite{darbandi2006multiblock,naderi2009developing,naderi2009phd} مراجعه نمود.

در این تحقیق از دیدگاه لاگرانژی-اویلری برای حل میدان جریان و با روش اجزای محدود بر مبنای حجم محدود استفاده می‌شود. از رویکرد مرتبه دوم زمانی در تقریب جملات زمانی استفاده شده است. برای تقریب جملات پخش در معادلات بخش هیدرولیک از توبع شکلی اجزای محدود استفاده شده می‌شود. جملات فشار نیز با استفاده از آن تقریب زده می‌شوند. در تمام معادلات حاکم بخشی از حملات جابجایی که به‌طور صریح متاثر از حرکت شبکه است با استفاده از قوانین بقای هندسی و رویکردهای زمانی ایجاد شده به دست می‌آید. بخش دوم جملات جابجایی هم به‌طور یکپارچه در تمام معادلات حاکم با روش بالادست و با استفاده از اثرات فیزیک جریان تقریب زده می‌شوند. با حل جریان در دامنه‌های پیچیده نشان داده می‌شود که روش ترکیبی به‌طور یکسان و دقیقی قادر به محاسبه شارها در المان‌های ثابت یا متحرک سه ضلعی و چهار ضلعی خواهیم بود. 

\begin{equation}
\frac{d}{dt}\iiint_{\vartheta (t)} \rho d\vartheta+\iint_{S(t)}\rho(V-V^{\prime}).dS=0
\end{equation}
\begin{equation}
\frac{d}{dt}\iiint_{\vartheta (t)} \rho V d\vartheta+\iint_{S(t)}\rho V(V-V^{\prime}).dS=\iint_{S(t)}\sigma \hat{n}dS+\iiint_{\vartheta (t)}F d\vartheta
\end{equation}
که در آن $V=u\hat{i}+v\hat{j}$ بردار سرعت سیال، $V^{\prime}=u^{\prime}\hat{i}+v^{\prime}\hat{j}$ بردار سرعت سطوح حجم محدود، $dS=(dS_x)\hat{i}+(dS_y)\hat{j}$ بردار عمود بر سطح سلول با بزرگی $dS=[(dS_x)^2+(dS_y)^2]$، $\hat{n}=(n_x)\hat{i}+(n_y)\hat{j}$ بردار یکه عمود بر سطح سلول و $\vartheta$ حجم مربوط به حجم محدود است. همچنین در این معادله $\sigma$ تانسور تنش کل است. تانسور تنش کل شامل فشار هیدرواستاتیک و تانسور تنش برشی است که به صورت $\sigma_{ij}=-p\delta_{ij}+\tau_{ij}$ تعریف می‌شود و \lr{F} بیانگر نیروهای حجمی می‌باشد. معادله مربوط به بقای هندسی در معادله بقای جرم نهفته است. در صورتی که سرعت ذرات سیال صفر بوده و چگالی واحد باشد از معادله بقای جرم برای هر سطح از حجم محدود خواهیم داشت:
\begin{equation}
V^{\prime}.dS\vert_m=\frac{d\vartheta}{dt}\vert_m
\label{2}
\end{equation}
که در آن \lr{m} شمارنده سطوح هر حجم محدود و $\upsilon$ حجم مربوط به هر زیر حجم محدود است. شکل \ref{CVFEM} یک حجم محدود را نشان می‌دهد که دارای 5 زیر حجم محدود \lr{SCV} بوده و محیط آن شامل 10 سطح \lr{SS} است.
\begin{figure}
\centering
\includegraphics[scale=0.6]{CVFEM}
\label{CVFEM}
\caption{یک حجم محدود شامل 5 زیر حجم محدود و 10 سطح اطراف آن}
\end{figure}

در صورتی که حرکت شبکه به نحوی انجام شود که همیشه رابطه \ref{2} برقرار باشد، بقای حجم زیرحجم محدود برقرار خواهد بود. به عبارت دیگر حرکت سطوح زیر حجم محدود باید به نحوی محاسبه شود که با قرار گرفتن در رابطه \ref{2} مقدار دقیق $d\vartheta /dt$ حاصل شود. برای تقریب مشتق زیر حجم محدود نیز می‌توان از رویکرد مرتبه دوم به صورت زیر استفاده کرد.
\begin{equation}
\frac{d\vartheta}{dt}=\frac{(3/2)\vartheta^{n+1}-2\vartheta^{n}+(1/2)\vartheta^{n-1}}{\Delta t}
\end{equation}
که در آن \lr{n} شمارنده گام‌های زمانی است. با انجام یکسری عملیات ریاضی، معادله بالا را می‌توان به صورت زیر بازنویسی کرد.
\begin{equation}
\frac{d\vartheta}{dt}=\frac{(3/2)\Delta\vartheta^{n+1}-(1/2)\Delta\vartheta^{n}}{\Delta t}
\end{equation}
که در آن $\Delta\vartheta^{n+1}=\vartheta^{n+1}-\vartheta^{n}$ و $\Delta\vartheta^{n}=\vartheta{n}-\vartheta{n-1}$ می‌باشد. در صورتی که حرکت سطوح منجر به افزایش مقدار حجم مربوط به زیر حجم محدود مورد بررسی شود مقدار $\Omega=d\vartheta /dt$ مثبت خواهد شد و به‌طور عکس منفی خواهد شد. با استفاده از این رابطه می‌توان قسمتی از معادله اندازه حرکت که ضریبی از اثرات حرکت شبکه را در خود دارد بازنویسی کرد. یعنی آنکه می‌توان نوشت:
\begin{equation}
V(V^{\prime}).dS\vert_m=V\Omega\vert_m
\end{equation}
گسسته‌سازی و تقریب‌زنی سایر جملات معادلات حاکم در فضای لاگرانژی-اویلری همانند گسسته‌سازی و تقریب‌زنی آن‌ها در فضای اویلری صورت می‌گیرد.

در این روش روی محیط ایرفویل با شبکه چهار ضلعی و کمی دورتر شبکه بی‌سازمان سه ضلعی تا محل مرزهای آزاد، ورودی و خروجی استفاده شده است. فاصله مرکز ایرفویل از مرزهای آزد و ورودی و خروجی 30 است. سرعت جریان در مرز ورودی و مرزهای آزد واحد است. مقدار فشار در خروجی نیز واحد است. عدد رینولدز جریان بر مبنای طول وتر ایرفویل و سرعت ورودی برابر 1000 است. به‌منظور بررسی حساسیت حل به شبکه محاسباتی و گام زمانی از دو شبکه ریز و درشت به ترتیب با حدود 5700 و 10600 گره استفاده شده است که، شبکه نهایی در همه حالات بیش از 9000 گره در جریان دارد و گام زمانی مورد بررسی در این تحقیق با توجه به نتایج \cite{naderi2009developing} برابر 0/05 انتخاب شده است. شکل \ref{FVEMgrid} شبکه عددی به کار رفته در روش حجم المان محدود نشان می‌دهد. همچنین برای اعتبارسنجی نتایج حاصل از جریان عبوری از روی یک سیلندر نوسانی با استفاده از روش حجم المان محدود با نتایج \cite{GUILMINEAU2002773,Yang200612} مقایسه و در شکل \ref{validcode} نشان داده شده است.
\begin{figure}
\centering
\includegraphics[scale=0.6]{validcode}
\label{validcode}
\caption{مقایسه پروفیل‌های ضریب فشار و ضریب اصطکاک پوسته سیلندر نوسانی روش عددی حاضر با نتایج جیلمن و کوتی\cite{GUILMINEAU2002773} و یانگ و بالاراس\cite{Yang200612} در فرکانس 0/8 }
\end{figure}


\begin{figure}
\centering
\subfigure{
\includegraphics[scale=0.5,angle=270]{FVEMgrid2}
}
\hspace*{0.5mm}
\subfigure{
\includegraphics[scale=0.4,angle=270]{FVEMgrid}
}
        \caption{شبکه ترکیبی تولید شده در روش حجم المان محدود}
        \label{FVEMgrid}
\end{figure}

\subsection{حل معادلات توسط نرم‌افزار فلوئنت}
برای شبیه‌سازی جریان اطراف ایرفویل نوسانی و با حضور نیروهای لورنتز، از نرم‌افزار تجاری فلوئنت استفاده شده است. برای اعمال نیروی لورنتز، یک ترم چشمه با استفاده توابع تعریف شده توسط کاربر به نرم‌افزار فلوئنت اضافه شده است. در این نرم‌افزار معادلات با استفاده از روش حجم محدود حل می‌شوند.
\subsubsection*{حل معادلات بر روی شبکه متحرک}
نرم‌افزار فلوئنت قابلیت حل جریان بر روی شبکه‌های متحرک را دارد. یکی از روش‌های متداول در این نرم‌افزار استفاده از شبکه‌های دینامیکی می‌باشد که نیاز به استفاده از تکنیک‌هایی مانند روش فنر خطی می‌باشد تا در زمان‌های بعدی، نقاط مناسب شبکه به دست آید. اما در این تحقیق با استفاده از یک استراتژی مشابه با \cite{kinsey2008parametric}، شبکه به دو ناحیه تقسیم شده است. همان‌طور که در شکل \ref{fluentgrid} دیده می‌شود، ناحیه داخلی در فاصله‌‌ای، که 5 برابر طول وتر ایرفویل می‌باشد، قرار گرفته است. قسمتی از شبکه که خارج از این ناحیه قرار داد، ثابت بوده و آن قسمتی از شبکه که درون این ناحیه قرار دارد با استفاده از یک تابع سینوسی نوسان می‌کند. همچنین فاصله ایرفویل از مرز جلویی و عقبی و بالایی (پایینی) به ترتیب 35 و 40 و 35 برابر طول وتر ایرفویل می‌باشد.

\begin{figure}[h]
\centering
\subfigure{
\includegraphics[scale=0.8,angle=270]{fluentgrid1}
}
\subfigure{
\includegraphics[scale=0.5,angle=270]{fluentgrid2}
}
\hspace*{0.5mm}
\subfigure{
\includegraphics[scale=0.5,angle=270]{fluentgrid3}
}
\caption{شبکه ترکیبی تولید شده برای نرم‌افزار فلوئنت و شرایط مرزی}
\label{fluentgrid}
\end{figure}


برای شبکه‌های متحرک فرم انتگرالی معادله بقا برای کمیت $\phi$ و یک حجم کنترل دلخواه \lr{V} که مرزهای آن متحرک می‌باشد، به صورت زیر نوشته می‌شود:
\begin{equation}
\frac{d}{dt}\int_V \rho\phi dV+\int_{\partial V}\rho\phi (\overrightarrow{u}-\overrightarrow{u_g})\cdot d\overrightarrow{A}=\int_{\partial V}\Gamma\nabla\phi\cdot\overrightarrow{A}+\int_V S_\phi\cdot dV
\label{ALE}
\end{equation}

جایی که $\rho$ چگالی سیال، $\overrightarrow{u}$ بردار سرعت سیال و $\overrightarrow{u_g}$ سرعت شبکه، شبکه متحرک می‌باشد. همچنین $\Gamma$ و $S_\phi$ به ترتیب ضریب دیفیوژن و ترم چشمه می‌باشند. این معادله در فضای اولری-لاگرانژی نوشته شده است.  اگر شبکه ثابت باشد، شکل اولرین معادله و اگر سرعت شبکه با سرعت جریان برابر باشد، شکل لاگرانژین معادله به دست می‌آید.
\begin{enumerate}
\item[1]
گسسته‌سازی ترم زمانی
\end{enumerate}
ترم زمانی در معادله \ref{ALE} با استفاده از تقریب اختلاف محدود پسرو مرتبه اول\LTRfootnote{first-order backward difference} تقریب زده می‌شود.
\begin{equation}
\frac{d}{dt}\int_V \rho\phi dV=\frac{(\rho\phi V)^{n+1}-(\rho\phi V)^n}{\Delta t}
\label{timedisc}
\end{equation}
جایی که \lr{n} و \lr{n+1} به ترتیب بیانگر کمیت‌ها در زمان جاری و زمان بعدی می‌باشند. همچنین حجم $V^{n+1}$ از رابطه زیر به دست می‌آید:

\begin{equation}
V^{n+1}=V^n+\frac{dV}{dt} \Delta t
\end{equation}
جایی که $dV/dt$ تغییرات زمانی حجم کنترل می‌باشد. برای ارضای قانون بقای هندسی، مشتق حجم نسبت به زمان یک حجم کنترل از رابطه زیر به دست می‌آید:
\begin{equation}
\frac{dV}{dt}=\int_{\partial V} \overrightarrow{u_g}\cdot d\overrightarrow{A}=\sum_{j}^{n_f}\overrightarrow{u_{g,j}}\cdot\overrightarrow{A_j}
\end{equation}
که در این رابطه، $n_f$ تعداد سطوح حجم کنترل می‌باشد. در روش شبکه‌های لغزان که در این پایان‌نامه برای حرکت شبکه در نظر گرفته شده است، حجم، حجم کنترل‌ها در زمان ثابت باقی می‌ماند،$\frac{dV}{dt}=0$، در نتیجه $V^{n+1}=V^n$ می‌شود و رابطه \ref{timedisc} به رابطه زیر تبدیل می‌شود:
\begin{equation}
\frac{d}{dt}\int_V \rho\phi dV=\frac{[(\rho\phi )^{n+1}-(\rho\phi)^n]V}{\Delta t}
\label{timedisc2}
\end{equation}
\begin{enumerate}
\item[2]
گسسته‌سازی ترم جابجایی
\end{enumerate}
نرم‌افزار فلوئنت این امکان را برای کاربر فراهم می‌سازد تا روش گسسته‌سازی را برای ترم‌های جابجایی هر یک از معادلات حاکم انتخاب نماید (دقت مرتبه دوم به طور خودکار برای ترم‌های لزج استفاده می‌شود.). وقتی که جریان با شبکه همراستا نباشد، گسسته‌سازی مرتبه اول، خطای عددی گسسته‌سازی را افزایش می‌دهد. به همین دلیل مشابه با مرجع\cite{kinsey2008parametric} از گسسته‌سازی مرتبه دوم بالادست، برای ترم‌های جابجایی استفاده شده است.

\subsubsection*{ارتباط بین میدان سرعت و فشار}
یکی از مشکلات موجود در روش‌های عددی تراکم‌ناپذیر، ارتباط بین میدان سرعت و فشار می‌باشد.  حل معادلات ناویر-استوکس، به خاطر نبودن یک معادله مستقل برای فشار (که گرادیان آن در هر سه معادله مومنتوم وجود دارد) پیچیده است. از این گذشته، معادله پیوستگی، در جریان‌های تراکم‌ناپذیر یک متغیر حاکم ندارد. پایستاری جرم، پیش از آنکه یک معادله دینامیکی باشد، یک قید سینماتیک روی میدان سرعت است. یک راه خارج شدن از این مشکل، طراحی میدان فشار با هدف تضمین ارضای معادله پیوستگی است\cite{ferziger1388computatopnal}.

در این‌ پایان‌نامه از روش تصویری (گام کسری)\LTRfootnote{Projection Method (Fractional Step)} برای رفع مشکل فوق استفاده شده است\cite{ferziger1388computatopnal,Guermond2006an}. این روش تقریباً مرسوم‌ترین روش حل‌های مستقیم (\lr{DNS}) و شبیه‌سازی گردابه‌های بزرگ می‌باشد، که در طی سال‌های گذشته به‌طور وسیعی استفاده شده است. روش به صورت رسمی در دهه 60 توسط تِمام و چورین\LTRfootnote{Temam and Chorin} معرفی شد. اساس روش در محاسبات گذرا با زمان است که این روش در نرم‌افزار فلوئنت موجود می‌باشد.

این روش که توسط کیم و معین\LTRfootnote{Kim and Moin}\cite{Kim1985application} معرفی شد، دیدگاهی است که از فشار در گام اول استفاده نمی‌کند. نقش فشار در یک جریان تراکم‌ناپذیر، اعمال پیوستگی است؛ به عبارت دیگر، فشار بیش‌تر یک متغیر ریاضی است تا یک ویژگی فیزیکی.

اساس روش بر قضیه هلمهولتز قرار دارد که بیان می‌کند هر میدان برداری مربعی انتگرال‌پذیر ($\overrightarrow{V}=L_2(\Omega)$) را می‌توان به صورت حاصل جمع دو میدان برداری نوشت که یکی بقایی و دیگری غیرچرخشی است.
\begin{equation}
\overrightarrow{V}=\overrightarrow{u}+\overrightarrow{w}
\end{equation}
بر این اساس هر میدان سرعت دلخواه را می‌توان به صورت زیر نوشت:
\begin{equation}
\overrightarrow{V}=\nabla\times\overrightarrow{q}+\nabla\times\overrightarrow{r}
\end{equation}
با مفروض کردن این قضیه، معادله مومنتوم را می‌توان به صورت زیر نوشت:
\begin{equation}
\frac{\partial u}{\partial t}+\nabla P=F(u)
\end{equation}
بر این اساس کمیت شبه فشاری به نام $\pi$ را به نحوی تعریف می‌کنیم که:
\begin{equation}
(\frac{\partial \hat{u}}{\partial t})_n=-(u^n\cdot\nabla)u^n+\nu\nabla^2u^n \label{1}
\end{equation}
\begin{equation}
(\frac{\partial u}{\partial t})_n=-\nabla \pi \label{2}
\end{equation}

ابتدا از یک میدان سرعت بقایی شروع می‌کنیم و معادله برگرز \ref{1} را در زمان انتگرال‌گیری می‌کنیم. با توجه به اینکه گرادیان فشار حذف شده بود $\hat{u}$ پیدا شده، بقایی نیست.

اگر فرض کنیم شبه فشار مناسب $\pi$ موجود است، با انتگرال‌گیری زمانی در معادله \ref{2} سرعت‌های بقایی گام \lr{n+1} پیدا خواهد شد. 

برای بقایی کردن سرعت‌ها از طرفین معادله \ref{1} دیورژانس می‌گیریم.
\begin{equation*}
\nabla \cdot{\frac{\partial {u}}{\partial t}}=-\nabla \cdot{\nabla \pi}
\end{equation*}
\begin{equation*}
\frac{\partial {\nabla \cdot u}}{\partial t}=-\nabla^2 \pi
\end{equation*}
\begin{equation*}
\frac{\nabla \cdot u^{n+1}-\nabla \cdot \hat{u}^n}{\Delta t}=-\nabla^2 \pi
\end{equation*}
\begin{equation*}
\nabla^2 \pi=\frac{\nabla \cdot \hat{u}^n}{\Delta t}
\end{equation*}
بر این اساس یک گام روش تصویری شامل مراحل زیر است:
\begin{enumerate}
\item
معادله برگرز در زمان انتگرال‌گیری می‌شود تا $u^{n+1}$ به دست آید.
\item
$\hat{u}$ به دست آمده بقایی نیست و لذا با دیورژانس گرفتن از آن معادله پوآسون حل می‌شود تا $\pi^{n+1}$ به دست آید.
\item
سرعت‌های فیزیکی (بقایی) گام \lr{n+1} از رابطه اصلاحی به دست می‌آید.
\item
در صورت نیاز معادله پوآسون فشار حل می‌شود تا مقادیر فشار نیز محاسبه شود در غیر این صورت به گام قبل بر می‌گردیم.
\end{enumerate}

تفاوت اصلی بین روش تصویری و روش تصحیح فشار از نوع \lr{SIMPLE} این است که در اولی، معادله فشار (یا تصحیح فشار) در هر گام زمانی یک بار حل می‌شود، در حالی که در دومی، هر دو دسته معادلات مومنتوم و تصحیح فشار در هر گام زمانی بارها حل می‌شوند (تکرارهای خارجی). این اختلاف بیش‌تر به این دلیل است که روش‌های تصویری اصولاً برای شبیه‌سازی جریان‌های ناپایدار استفاده می‌شوند، در حالی که روش‌های تصحیح فشار برای محاسبه جریان‌های پایدار به کار می‌روند. به دلیل اینکه در روش‌هایی از نوع \lr{SIMPLE}، پایستاری جرم فقط در انتهای گام زمانی اعمال می‌شود، لازم نیست که معادله تصحیح فشار در هر مرحله از تکرارهای خارجی به‌طور دقیق حل شود (در بیش‌تر موارد کاهش باقیمانده تا مرتبه اول کفایت می‌کند). در حقیقت در جریان‌های پایدار، ارضای کامل شرط پیوستگی فقط در مرحله همگرایی لازم است. در شبیه‌سازی جریان‌های ناپایدار به منظور ارضای بقای جرم در هر گام زمانی، معادله فشار (یا تصحیح فشار) باید با دامنه خطایی محدودتر حل شود. معمولاً برای حل معادله پوآسون فشار در شبیه‌سازی جریان‌های ناپایدار یا هندسه‌های پیچیده، معادلات خطی با استفاده از روش‌های گرادیان مزدوج حل می‌شوند\cite{ferziger1388computatopnal}.

در بررسی حاضر همچون \cite{kinsey2008parametric}، گام زمانی از رابطه زیر به دست آمده است که بسیار کوچک است. این در حالی است که در روش حجم المان محدود مقدار گام زمانی را 0/05 در نظر گرفته‌ایم. گفتنی است که در اینجا هم سرعت ورودی جریان و طول وتر ایرفویل واحد است و متغیرهای زمانی با این دو بی‌بعد می‌شوند.
\begin{equation}
\Delta t=\frac{T}{2000}
\end{equation}

همچنین به منظور اعتبارسنجی روش عددی مورد نظر، نتایج حاصله از این روش با نتایج عددی به دست آمده از روش حجم المان محدود که اعتبار آن‌ قبلاً ثابت شده است، در شکل \ref{validatefluent} مقایسه شده است.
\begin{figure}
\centering
\includegraphics[scale=0.5]{validate}
\caption{اعتبارسنجی نتایج به دست آمده از نرم‌افزار فلوئنت}
\label{validatefluent}
\end{figure}