Commit a120ddb2 authored by Matthieu Constant's avatar Matthieu Constant
Browse files

add missing term + doc equations

parent be58d5fd
Pipeline #3712 failed with stage
in 14 seconds
\documentclass[12pt]{paper}
\usepackage[utf8]{inputenc}
\usepackage{amsmath,amssymb}
\usepackage{bm}
\usepackage[a4paper,left=2cm,right=2cm,top=2cm,bottom=2cm]{geometry}
\usepackage[dvipsnames]{xcolor}
%\makeatletter
%\newcount\bracketnum
%\newcommand\makecolorlist[1]{%
% \bracketnum0\relax
% \makecolorlist@#1,.%
% \bracketnum0\relax
%}
%\def\makecolorlist@#1,{%
% \advance\bracketnum1\relax
% \expandafter\def\csname bracketcolor\the\bracketnum\endcsname{\color{#1}}%
% \@ifnextchar.{\@gobble}{\makecolorlist@}%
%}
%\let\oldleft\left
%\let\oldright\right
%\def\left#1{%
% \global\advance\bracketnum1\relax
% \colorlet{temp}{.}%
% \csname bracketcolor\the\bracketnum\endcsname
% \oldleft#1%
% \color{temp}%
%}
%\def\right#1{%
% \colorlet{temp}{.}%
% \csname bracketcolor\the\bracketnum\endcsname
% \oldright#1%
% \global\advance\bracketnum-1\relax
% \color{temp}%
%}
%\makeatother
\usepackage[english]{babel}
\setlength{\parindent}{0pt}
%\makecolorlist{black,red,blue,ForestGreen,orange,Magenta}
\usepackage{graphics}
\newcommand{\langlC}{\thicklines\begin{picture}(7,7)
\put(1.1,2.5){\rotatebox{60}{\line(1,0){15}}}
\put(1.1,2.5){\rotatebox{300}{\line(1,0){15}}}
\end{picture}}
\newcommand{\ranglC}{\thicklines\begin{picture}(7,7)
\put(-2,2.5){\rotatebox{120}{\line(1,0){15}}}
\put(-2,2.5){\rotatebox{240}{\line(1,0){15}}}
\end{picture}}
\newcommand{\langlCp}{\thicklines\begin{picture}(7,7)
\put(1.1,2.5){\rotatebox{60}{\line(1,0){7}}}
\put(1.1,2.5){\rotatebox{300}{\line(1,0){7}}}
\end{picture}}
\newcommand{\ranglCp}{\thicklines\begin{picture}(7,7)
\put(.2,2.5){\rotatebox{120}{\line(1,0){7}}}
\put(.2,2.5){\rotatebox{240}{\line(1,0){7}}}
\end{picture}}
\newcommand{\langlen}{\Huge\langle}
\newcommand {\dphi}[2]{\dfrac{\partial\phi_{#1}}{\partial #2}}
\newcommand{\dcdt}{ \langlC\dfrac{C_k^{old}-C_k}{\Delta t}\phi_k\ranglC}
\newcommand{\dudt}{ \langlC\dfrac{U_k^{old}-U_k}{\Delta t}\phi_k\ranglC}
\newcommand{\dvdt}{ \langlC\dfrac{V_k^{old}-V_k}{\Delta t}\phi_k\ranglC}
\newcommand{\cn}{\langlCp C_k\phi_k\ranglCp}
\newcommand{\dcdx}{\langlC\dphi{k}{x}C_k\ranglC}
\newcommand{\dcdy}{\langlC\dphi{k}{y}C_k\ranglC}
\newcommand{\qn}{\langlCp Q_k\phi_k\ranglCp}
\newcommand{\pn}{\langlCp P_k\phi_k\ranglCp}
\newcommand{\dpn}[1]{\langlC\dphi{k}{#1}P_k\ranglC}
\newcommand{\un}{\langlCp U_k\phi_k\ranglCp}
\newcommand{\vn}{\langlCp V_k\phi_k\ranglCp}
\newcommand{\dudx}{\langlC \dphi{k}{x}U_k\ranglC}
\newcommand{\dudy}{\langlC \dphi{k}{y}U_k\ranglC}
\newcommand{\dvdx}{\langlC \dphi{k}{x}V_k\ranglC}
\newcommand{\dvdy}{\langlC\dphi{k}{y}V_k\ranglC}
\newcommand{\divu}{\langlC \dphi{k}{x}U_k+\dphi{k}{y}V_k\ranglC}
\newcommand{\utaux}{\un \left(\dudx+ \dcdx \dfrac{\un}{\cn}\right)+\vn \left(\dudy+ \dcdy \dfrac{\un}{\cn}\right)}
\newcommand{\utauy}{\un \left(\dvdx+ \dcdx \dfrac{\vn}{\cn}\right)+\vn \left(\dvdy+ \dcdy \dfrac{\vn}{\cn}\right)}
\newcommand{\nb}{\bm{\nabla}}
\newcommand{\ub}{\bm{u}_k}
\newcommand{\ib}{\bm{I}}
\newcommand{\sigmab}{\bm{\sigma}_k}
\newcommand{\gb}{\bm{g}}
\newcommand{\dpartial}[2]{\dfrac{\partial#1}{\partial#2}}
\begin{document}
Balance equations resulting from applying Navier-Stokes equations to a control volume containing $N_f$ fluid phases are the following:
\begin{align*}
N_f\times\left\lbrace\begin{aligned}
\rho_k\dpartial{\ub}{t}+\rho_k\nb\cdot\left(\dfrac{\ub\ub}{q_k}\right)\\
\dpartial{q}{t}+\nb\cdot\ub+\varepsilon\nb\cdot \bm{r}_k\end{aligned}\right.&\begin{aligned}
&=\nb\cdot\sigmab q_k+q_k\rho_k \gb\phantom{\dpartial{\ub}{t}}\\
&=0\end{aligned}&\begin{aligned}\phantom{\dpartial{\ub}{t}}\text{phase momentum equations}\\
\phantom{\dpartial{q}{t}}\text{phase mass equations}\end{aligned}\\
1-c-\sum_{k=1}^{N_f}q_k&\begin{aligned}&=0\phantom{\dpartial{q}{t}}\end{aligned}&\begin{aligned}\text{global volume conservation equation}\end{aligned}
\end{align*}
where $\bm{r}_k$ is the residual of the momentum equation used to stabilised the $\mathbb{P}_1-\mathbb{P}_1$ finite element formulation and the Cauchy stress tensor of the fluid phase $k$ can be written:
\[\sigmab=-p\ib+\mu_k\left(\nb\dfrac{\ub}{q}+\left(\nb\dfrac{\ub}{q}\right)^T\right)\]
It is important to note that the velocity appearing in the above equations is the mean velocity of the fluid phase under consideration and can be written as the product of the point phase velocity and the volume fraction of the fluid phase. The resiudal stabilisation introduces an unknown parameter $\varepsilon$ that could change the solution. The Pressure Stabilisation Petrov-Galerkin (PSPG) is the simplest way to stabilise the finite element formulation by linking all the unknown together to prevent the development of high frequency pressure modes. PSPG method of order one consists in only keeping the normal stress part of the Cauchy stress tensor:
\[\bm{r}_k=-\nb\cdot q_k p\ib \]
The weak form of these equations is obtained by integrating it multiplied by some test function over the whole domain:
\begin{align*}
N_f\times\left\lbrace\begin{aligned}
\langlC\rho_k\dpartial{\ub}{t}\psi_i\ranglC+\langlC\rho_k\nb\cdot\left(\dfrac{\ub\ub}{q_k}\right)\psi_i\ranglC\\
\langlC\dpartial{q}{t}\psi_i\ranglC+\langlC\nb\cdot\ub\psi_i\ranglC+\langlC\varepsilon\nb\cdot\bm{r}_k\ranglC\end{aligned}\right.&\begin{aligned}
&=\langlC\nb\cdot\sigmab q_k\psi_i\ranglC+\langlC q_k\rho_k \gb\psi_i\ranglC\phantom{\dpartial{\ub}{t}}\\
&=0\phantom{\dpartial{q}{t}}\end{aligned}\\
\langlC \left(1-c-\sum_{k=1}^{N_f}q_k\right)\psi_i\ranglC&\begin{aligned}&=0\end{aligned}
\end{align*}
These equations can be simplified by integrating by part and applying the divergence theorem to the flux terms:
\begin{align*}
N_f\times\left\lbrace\begin{aligned}
\langlC\rho_k\dpartial{\ub}{t}\psi_i\ranglC+\langlC\rho_k\nb\cdot\left(\dfrac{\ub\ub}{q_k}\right)\psi_i\ranglC\\
\langlC\dpartial{q}{t}\psi_i\ranglC+\langlC\nb\cdot\ub\psi_i\ranglC\end{aligned}\right.&\begin{aligned}
&=-\langlC\sigmab\cdot\nb q_k\psi_i\ranglC+\langlC\langlC\sigmab q_k\psi_i\cdot\bm{n}\ranglC\ranglC+\langlC q_k\rho_k \gb\psi_i\ranglC\phantom{\dpartial{\ub}{t}}\\
&=\langlC\varepsilon\bm{r}_k\cdot\nb\psi_i\ranglC+\langlC\langlC\varepsilon\bm{r}_k\psi_i\cdot\bm{n}\ranglC\ranglC\phantom{\dpartial{q}{t}}\end{aligned}\\
\langlC \left(1-c-\sum_{k=1}^{N_f}q_k\right)\psi_i\ranglC&\begin{aligned}&=0\end{aligned}
\end{align*}
The fully developed 2D equations considering that the vertical velocity is $v$ while the horizontal velocity is $u$ are given by:
\renewcommand{\arraystretch}{2}
\begin{center}
\begin{tabular}{c|c|c|c}
Source terms&Flux terms&Layer terms&\\
\hline
$\langlC\bm{f}_{\bm{u}_k}^0,\psi_i\ranglC$ &$+\langlC\bm{f}_{\bm{u}_k}^1,\nb\psi_i\ranglC$ &$+\langlC\langlC\bm{f}_{\bm{u}_kb},\psi_i\ranglC\ranglC $&$=0$\\[1em]
$\langlC\bm{f}_{q_k}^0,\psi_i\ranglC$&$+\langlC\bm{f}_{q_k}^1,\nb\psi_i\ranglC$&$+\langlC\langlC\bm{f}_{q_k b},\psi_i\ranglC\ranglC$&$=0$\\[1em]
$\langlC\bm{f}_{q_p}^0,\psi_i\ranglC$&&&$=0$\\[1em]
\end{tabular}
\end{center}
\renewcommand{\arraystretch}{1}
\begin{align*}
\bm{f}_{\bm{u}}^0&=\rho_k\left( \ub+ \ub\dfrac{\nb\cdot\ub}{q}+\dfrac{\bm{A}_k\cdot\ub}{q}\right)-p\ib\cdot\nb q_k+\textcolor{red}{\mu\dfrac{\bm{B}_k\cdot\nb q_k}{q_k}}-q_k\rho_k\bm{g}\\[1em]
\bm{f}_{\bm{u}}^1&= -q_kp\ib+\mu \bm{B}_k\\[1em]
\bm{f}_{\bm{ub}}&=\bm{B}_k\cdot\bm{n}_k\\[1em]
\bm{f}_{q_k}^0&=\dpartial{q}{t}+\nb\cdot\ub\\[1em]
\bm{f}_{q_k}^1&=\varepsilon\left(\textcolor{green!70!black}{p\ib\cdot\nb q_k}+q_k\nb p\cdot\ib\right)\\[1em]
\bm{f}_{q_k b}&=\varepsilon\left(p\ib\cdot\nb q_k+q_k\nb p\cdot\ib\right)\cdot\bm{n}_k\\[1em]
\bm{f}_{q_p}^0&=1-c-\sum_{k=1}^{N_f}q_k
\end{align*}
where the matrices are defined as follow:
\begin{equation*}
\ub=\begin{pmatrix}
u_k\\
v_k
\end{pmatrix},\qquad\bm{A}_k=\begin{pmatrix}
\dpartial{u_k}{x}-\dfrac{u_k}{q}\dpartial{q}{x}&\dpartial{u_k}{y}-\dfrac{u_k}{q}\dpartial{q}{y}\\
\dpartial{v_k}{x}-\dfrac{v_k}{q}\dpartial{q}{x}&\dpartial{v_k}{y}-\dfrac{v_k}{q}\dpartial{q}{y}\\
\end{pmatrix},\qquad \bm{B}_k=\bm{A}_k+\bm{A}_k^T
\end{equation*}
%\begin{align*}
%\rho_k\dpartial{u_k}{t}\psi_i&+\dfrac{\rho_k}{q_k}\left(\dfrac{u_k}{q_k}\dpartial{u}{x}+\dfrac{u_k}{q_k}\dpartial{v}{y}+u_k\dpartial{u_k}{x}+v_k\dpartial{u_k}{y}-\dfrac{u_ku_k}{q_k}\dpartial{q_k}{x}-\dfrac{v_ku_k}{q_k}\dpartial{q_k}{y}\right)\psi_i\\
%&=q_kp\dpartial{\psi_i}{x}+p\dpartial{q}{x}\psi_i-2\mu \left(\dpartial{u}{x}-\dfrac{u}{q}\dpartial{q}{x}\right)\dpartial{\psi_i}{x}-\mu\left(\dpartial{u}{y}+\dpartial{v}{x}-\dfrac{u}{q}\dpartial{q}{y}-\dfrac{v}{q}\dpartial{q}{x}\right)\dpartial{\psi_i}{y}\\
%&-2\mu \dfrac{\psi_i}{q}\left(\dpartial{u}{x}-\dfrac{u}{q}\dpartial{q}{x}\right)\dpartial{q}{x}-\mu\dfrac{\psi_i}{q}\left(\dpartial{u}{y}+\dpartial{v}{x}-\dfrac{u}{q}\dpartial{q}{y}-\dfrac{v}{q}\dpartial{q}{x}\right)\dpartial{q}{y}\\
%%&-2\mu\dfrac{\psi_i}{q}\left(\dpartial{u}{x}-\dfrac{u}{q}\dpartial{q}{x}\right)\dpartial{q}{x}\\
%%&-\mu\dfrac{\psi_i}{q}\left(\dpartial{u}{y}+\dpartial{v}{x}-\dfrac{u}{q}\dpartial{q}{y}-\dfrac{v}{q}\dpartial{q}{x}\right)\dpartial{q}{y}
%&+2\mu \left(\dpartial{u}{x}-\dfrac{u}{q}\dpartial{q}{x}\right)\psi_in_x+\mu\left(\dpartial{u}{y}+\dpartial{v}{x}-\dfrac{u}{q}\dpartial{q}{y}-\dfrac{v}{q}\dpartial{q}{x}\right)\psi_i n_y
%\end{align*}
\end{document}
\ No newline at end of file
......@@ -357,6 +357,7 @@ static void fluid_problem_f(FluidProblem *problem, double *f0, double *f1, doubl
//f0[Q] += epsilonp*dq[i]*dp[i];
f0[U+i] = rho*((u[i]-uold[i])/dt + u[i]*oq2*divu + utau[i]*oq2) -p*dq[i]+ (i==(D-1) ? -q*rho*problem->g : 0);
for (int j = 0; j < D; ++j) {
f0[U+i] += mu*(tau[i][j]+tau[j][i])*dq[j];
f1[(U+j)*D+i] = mu*(tau[i][j]+tau[j][i]) + (i==j ? -q*p : 0);
}
f1[Q*D+i] = epsilonq*dq[i]-u[i]+q*epsilonp*dp[i];
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment