Commit fca8368d authored by Jonathan Lambrechts's avatar Jonathan Lambrechts

imex doc

parent 1370259c
File added
%% Creator: Inkscape inkscape 0.91, www.inkscape.org
%% PDF/EPS/PS + LaTeX output extension by Johan Engelen, 2010
%% Accompanies image file 'dessin.pdf' (pdf, eps, ps)
%%
%% To include the image in your LaTeX document, write
%% \input{<filename>.pdf_tex}
%% instead of
%% \includegraphics{<filename>.pdf}
%% To scale the image, write
%% \def\svgwidth{<desired width>}
%% \input{<filename>.pdf_tex}
%% instead of
%% \includegraphics[width=<desired width>]{<filename>.pdf}
%%
%% Images with a different path to the parent latex file can
%% be accessed with the `import' package (which may need to be
%% installed) using
%% \usepackage{import}
%% in the preamble, and then including the image with
%% \import{<path to file>}{<filename>.pdf_tex}
%% Alternatively, one can specify
%% \graphicspath{{<path to file>/}}
%%
%% For more information, please see info/svg-inkscape on CTAN:
%% http://tug.ctan.org/tex-archive/info/svg-inkscape
%%
\begingroup%
\makeatletter%
\providecommand\color[2][]{%
\errmessage{(Inkscape) Color is used for the text in Inkscape, but the package 'color.sty' is not loaded}%
\renewcommand\color[2][]{}%
}%
\providecommand\transparent[1]{%
\errmessage{(Inkscape) Transparency is used (non-zero) for the text in Inkscape, but the package 'transparent.sty' is not loaded}%
\renewcommand\transparent[1]{}%
}%
\providecommand\rotatebox[2]{#2}%
\ifx\svgwidth\undefined%
\setlength{\unitlength}{480.39390384bp}%
\ifx\svgscale\undefined%
\relax%
\else%
\setlength{\unitlength}{\unitlength * \real{\svgscale}}%
\fi%
\else%
\setlength{\unitlength}{\svgwidth}%
\fi%
\global\let\svgwidth\undefined%
\global\let\svgscale\undefined%
\makeatother%
\begin{picture}(1,0.29480557)%
\put(0,0){\includegraphics[width=\unitlength,page=1]{dessin.pdf}}%
\put(0.50000611,0.26107532){\color[rgb]{0,0,0}\makebox(0,0)[lb]{\smash{$u_s$}}}%
\put(0.50000611,0.02793328){\color[rgb]{0,0,0}\makebox(0,0)[lb]{\smash{$u$}}}%
\end{picture}%
\endgroup%
This diff is collapsed.
\documentclass{paper}
\usepackage[mathletters]{ucs}
\usepackage[utf8x]{inputenc}
\usepackage{amsmath}
\DeclareUnicodeCharacter{183}{{\,\cdot\,}}
\newcommand{\s}{_{\text{s}}}
\renewcommand{\ss}[1]{_{\text{s}#1}}
\newcommand{\re}{\mathit{R\kern-.15em e}}
\newcommand{\acc}{a}
\newcommand{\tia}{\tilde {\acc}}
\newcommand{\tiv}{\tilde v}
\newcommand{\aim}{\tia_{\text{im}}}
\newcommand{\aex}{\tia_{\text{ex}}}
\newcommand{\aext}{\acc_{\mathrm{ext}}}
\newcommand{\af}{\acc_{\text{f}}}
\newcommand{\haf}[1]{\hat\acc_{\text{f#1}}}
\newcommand{\vf}{v_{\text{f}}}
\newcommand{\vfree}{v_{\text{free}}}
\usepackage{xcolor}
\newcommand{\nb}[1]{{\color{green!30!black}\begin{itemize}#1\end{itemize}}}
\usepackage{a4}
\usepackage{graphicx}
\usepackage[round]{natbib}
\begin{document}
\bibliographystyle{plainnat}
\section{Avancement des particules}
Pour le moment, on résout en déplacement (si j'ai bien compris ce qu'on entendait par là). C'est à dire que:
\begin{align*}
a^i &= a(x^i, v^i)\\
\vfree^{i+1} &= v^i + a^i\,Δt\\
v^{i+1} &= \mathrm{Contacts}\left(x^i, v^i, \vfree^{i+1}, θ\right)\\
x^{i+1} &= x^i + ((1-θ)\,v^i + θ\,v^{i+1})\,Δt
\end{align*}
où les contacts sont résolus de manière à ce qu'il n'y ait pas d'interpénétration dans $x^{i+1}$.
Jusqu'ici, j'ai pris $θ = 0$ mais passer à $θ = 0.5$ ne devrait pas poser de problème.
Ce schéma fonctionne mais on est limité à des petits pas de temps (autant pour la précision que pour la stabilité semble-t-il). Afin de s'affranchir de cette limite, un schéma multi-pas partiellement implicite (IMEX RK) sera utilisé pour calculer $\vfree$.
\nb{
\item À la place des deux dernières lignes, est-ce ceci ?
\begin{align*}
v^{i+1} &= \mathrm{Contacts}\left(x^i, (1-θ)\,v^i + θ\, \vfree^{i+1}\right)\\
x^{i+1} &= x^i + v^{i+1}\, Δt
\end{align*}
ça ne change rien sur un pas de temps mais le $v^{i+1}$ final est différent ce qui aura une influence sur le pas de temps suivant.
\item Je m'interroge sur l'opportunité d'intégrer $x$ dans le schéma multi pas mais la résolution des contacts pose problème, à moins de l'appliquer à tous les pas explicites, puis de nouveau sur la combinaison linéaire finale.
}
\section{Accélération partiellement implicite}
\begin{align*}
\acc(x, v) &= \af(x, v) + \aext\\
\af(x, v) &= -β∇p -α(v-\vf)\,(v-\vf)
\end{align*}
Il est nécessaire de résoudre le problème de (Navier-)Stokes pour trouver $∇p$ et $\vf$.
\begin{align*}
0 &= ∇·(v_f + P(v))\\
0 &= D(v_f) -∇p + P(ρ_s\,\af)
\end{align*}
$P$ est l'opérateur qui projette sur le maillage une variable définie sur les particules.
Comme ce problème fait intervenir $\af$, on a un système à résoudre. Ce système est non linéaire en $\vf$, même si on traite $v$ de manière explicite.
Lorsque $α$ est grand, le pas de temps stable explicite devient très petit, on va donc chercher à traiter le terme $α(v-\vf)\,(v-\vf)$ de manière partiellement implicite (en $v$). Pour ce faire, on choisit une approximation:
\[
\tia(x, v, \tiv) = -β∇p -α(\tiv-\vf)\,(v-\vf) + \aext
\]
telle que $\tia(x, v, v) = a(x, v)$ et
\[
\acc(x,v) = \underbrace{\tia(x, v, \tiv)}_{\displaystyle \aim(x, v, \tiv)} + \underbrace{\tia(x, v, v) - \tia(x, v, \tiv)}_{\displaystyle \aex(x, v, \tiv)}
\]
$v$ est traité implicitement dans l'approximation $\aim$ et explicitement dans la correction $\aex$. Comme $\tiv$, on prend la vitesse au début du pas de temps.
\nb{\item si on applique RK à $x$, il faut aussi introduire un $\tilde x$ pour la partie implicite}
\section{Formalisme IMEX}
Je suis ici les notations d'\cite{ascher_implicit-explicit_1995}: $u$ est le vecteur inconnu, $g$ la partie implicite et $f$ la partie explicite.
Pour un pas de temps, on a:
\begin{description}
\item évaluer $\hat K_1 = f(u^n)$
\item pour $i=1,\dots,s $ :
\begin{description}
\item résoudre
$ K_i = g(u^n_i) $
et
$ u^n_i = u^n + Δt \sum_{j=1}^i\left(a_{i,j}\,K_j+\hat a_{i+1,j}\,\hat K_j\right)$
\item évaluer $\hat K_{i+1} = f(u_i)$
\end{description}
\item $u^{n+1} = u^n + Δt\left(\sum_{j=1}^sb_j\,K_j + \sum_{j=1}^{s+1}\hat b_j\,\hat K_j\right)$
\end{description}
NB: Lorsque $\hat b_{s+1} = 0$, l'évaluation de $\hat K_{s+1}$ n'est pas nécessaire.
\section{Application}
Dans notre cas, vu que la linéarisation se fait autour de $v^n$ et que $x$ ne rentre pas dans le schéma RK, on a $\hat K_1 = 0$, ce qui réduit de 1 le nombre d'évaluation de la partie explicite.
Le $u^{n+1}$ désigne en fait $\vfree^{n+1}$ et la seule résolution des contacts a lieu, à la fin de l'itération.
\subsection{Forward-backward Euler (1,1,1)}
\[
\text{implicite: }\begin{array}{c|cc}
0 & 0 & 0 \\
1 & 0 & 1 \\
\hline
& 0 & 1
\end{array},
\hspace{1cm} \text{explicite: } \begin{array}{c|cc}
0 & 0 & 0 \\
1 & 1 & 0 \\
\hline
& 1 & 0
\end{array}
\]
Parmi les $\hat K_i, K_i$, le seul non nul utile est $K_2 = g(u^n_2)$ et $u^{n+1} = u^n_2$. Pour finir, il reste simplement:
\[
\vfree^{n+1} = v^n + Δt\,g(\vfree^{n+1})
\]
\[g(\vfree^{n+1}) = \aim(x^n,\vfree^{n+1},v^n)= -β∇p -α(v^n-\vf)\,(\vfree^{n+1}-\vf) + \aext\]
On peut isoler $\vfree^{n+1}$
\[
\vfree^{n+1} = \frac{v^n + Δt(-β∇p +α(v^n-\vf)\,\vf + \aext)}{1 + Δt\,α(v^n-\vf)}
\]
On résout le problème de Stokes en utilisant:
\[\af = -β∇p - α(v^n-\vf)\,(\vfree^{n+1}-\vf)
= \frac{-β∇p-α(v^n-\vf)(v^n - \vf + Δt\,\aext)}{1 + Δt\,α(v^n-\vf)}\]
Et on obtient la vitesse libre:
\[\vfree^{n+1} = v^{n} + Δt\, (\af + \aext)\]
\subsection{Forward-backward Euler (1,2,1)}
\[
\text{implicite: }\begin{array}{c|cc}
0 & 0 & 0 \\
1 & 0 & 1 \\
\hline
& 0 & 1
\end{array},
\hspace{1cm} \text{explicite: } \begin{array}{c|cc}
0 & 0 & 0 \\
1 & 1 & 0 \\
\hline
& 0 & 1
\end{array}
\]
On ajoute $\hat K_2$.
\begin{description}
\item $u^n_2 = u^n + Δt\,g(u^n_2)$
\item $u^{n+1} = u^n + Δt\left[g(u^n_2) + f(u^n_2)\right]$
\end{description}
Comme dans le cas (1,1,1), on résout le problème de Stokes en utilisant:
\[{\af}_2 = \frac{-β∇p-α(v^n-\vf)(v^n - \vf + Δt\,\aext)}{1 + Δt\,α(v^n-\vf)}\]
On obtient la vitesse intermédiaire $v_2$:
\[v_2 = v^{n} + Δt\, ({\af}_2 + \aext)\]
À partir de laquelle on peut calculer la correction explicite:
\begin{align*}
\aex(x, v_2, v^n)
&= \tia(x, v_2,v_2) - \tia(x, v_2,v^n)\\
&= \haf{2} - {\af}_2
\end{align*}
$\haf{2}$ est obtenu en résolvant le problème de Stokes en utilisant $v = v_2$.
\[\haf{2} = -β∇p-α(v_2-\vf)(v_2 - \vf)\]
On trouve alors $\vfree^{n+1}$ à partir du {$\haf{2}$} explicite :
\[\vfree^{n+1} = v^n + Δt\,({\af}_2 + \aext + \haf{2} - {\af}_2) = v^n + Δt\,(\haf{2} + \aex)\]
\subsection{Implicit-Explicit midpoint (1,2,2)}
\[
\text{implicite: }\begin{array}{c|cc}
0 & 0 & 0 \\
1/2 & 0 & 1/2 \\
\hline
& 0 & 1
\end{array},
\hspace{1cm} \text{explicite: } \begin{array}{c|cc}
0 & 0 & 0 \\
1/2 & 1/2 & 0 \\
\hline
& 0 & 1
\end{array}
\]
Même coût que le précédent mais second ordre de précision.
\begin{description}
\item $u^n_2 = u^n + \frac{Δt}{2}\,g(u^n_2)$
\item $u^{n+1} = u^n + Δt\left[g(u^n_2) + f(u^n_2)\right]$
\end{description}
Les calculs sont identiques, il faut juste utiliser $\frac{Δt}{2}$ au lieu de $Δt$ pour le calcul de ${\af}_2$ et $v_2$.
\subsection{plus haut ordre de précision}
On va déjà essayer les précédents ...
\bibliography{zotero}
\end{document}
\documentclass{paper}
\usepackage[mathletters]{ucs}
\usepackage[utf8x]{inputenc}
\usepackage{amsmath}
\DeclareUnicodeCharacter{183}{{\,\cdot\,}}
\newcommand{\s}{_{\text{s}}}
\renewcommand{\ss}[1]{_{\text{s}#1}}
\newcommand{\re}{\mathit{R\kern-.15em e}}
\usepackage{a4}
\usepackage{color}
\usepackage{graphicx}
\begin{document}
\paragraph{Instable}
\begin{align*}
&\left\{
\begin{aligned}
0 &= ∇·(u + u^n\s)\\
0 &= ∇p + D(u, c\s) + F(u, u^n\s, c^n\s, p)
\end{aligned} \right.\\
&u\s^{n+1}, c\s^{n+1} = C(u\s^{n}, c\s^{n}, V F(u, u^n\s, c^n\s, p))
\end{align*}
$F$ sont les forces d'interaction fluides-particules par unité de volume, $V$ le volume des particules.
\paragraph{Stable mais inconsistant}(Code actuel)
\begin{align*}
&\left\{
\begin{aligned}
0 &= ∇·(u + u^n\s + \frac{Δt}{ρ\s}c\s F(u, u^n\s, c^n\s, p))\\
0 &= ∇p + D(u, c\s) + F(u, u^n\s, c^n\s, p)
\end{aligned} \right.\\
&u\s^{n+1}, c\s^{n+1} = C(u\s^{n}, c\s^{n}, V F(u, u^n\s, c^n\s, p))
\end{align*}
$u^n\s + \frac{Δt}{ρ\s}c\s F$ est une approximation de $u_{\text{s free}}^{n+1}$. Cette approximation ignore les projections. En outre, elle ne tient compte que des forces du fluides, en particulier les contacts ainsi que les forces extérieures (la gravité) sont ignorées. C'est pourquoi, cette formulation est stable mais inconsistante. (On pourrait sans doute la rendre consistante en incluant les forces extérieurs et les réactions de contacts au temps $n$.)
\paragraph{Stable et consistant}(mais Newton inexacte)
\[ \left\{
\begin{aligned}
u\s^{n+1}, c\s^{n+1} &= C(u\s^{n}, c\s^{n}, F(u, u\s, c\s, p))\\
&= ∇·(u + u^{n+1}\s)\\
0 &= ∇p + D(u, c\s) + F(u, u^n\s, c^n\s, p)
\end{aligned} \right.\]
Mais $C$ doit être résolu à chaque itération du Newton (ce qui peut être cher avec LMGC) et sa linéarisation est problématique.
Cependant, cette linéarisation peut être inexacte, au prix d'une moins bonne convergence du Newton.
L'idée est donc d'utiliser $u\s^{n+1}$ exact mais une linéarisation approchée (qui ignore en tout cas les contacts et sans doute les projections).
\[
\frac{}{∂U_i}\big<u\s^{n+1}∇ψ_j\big>
= \sum_k\frac{∂u^{n+1}\s}{∂u}\bigg|_k ψ_{ik}∇ψ_{jk}w_k
\]
En négligeant les contacts et les projections (et le déplacement des particules $c\s$ en $u$):
\[
\frac{∂u\s^{n+1}}{∂u} \approx \frac{Δt}{ρ\s} c\s\frac{∂F}{∂u}
\]
TODO : il faudrait garder les projections de manière à réutiliser $\frac{∂F}{∂u}\big|_p$ et à obtenir une meilleure approximation à moindre coût.
\paragraph{Approche Directe}
Je pense que ça ne peut pas marcher sans tenir compte des contacts dans la linéarisation.
\newpage
\paragraph{condition de stabilité dans le cas 1d homogène}\ \\
Dans un cas où la vitesse et la compacité seraient homogènes, on aurait :
{\centering \def\svgwidth{\columnwidth} \input{dessin.pdf_tex}}
par conservation, on a :
\[u = - u_s\]
si on ne considère que le terme de trainée :
\[G_p = C_d|u_p|u_p\frac{πd²__f}{2ε^{1.8}}\]
$C_d$ est une fonction du nombre de Reynolds de la particule
\[\re_p =d__f(u_p-u_f)εμ_f^{-1}\]
par exemple:
\begin{align*}
C_d &= \bigg(0.63 + \frac{4.8}{\re_p^{0.5}}\bigg
\end{align*}
si $\re_p < 25$, le second terme domine et on a :
\[G_p ≃ \frac{5μ}{d__f2|u_p|ε}|u_p|u_p\frac{πd²__f}{2ε^{1.8}}4 μ ε^{-2.8}d_p u_p \]
prenons $ε=0.3$
\[G_p ≃ 150 μd_pu_p\]
avec un schéma d'Euler explicite, ça donne :
\[ u_p^{n+1} = u_p^n - Δt\frac{G_p^n}{m_p}\left(1 - \frac{150μd_p}{m_p}Δt\right)u^n_p\]
on trouve donc comme condition de stabilité
\[Δt < \frac{m_p}{150μd_p}\]
Ce qui se vérifie expérimentalement (eg pour depot 2d, malgré les nombreuses approximations, on trouve 0.023 pour dt max et en pratique, ça explose à partir de 0.016, ce qui est proche).
\paragraph{Solution : utiliser un schéma de Patankar (orth ?)} \ \\
sans perte de généralité, on peut écrire :
\[F_p = -V∇P - α_p|u_p|u_p \]
de sorte qu'on a :
\[\hat{u}_{p}^{n+1} = u_p^n -\frac{Δt}{m_p}\left(V∇P^n + α^n_p|u_p^n|\hat{u}_{p}^{n+1}\right)\]
$\hat{u}$ désigne la vitesse libre
on peut facilement isoler $\hat u^{n+1}$
\[\left(1 +\frac{Δt}{m_p}α^n_p|u_p^n|\right) \hat{u}_{p}^{n+1} = u_p^n -\frac{Δt}{m_p}V∇P^n\]
et itérer de cette manière les particules (on pourrait aussi envisager un schéma $Θ$).
Pour ne pas changer l'implémentation, on peut écrire
\begin{align*}
\hat{u}_{p}^{n+1}
&= \frac{u_p^n -\frac{Δt}{m_p}V∇P^n}{1 +\frac{Δt}{m_p}α^n_p|u_p^n|}\\
&= u_p^n + \frac{Δt}{m_p} \left(-u_p^n\frac{m_p}{Δt} + \frac{m_p}{Δt} \frac{u_p^n -\frac{Δt}{m_p}V∇P^n}{1 +\frac{Δt}{m_p}α^n_p|u_p^n|}\right)\\
&= u_p^n + \frac{Δt}{m_p} \left(-u_p^n\frac{m_p}{Δt} + \frac{\frac{m_p}{Δt} u_p^n -V∇P^n}{1 +\frac{Δt}{m_p}α^n_p|u_p^n|}\right)\\
&= u_p^n - \frac{Δt}{m_p} \left(\frac{α^n_p|u_p^n|u_p^n + V∇P^n}{1 +\frac{Δt}{m_p}α^n_p|u_p^n|}\right)\\
\end{align*}
Ce qui revient à réduire le pas de temps d'un facteur ridicule, il y a un problème qqpart. En plus, il faut garder $u_p - u_f$ et pas jute $u_p$.
Une autre solution (plus précise ?) serait de sous-itérer le schéma de particule en updatant la force.
\end{document}
This diff is collapsed.
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