P301191
/
iNSE-ALE-Article
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity
iNSE-ALE-Article/tex/1_submission/paper.tex

508 lines
37 KiB
TeX
Raw Permalink Blame History

%%
%% Copyright 2019-2020 Elsevier Ltd
%%
%% This file is part of the 'CAS Bundle'.
%% --------------------------------------
%%
%% It may be distributed under the conditions of the LaTeX Project Public
%% License, either version 1.2 of this license or (at your option) any
%% later version. The latest version of this license is in
%% http://www.latex-project.org/lppl.txt
%% and version 1.2 or later is part of all distributions of LaTeX
%% version 1999/12/01 or later.
%%
%% The list of all files belonging to the 'CAS Bundle' is
%% given in the file `manifest.txt'.
%%
%% Template article for cas-sc documentclass for
%% single column output.
%\documentclass[a4paper,fleqn,longmktitle]{cas-sc}
\documentclass[a4paper,fleqn]{cas-sc}
\usepackage[numbers]{natbib}
%\usepackage[authoryear]{natbib}
%\usepackage[authoryear,longnamesfirst]{natbib}
\usepackage[normalem]{ulem}
%%%Author macros
\def\tsc#1{\csdef{#1}{\textsc{\lowercase{#1}}\xspace}}
\tsc{WGM}
\tsc{QE}
\tsc{EP}
\tsc{PMS}
\tsc{BEC}
\tsc{DE}
%%%
\newtheorem{proposition}{Proposition}
\newtheorem{corollary}[proposition]{Corollary}
\newtheorem{theorem}{Theorem}
\newtheorem{lemma}[theorem]{Lemma}
\newdefinition{remark}{Remark}
\newproof{proof}{Proof}
\begin{document}
\let\WriteBookmarks\relax
\def\floatpagepagefraction{1}
\def\textpagefraction{.001}
\shorttitle{Incompressible flows in moving domains}
\shortauthors{Ar\'ostica and Bertoglio}
%\begin{frontmatter}
\title[mode = title]{On monolithic and Chorin-Temam schemes for incompressible flows in moving domains}
\tnotetext[1]{This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No 852544 - CardioZoom).}
\author[1]{Reidmen Ar\'{o}stica}
\ead{r.a.arostica.barrera@rug.nl}
%\ead[url]{https://sites.google.com/view/cvmath/people/ar\'{o}stica?authuser=0}
%\credit{Conceptualization, Analysis, Implementation}
\author[1]{Crist\'{o}bal Bertoglio}
\ead{c.a.bertoglio@rug.nl}
%\ead[url]{https://sites.google.com/view/cvmath/people/bertoglio?authuser=0}
%\credit{Conceptualization, Methodology, Guidance}
\address[1]{Bernoulli Institute, University of Groningen, Groningen, 9747AG, The Netherlands}
%\cortext[cor1]{Corresponding author}
%\cortext[cor2]{Principal corresponding author}
\begin{abstract}
%When working with numerical schemes for the incompressible Navier-Stokes equations (iNSE), several schemes can be proposed and diverse techniques be used, from monolithic approaches to Chorin-Temam methods. We study a standard but general monolithic scheme that characterizes several approaches used in literature, where unconditional energy stability will holds under constrains with the addition of consistent terms. With such finding, we derive a Chorin-Temam scheme and show the unconditional stability of it.
Several time discretization schemes for the incompressible Navier-Stokes equations (iNSE) in moving domains have been proposed. Here we introduce them in a unified fashion, allowing a common well possedness and time stability analysis. It can be therefore shown that only a particular choice of the numerical scheme ensures such properties. The analysis is performed for monolithic and Chorin-Temam schemes. Results are supported by numerical experiments.
\end{abstract}
%\begin{graphicalabstract}
%%\includegraphics{figs/grabs.pdf}
%\end{graphicalabstract}
%\begin{highlights}
%\item Research highlights item 1
%\item Research highlights item 2
%\item Research highlights item 3
%\end{highlights}
\begin{keywords}
incompressible flows \sep moving domains \sep monolithic coupling \sep Chorin-Temam \sep stability analysis %\\
\end{keywords}
% EDITING COMMANDS
\newcommand{\del}[1]{\sout{#1}}
\renewcommand{\mod}[2]{\sout{#1}$\mapsto$#2}
\newcommand{\RA}[1]{{\color{magenta}RA: #1}}
\newcommand{\CB}[1]{{\color{blue}CB: #1}}
\maketitle
\section{Introduction}
%When working with flows from a simulation point of view, several schemes are proposed depending on its applications, suitable for specific requirement e.g. high spatio-temporal resolution, fidelity, stability or fast implementation/simulation times.
%The literature is vast, but to our knowledge there is not a overview trying to summarize their approaches in a single scheme.
Several works have been reported dealing with the numerical solution of the iNSE in moving domains within an Arbitrary Lagrangian Eulerian formulation (ALE), primarily in the context of fluid-solid coupling.
In particular different choices of time discretization have been reported on \cite{Basting2017, Murea2016, Landajuela2016,Lozovskiy2018,smaldone2014,langer2014numerical,LeTallec2001,Liu2018,failer2020impact,Hessenthaler2017}.
% In \cite{Basting2017, Murea2016, Landajuela2016} the FSI within Arbitrary Lagrangian-Eulerian (ALE) formalism in a linearized approach is studied, and \cite{Lozovskiy2018,smaldone2014} also include mass conservation terms. \cite{langer2014numerical,LeTallec2001,Liu2018} propose fully non-linear FSI schemes, and extensions using Crank-Nicholson are used by \cite{failer2020impact,Hessenthaler2017}. \CB{No entiendo la clasificacion de los papers, me parece que le falta un poco de coherencia. No deberian tambien haber mas papers, o estos son todos? }
%Thougthful
To the best of the authors knowledge, only a few monolithic schemes have been thoroughly analyzed, e.g. in \cite{Lozovskiy2018, Burtschell2017, LeTallec2001, smaldone2014}, while no analysis has been reported for Chorin-Temam (CT) methods.
The goal of this work is therefore to assess well-posedness and unconditional energy balance of the iNSE-ALE for all reported monolithic and CT discretization schemes within a single formulation.
% Maybe we need to add the works of \cite{Boffi2004} for general surveys of ALE schemes.
%State-of-the art: monolithic (what is proven), CT (not proven, just "used") \\
%, which as it will be seen later on, holds under some restrictions with help of consistent stabilization terms.
%The findings for the monolithic case are then used to introduce a Chorin-Temam scheme, for which unconditional energy stability holds.
The reminder of this paper is structured as follows:
Section \ref{sec:continuous_problem} provides the continuous problem that will be studied. Section \ref{sec:monolithic_schemes} introduces a general monolithic scheme that characterizes several approaches used in literature, well-posedness and energy stability are studied and discussed. Section \ref{sec:chorin_temam_schemes} introduces the Chorin-Temam schemes where time stability is analyzed. Finally, Section \ref{sec:numerical_examples} provides numerical examples testing our results.
\section{The continuous problem}
\label{sec:continuous_problem}
%\CB{propongo sacar el $$, en el paper no tiene mucho sentido, solo tiene sentido en FSI. Ademas asi se va a achicar el texto}
In the following, let us consider a domain $\Omega_0 \subset \mathbb{R}^d$ with $d = 2,3$ and a deformation mapping $\mathcal{X}: \mathbb{R}^d\times \mathbb{R}_{+}\mapsto\mathbb{R}^d$ that defines the time evolving domain $\Omega_t := \mathcal{X}(\Omega_0, t)$. We assume $\mathcal{X}$ a continuous mapping, 1-to-1 with continuous inverse. We denote $\mathbf{X} \in \mathbb{R}^d$ the cartesian coordinate system in $\Omega_0$ and $\mathbf{x}_t := \mathcal{X}(\mathbf{X}, t)$ the one in $\Omega_t$, by $F_t := \frac{\partial \mathbf{x}_t}{\partial \mathbf{X}}$ the deformation gradient, $F^{-1}_t$ its inverse and $J_t := det(F_t)$ its jacobian.
Similarly, $Grad(\mathbf{\mathbf{g}}) := \frac{\partial \mathbf{g}}{\partial \mathbf{X}}$, $Div(\mathbf{g}) := \frac{\partial}{\partial \mathbf{X}} \cdot \mathbf{g}$ denote the gradient and divergence operators, respectively, and $\epsilon^t(\mathbf{g}) := \frac{1}{2}( Grad(\mathbf{g})F_t^{-1} + F_t^{-T} Grad(\mathbf{g})^{T})$ the symmetric gradient, for $\mathbf{g}$ a well-defined vector function. We consider the weak form of the iNSE in ALE form \cite[Ch. 5]{Richter2017}: Find $(\mathbf{u}(t), p(t)) \in \mathbf{H}^{1}_0(\Omega_0) \times L^2_0(\Omega_0)$ with $\mathbf{u}(0) = \mathbf{u}_{init}$ s.t.
\begin{equation}
\label{eq:continuous_formulation}
\begin{aligned}
% \int_{\Omega_0} \rho J^t \partial_{t} \mathbf{u} \cdot \mathbf{v} + \int_{\Omega_0} \rho J^t Grad(\mathbf{u}) F_t^{-1} (\mathbf{u} - \mathbf{w}) \cdot \mathbf{v} + \int_{\Omega_0} J^t 2\mu \, \epsilon^t(\mathbf{u}):\epsilon^t(\mathbf{v}) - \int_{\Omega_0} Div(J^t F_t^{-1} \mathbf{v}) p & \\
% \textcolor{red}{ + \alpha \int_{\Omega_0} \frac{\rho}{2} \left(\partial_{t} J^t - Div\left(J^t F_t^{-1} \mathbf{w} \right) \right) \mathbf{u} \cdot \mathbf{v} + \beta \int_{\Omega_0} \frac{\rho}{2} Div\left(J^t F_t^{-1} \mathbf{u}\right) \mathbf{u} \cdot \mathbf{v} }
% + \int_{\Omega_0} Div(J^t F_t^{-1} \mathbf{u})q
% & = 0 \\
\int_{\Omega_0} \rho J^t \partial_{t} \mathbf{u} \cdot \mathbf{v} + \rho J^t Grad(\mathbf{u}) F_t^{-1} (\mathbf{u} - \mathbf{w}) \cdot \mathbf{v} + J^t 2\mu \, \epsilon^t(\mathbf{u}):\epsilon^t(\mathbf{v}) - Div(J^t F_t^{-1} \mathbf{v}) p + Div(J^t F_t^{-1} \mathbf{u})q = 0%& \\
\end{aligned}
\end{equation}
for all $(\mathbf{v}, q) \in \mathbf{H}^1_0(\Omega_0) \times L^2_0(\Omega_0),\, t > 0$, $\mathbf{u}_{init}, \mathbf{w}(t) \in \mathbf{H}^1_0(\Omega_0)$ given initial and domain velocities, where $\mathbf{H}^1_0(\Omega_0)$ is the standard Sobolev space of vector fields $\mathbf{u}$ defined on $\Omega_0$ with values in $\mathbb{R}^d$ s.t. $\mathbf{u} = \mathbf{0}$ on $\partial \Omega_0$, $L^2_0(\Omega_0)$ the standard square integrable space of functions $p$ s.t. $p(\mathbf{0}) = 0$.
Notice that the flow at time $t$ is given by $\mathbf{u} \circ \mathcal{X}^{-1}(\cdot, t)$.
%
%\begin{remark}
%Formulation \eqref{eq:continuous_formulation} includes strongly consistent \textit{Geometric Conservation Law} (GCL) and \textit{Temam} stabilization terms multiplied by scalars $\alpha$ and $\beta$ respectively. This allows for a unified analysis for all reported methods.
%\end{remark}
\begin{remark}
Although Dirichlet boundary conditions are used throughout this work, it can be extended straightforwardly to the Neumann case by including the so called \textit{backflow stabilizations}, see e.g. \cite{bertoglio2018benchmark}. Moreover, in the discrete case the extension of well-posedness results to the case with non-zero boundary conditions follow from trace theorem.
\end{remark}
%An energy estimate can be deduced from \eqref{eq:continuous_formulation} as done in \cite[Chap. 9]{Quarteroni2009Cardiovascular}.
\begin{proposition}
\label{prop:energy_continuous} \cite[Chap. 9]{Quarteroni2009Cardiovascular}
Provided $(\mathbf{u}(t), p(t))$ a solution of the formulation \eqref{eq:continuous_formulation}, the following energy balance holds:
\begin{equation}
\label{eq:continuous_energy_estimate}
%\begin{aligned}
\partial_{t} \int_{\Omega_0} \frac{\rho}{2} J^t \vert \mathbf{u} \vert^2 = - \int_{\Omega_0} J^t 2\mu \vert \epsilon^t (\mathbf{u}) \vert^2 %- (\alpha - 1) \int_{\Omega_0} \frac{\rho}{2} \partial_{t} J^t - \int_{\Omega_0} \frac{\rho}{2} Div \left( J^t F^{-1}_t \left( (\beta - 1) \mathbf{u} - (\alpha -1) \mathbf{w} \right)\right)
%\end{aligned}
\end{equation}
\end{proposition}
\begin{remark}
Proposition \ref{prop:energy_continuous} makes use of the \textit{Geometric Conservation Law} (GCL) $\partial_t J^t = Div\left( J^t F_{t}^{-1} \mathbf{w}\right)$.
\end{remark}
\section{Monolithic schemes (first order in time)}
\label{sec:monolithic_schemes}
%Different \CB{time} discretization approaches for Problem \eqref{eq:continuous_formulation} can be taken into account, nevertheless if we fix a first order discretization for the time-derivative term, the schemes used in the literature can be characterized within a single formulation.
Most of the numerical schemes for Problem \eqref{eq:continuous_formulation} reported in the literature are first order and can be written as follows.
Given a conforming finite element space $\mathbf{V} \times Q$ of $\mathbf{H}^1_0(\Omega_0) \times L^2_0(\Omega_0)$ for velocity and pressure fields, the discrete problem of interest reads:
Find $(\mathbf{u}^{n+1}, p^{n+1}) \in \mathbf{V} \times Q$ s.t.
\begin{equation}
\label{eq:discretized_monolithic_formulation}
\mathcal{A}(\mathbf{u}^{n+1},\mathbf{v}) - \mathcal{B}(\mathbf{v},p^{n+1}) + \mathcal{B}(\mathbf{u}^{n+1},q) = F(\mathbf{v}) \quad \forall (\mathbf{v}, q) \in \mathbf{V} \times Q
\end{equation}
being
\begin{equation}
\label{eq:lhs_bilinear_form_A}
\begin{aligned}
\mathcal{A} (\mathbf{u}, \mathbf{v}) := & \int_{\Omega_0} \rho \frac{J^{\star\star}}{\tau} \mathbf{u} \cdot \mathbf{v} + \int_{\Omega_0} \rho J^{\star} Grad(\mathbf{u}) F^{-1}_{\star} (\mathbf{u}^{\ast} - \mathbf{w}^{\ast\ast}) \cdot \mathbf{v} + \int_{\Omega_0} J^{\star} 2\mu \epsilon^{\star}(\mathbf{u}):\epsilon^{\star}(\mathbf{v}) \\
& + \alpha \int_{\Omega_0} \frac{\rho}{2} \left( \frac{J^{n+1} - J^{n}}{\tau} - Div\left( J^{\star} F^{-1}_{\star} \mathbf{w}^{\ast\ast}_{h} \right) \right) \mathbf{u} \cdot \mathbf{v} + \beta \int_{\Omega_0} \frac{\rho}{2} Div\left( J^{\star} F^{-1}_{\star} \mathbf{u}^{\ast} \right) \mathbf{u} \cdot \mathbf{v}
\end{aligned}
\end{equation}
with $\alpha, \beta \in \{0, 1\}$ given parameters, and
%The matrix $B \in \mathbb{R}^{m^{Q}\times m^{\mathbf{V}}}$ and vector $F_n \in \mathbb{R}^n$ correspond also to the discretization of the bilinear and linear forms:
\begin{equation}
\label{eq:remaining_forms}
\begin{aligned}
\mathcal{B}(\mathbf{u}, q) & := \int_{\Omega_0} Div\left( J^{\star} F^{-1}_{\star} \mathbf{u} \right) q \quad \forall q \in Q, \quad \mathcal{F}(\mathbf{v}) := \int_{\Omega_0} \rho \frac{J^n}{\tau} \mathbf{u}^n \cdot \mathbf{v} \quad \forall \mathbf{v} \in \mathbf{V}
\end{aligned}
\end{equation}
\begin{remark}
The term multiplying $\alpha$ is the discrete residual of GCL, while the one multiplying $\beta$ is a strongly consistent term vanishing for incompressible velocity fields.
\end{remark}
%\begin{equation}
%\label{eq:discretized_monolithic_formulation}
%\begin{aligned}
% \int_{\Omega_0} \rho J^{\star \star} \frac{\mathbf{u}^{n+1} - \mathbf{u}^n}{\tau} \cdot \mathbf{v} + \int_{\Omega_0} \rho J^{\star} Grad(\mathbf{u}^{n+1}) F_{\star}^{-1} (\mathbf{u}^{\ast} - \mathbf{w}^{\ast \ast}) \cdot \mathbf{v} +\int_{\Omega_0} J^{\star} 2\mu \epsilon^{\star} (\mathbf{u}^{n+1}):\epsilon^{\star} (\mathbf{v}) & \\
% - \int_{\Omega_0} Div(J^{\star} F_{\star}^{-1} \mathbf{v}) p^{n+1} + \int_{\Omega_0} Div(J^{\star} F_{\star}^{-1} \mathbf{u}^{n+1}) q & \\
% \color{red}{ + \alpha \int_{\Omega_0} \frac{\rho}{2} \left( \frac{J^{n+1} - J^n}{\tau} - Div\left( J^{\star} F_{\star}^{-1}\mathbf{w}^{\ast\ast} \right) \right) \mathbf{u}^{n+1} \cdot \mathbf{v} + \beta \int_{\Omega_0} \frac{\rho}{2} Div\left( J^{\star} F_{\star}^{-1} \mathbf{u}^{\ast} \right) \mathbf{u}^{n+1} \cdot \mathbf{v} } & = 0
%\end{aligned}
%\end{equation}
%Within \CB{Formulation \eqref{eq:discretized_monolithic_formulation}}, an explicit domain treatment is done by \cite{Basting2017} whom linearize the convective term with $(\star, \star\star, \ast, \ast\ast) = (n, n, k, n)$, where $k$ the Newton's iteration index while \cite{Murea2016} takes $(\star, \star\star, \ast, \ast \ast) = (n, n, n, n)$, both using $\alpha = \beta = 0$; similarly \cite{Lozovskiy2018} use an implicit coupling with $(\star, \star\star, \ast, \ast \ast) = (n, n+1, n, n+1)$ whereas \cite{smaldone2014} uses an explicit one with $(\star, \star\star,\ast, \ast\ast) = (n, n+1, n, n)$ where $\alpha = \beta = 1$.
%If the domain treatment is implicit \cite{Langer2016} takes the approach with $(\star, \star\star, \ast, \ast\ast) = (n+1, n+1, n+1, n+1), \, \alpha = \beta = 0$ by applying a Newton solver for the fully nonlinear problem, same formulation is also taken by \cite{LeTallec2001}, \cite{Wall2009} and \cite{Liu2018} in the context of space-dependent densities with different choices of $\alpha, \beta$; if the coupling is explicit \cite{Landajuela2016} takes $(\star, \star\star, \ast, \ast\ast) = (n+1, n+1, n, n+1)$ with $\alpha = \beta = 0$ providing a linearization of the convective term, while \cite{Wang2020} utilizes an alternative weak formulation where the domain velocity $\mathbf{w}$ is obtain as solution to a specific weak form with $(\star, \star\star, \ast) = (n+1, n, n+1)$, $\alpha = \beta = 1$.
Formulation \eqref{eq:discretized_monolithic_formulation} contains a wide family of reported methods:
\begin{itemize}
\item Using $\alpha = \beta = 0$: $(\star, \star\star, \ast, \ast\ast) = (n, n, n+1, n)$ is used in \cite{Basting2017}, $(\star, \star\star, \ast, \ast \ast) = (n, n, n, n)$ in \cite{Murea2016} and $(\star, \star\star, \ast, \ast\ast) = (n+1, n+1, n+1, n+1)$ in \cite{Langer2016}, and $(\star, \star\star, \ast, \ast\ast) = (n, n+1, n, n+1)$ in \cite{Landajuela2016}.
\item Using $\alpha = \beta = 1$: $(\star, \star\star, \ast, \ast \ast) = (n, n+1, n, n+1)$ in \cite{Lozovskiy2018}, $(\star, \star\star,\ast, \ast\ast) = (n, n+1, n, n)$ in \cite{smaldone2014} and $(\star, \star\star, \ast, \ast\ast) = (n+1, n, n+1, n+1)$ in \cite{LeTallec2001, Wang2020}.
\end{itemize}
\begin{remark}
Other methods such as second order approximations can be found in \cite{Tallec2003, Liu2018} and Crank-Nicolson approaches in \cite{failer2020impact, Hessenthaler2017}.
% and more general surveys in \cite{Wall2009}.
% Notice that, whenever we take $\star \star = n+1$ the energy equality associated to the problem cannot be easily stabilized with consistent terms.
\end{remark}
%\subsection{Well posedness at each time step}
%In the following, well-posedness of Problem \eqref{eq:discretized_monolithic_formulation} is assessed using the standard techniques for saddle-point problems. As standard in literature, the problem is studied based on its matrix formulation that reads for each time step:
%\begin{equation}
%\label{eq:matrix_formulation}
%\begin{bmatrix}
%A & -B^T \\
%B & 0
%\end{bmatrix}
%\begin{bmatrix}
%U^{n+1} \\
%P^{n+1}
%\end{bmatrix}
%=\begin{bmatrix}
%F_n \\
%0
%\end{bmatrix},\quad U^{0} = U_{init}
%\end{equation}
%
%where $(U^{n+1}, P^{n+1}) \in \mathbb{R}^{m^{\mathbf{V}}} \times \mathbb{R}^{m^{Q}}$ denote the coefficients associated to the finite element basis with $(m^{\mathbf{V}}, m^{Q})$ the degrees of freedom pair for velocity/pressure elements and $U_{init}$ the coefficients associated the discretization of $\mathbf{u}_{init}$.
%\CB{The matrix $A_n \in \mathbb{R}^{m^{\mathbf{V}}\times m^{\mathbf{V}}}$ corrresponds to the discretization of the bilinear form:} % \rightarrow \mathbb{R}$ is given by
%\begin{equation}
%\label{eq:lhs_bilinear_form_A}
%\begin{aligned}
%\mathcal{A} (\mathbf{u}, \mathbf{v}) = & \int_{\Omega_0} \rho \frac{J^{\star\star}}{\tau} \mathbf{u}^{n+1} \cdot \mathbf{v} + \int_{\Omega_0} \rho J^{\star} Grad(\mathbf{u}^{n+1}) F^{-1}_{\star} (\mathbf{u}^{n} - \mathbf{w}^{\ast\ast}) \cdot \mathbf{v} + \int_{\Omega_0} J^{\star} 2\mu \epsilon^{\star}(\mathbf{u}^{n+1}):\epsilon^{\star}(\mathbf{v}) \\
%& + \alpha \int_{\Omega_0} \frac{\rho}{2} \left( \frac{J^{n+1} - J^{n}}{\tau} - Div\left( J^{\star} F^{-1}_{\star} \mathbf{w}^{\ast\ast}_{h} \right) \right) \mathbf{u}^{n+1} \cdot \mathbf{v} + \beta \int_{\Omega_0} \frac{\rho}{2} Div\left( J^{\star} F^{-1}_{\star} \mathbf{u}^{n} \right) \mathbf{u}^{n+1} \cdot \mathbf{v}
%\end{aligned}
%\end{equation}
%The matrix $B \in \mathbb{R}^{m^{Q}\times m^{\mathbf{V}}}$ and vector $F_n \in \mathbb{R}^n$ correspond also to the discretization of the bilinear and linear forms:
%\begin{equation}
%\label{eq:remaining_forms}
%\begin{aligned}
%B(\mathbf{u}, q) & := \int_{\Omega_0} Div\left( J^{\star} F^{-1}_{\star} \mathbf{u} \right) q \quad \forall q \in Q, \quad F_n(\mathbf{v}) := \int_{\Omega_0} \rho \frac{J^n}{\tau} \mathbf{u}^n \cdot \mathbf{v} \quad \forall \mathbf{v} \in \mathbf{V}
%\end{aligned}
%\end{equation}
\begin{proposition}
\label{prop:monolithic_schemes}
By assuming well-posed, orientation-preserving deformation mappings, i.e. $(J^n)_{n \in \mathbb{N}}$ bounded in $L^{\infty}(\Omega_0)$, $J^n > 0$ for each $n \geq 0$, Problem \eqref{eq:discretized_monolithic_formulation} has unique solution for inf-sup stable finite element spaces if $\left( 2J^{\star\star} + J^{n+1} - J^{n} \right) > 0$ and $\alpha = \beta = 1$.
\end{proposition}
\begin{proof}
%The matrix system \eqref{eq:matrix_formulation} is solvable whenever $A_n^{-1}$ exists and $BA^{-1}_nB^{T}$ is invertible.
Since all operators are bounded, and inf-sup stable elements are used for velocity and pressure, it is enough to ensure that the bilinear form $\mathcal{A}$ is coercive.
%Thus its enough to ensure $\text{rank}(B) = m$ and $A_n$ positive definite.
%As standard in literature, let us evaluate \eqref{eq:lhs_bilinear_form_A} in $\mathbf{v} = \mathbf{u}$ and $q = p$. Integrating by parts the convective term and joining expressions, the following equality holds:
Indeed:
\begin{equation}
\label{eq:proof_monolithic_case}
\begin{aligned}
\mathcal{A}(\mathbf{u}, \mathbf{u}) = & \int_{\Omega_0} \frac{J^{\star}}{2\tau} \left( \frac{2 J^{\star\star}}{J^{\star}} + \alpha \frac{J^{n+1} - J^{n}}{J^{\star}} \right) \vert \mathbf{u} \vert^2 + % \int_{\Omega_0}
J^{\star} 2\mu \vert \epsilon^{\star}(\mathbf{u}) \vert^2 %\\
%&
+ %\int_{\Omega_0}
\frac{\rho}{2} Div\left( J^{\star} F^{-1}_{\star} \left( (\beta -1) \mathbf{u}^{\ast} - (\alpha-1) \mathbf{w}^{\ast\ast} \right) \right) \vert \mathbf{u} \vert^2
\end{aligned}
\end{equation}
being the last quantity strictly positive under the stated assumptions.
\qed
\end{proof}
\begin{corollary}
Assuming $\alpha = \beta = 1$, Problem \eqref{eq:discretized_monolithic_formulation} is well posed when:
\begin{itemize}
\item $3J^{n+1} - J^{n} > 0$ if $\star \star = n+1$, i.e. a restriction on the time step size.
\item $J^{n+1} + J^{n} > 0$ if $\star \star = n$, i.e. no restriction on the time step size.
\end{itemize}
No restrictions apply to $\star,\ast, \ast\ast$.
\end{corollary}
%\subsection{Energy balance}
\begin{proposition}
\label{prop:energy_estimate_monolithic}
Under assumptions of Proposition \ref{prop:monolithic_schemes} and $\alpha = \beta = 1, \star\star = n$, the scheme \eqref{eq:discretized_monolithic_formulation} is unconditionally stable.
\end{proposition}
\begin{proof}
By setting $\mathbf{v} = \mathbf{u}^{n+1}$ in the bi-linear form \eqref{eq:lhs_bilinear_form_A}, $q = p^{n+1}$ in forms \eqref{eq:remaining_forms} and manipulating terms as standard in literature, the energy equality follows:
\begin{equation}
\begin{aligned}
\int_{\Omega_0} \rho \frac{J^{n+1}}{2\tau} \vert \mathbf{u}^{n+1} \vert^2 - \int_{\Omega_0} \rho \frac{J^{n}}{2\tau} \vert \mathbf{u}^n \vert^2 = & \int_{\Omega_0} \frac{\rho}{2\tau} (J^{n+1} - J^{\star\star}) \vert \mathbf{u}^{n+1} \vert^2 + \int_{\Omega_0} \frac{\rho}{2\tau} (J^{\star\star} - J^n) \vert \mathbf{u}^n \vert^2 - \int_{\Omega_0} 2\mu J^{\star} \vert \epsilon^{\star} (\mathbf{u}^{n+1}) \vert^2 \\
& - \int_{\Omega_0} \frac{\rho}{2\tau} J^{\star\star} \vert \mathbf{u}^{n+1} - \mathbf{u}^{n} \vert^2 + \int_{\Omega_0} \frac{\rho}{2} Div(J^{\star} F^{-1}_{\star} (\mathbf{u}^{\ast} - \mathbf{w}^{\ast\ast})) \vert \mathbf{u}^{n+1} \vert^2 \\
& - \int_{\Omega_0} \frac{\rho}{2} \alpha \frac{J^{n+1} - J^{n}}{\tau} \vert \mathbf{u}^{n+1} \vert^2 + \int_{\Omega_0} \frac{\rho}{2} Div\left( J^{\star}F^{-1}_{\star} (\beta \mathbf{u}^{\ast} - \alpha \mathbf{w}^{\ast\ast}) \right) \vert \mathbf{u}^{n+1} \vert^2
\end{aligned}
\end{equation}
Thus, for $\alpha=\beta=1$ and $\star\star = n$ the result follows.
\qed
\end{proof}
%\begin{remark}
%In the general case where Dirichlet and Neumann boundary conditions appear, $\Gamma_N \cup \Gamma_D = \partial \Omega_0$ with $\Gamma_N \neq \emptyset$, the same result also holds true with the addition of backflow stabilization term.
%\end{remark}
\section{Chorin-Temam schemes}
\label{sec:chorin_temam_schemes}
In the following, we describe a family of Chorin-Temam (CT) schemes for the iNSE-ALE problem, as we did for the monolithic case. %Such description keeps the freedom of choice for certain coefficients, which must be restricted to ensure unconditional stability.
Given $\mathbf{\widetilde V} $ a conforming space of $\mathbf{H}^1_0 (\Omega_0)$ and $ \widetilde Q$ a conforming space of $ L^2_0 (\Omega_0) \cap H^1(\Omega_0)$, $\tilde{\mathbf{u}}^{0} \in \mathbf{\widetilde V}$, for $n\geq0$:% the proposed two-step CT schemes reads
\begin{enumerate}[1.]
\item \textbf{Pressure-Projection Step $(\text{PPS})_{n}$}
%Seeks a pressure term allowing the projection of $\tilde{\mathbf{u}}^{n+1}$ in the space of divergence-free solutions, in the form:
Find $p^{n} \in \widetilde{Q}$ s.t.
\begin{equation}
\label{eq:chorin_temam_pps}
\int_{\Omega_0} \frac{\tau}{\rho} J^{\circ} Grad(p^{n}) F^{-1}_{\circ} : Grad(q) F^{-1}_{\circ} = - \int_{\Omega_0} Div\left( J^{\circ} F^{-1}_{\circ} \tilde{\mathbf{u}}^{n}\right) q \quad \forall q \in \widetilde{Q}
\end{equation}
\item \textbf{Fluid-Viscous Step $(\text{FVS})_{n+1}$} %Seeks a tentative velocity field $\tilde{\mathbf{u}}^{n+1}$ with pressure given explicitly:
Find $\tilde{\mathbf{u}}^{n+1} \in \mathbf{\widetilde{V}}$ s.t.
\begin{equation}
\label{eq:chorin_temam_fvs}
\begin{aligned}
\int_{\Omega_0} \rho J^{\star\star} \frac{\tilde{\mathbf{u}}^{n+1} - \tilde{\mathbf{u}}^n}{\tau} \cdot \mathbf{v} + \int_{\Omega_0} \rho J^{\star} Grad(\tilde{\mathbf{u}}^{n+1}) F_{\star}^{-1} (\tilde{\mathbf{u}}^{n} - \mathbf{w}^{\ast\ast}) \cdot \mathbf{v} +\int_{\Omega_0} J^{\star} 2 \mu \epsilon^{\star} (\tilde{\mathbf{u}}^{n+1}) : \epsilon^{\star} (\mathbf{v}) & \\
-\int_{\Omega_0} Div(J^{\circ \circ} F_{\circ \circ}^{-1} \mathbf{v}) p^n
+ \int_{\Omega_0} \frac{\rho}{2} \frac{J^{n+1} - J^{n}}{\tau} \tilde{\mathbf{u}}^{n+1} \cdot \mathbf{v} + \int_{\Omega_0} \frac{\rho}{2} Div\left( J^{\star} F_{\star}^{-1}(\tilde{\mathbf{u}}^{n} - \mathbf{w}^{\ast\ast})\right) \tilde{\mathbf{u}}^{n+1} \cdot \mathbf{v} & = 0 \quad \forall \mathbf{v} \in \mathbf{\widetilde{V}}
\end{aligned}
\end{equation}
% In particular, the corrected velocity is obtained through the update:
% \begin{equation}
% \begin{aligned}
% \mathbf{u}^{n+1} & = \tilde{\mathbf{u}}^{n+1} - \frac{\tau}{\rho} Grad(p^{n+1}) F^{-1}_{\circ + 1} \text{ in } \Omega_0
% \end{aligned}
% \end{equation}
\end{enumerate}
%\begin{remark}
%The CT scheme iNSE above, uses a linearization of the convective term $\ast = n$, easily extended to the case $\ast = n+1$.. In particular, the stabilization terms are active $\alpha = \beta = 1$, with freedom in the choice of $\circ\circ$. The objective in the following will be to show which values achieve unconditional stability in our formulation.
%\end{remark}
%\subsection{Time Stability}
The following energy estimate can be obtained under suitable conditions:
\begin{proposition}
\label{prop:energy_estimate_chorin_temam}
Under assumptions $\circ = \circ\circ = \star \star = n$, the solution to scheme \eqref{eq:chorin_temam_pps}-\eqref{eq:chorin_temam_fvs} is unconditionally stable with
%is unconditionally energy stable if $\circ = \circ \circ = \star\star = n$, without condition over $\star$. Moreover, the energy estimate is given by
\begin{equation}
\begin{aligned}
\int_{\Omega_0} \rho \frac{J^{n+1}}{2\tau} \vert \tilde{\mathbf{u}}^{n+1} \vert^2 - \int_{\Omega_0} \rho \frac{J^{n}}{2\tau} \vert \tilde{\mathbf{u}}^{n} \vert^2 \leq & - \int_{\Omega_0} J^{\star} 2\mu \vert \epsilon^{\star} (\tilde{\mathbf{u}}^{n+1}) \vert^2 - \int_{\Omega_0} J^{n} \frac{\tau}{2\rho} \vert Grad(p^n) F^{-1}_{n} \vert^2 .
\end{aligned}
\end{equation}
\end{proposition}
\begin{proof}
As standard in literature, let us take $\mathbf{v} = \tilde{\mathbf{u}}^{n+1}$ in $(\text{FVS})_{n+1}$, and $q = p^n$ in $(\text{PPS})_{n}$. Adding both equalities and rewriting expressions, it follows:
\begin{equation}
\begin{aligned}
\int_{\Omega_0} \rho \frac{J^{n+1}}{2\tau} \vert \tilde{\mathbf{u}}^{n+1} \vert^2 - \int_{\Omega_0} \rho \frac{J^{n}}{2\tau} \vert \tilde{\mathbf{u}}^{n} \vert^2 = & \int_{\Omega_0} \frac{\rho}{2\tau}(J^{n+1} - J^{\star\star}) \vert \tilde{\mathbf{u}}^{n+1} \vert^2 + \int_{\Omega_0} \frac{\rho}{2\tau}(J^{\star\star} - J^{n}) \vert \tilde{\mathbf{u}}^{n} \vert^2 \\
& - \int_{\Omega_0} \frac{\rho}{2\tau} J^{\star\star} \vert \tilde{\mathbf{u}}^{n+1} - \tilde{\mathbf{u}}^{n} \vert^2 - \int_{\Omega_0} J^{\star} 2\mu \vert \epsilon^{\star} (\tilde{\mathbf{u}}^{n+1}) \vert^2 \\
& + \int_{\Omega_0} Div\left( J^{\circ\circ} F^{-1}_{\circ\circ} (\tilde{\mathbf{u}}^{n+1} - \tilde{\mathbf{u}}^{n}) \right) p^n + \int_{\Omega_0} Div\left( (J^{\circ\circ}F^{-1}_{\circ\circ} - J^{\circ}F^{-1}_{\circ}) \tilde{\mathbf{u}}^n \right) p^n \\
& - \int_{\Omega_0} \frac{\tau}{\rho} J^{\circ} \vert F^{-T}_{\circ} Grad(p^n) \vert^2 - \int_{\Omega_0} \frac{\rho}{2\tau} (J^{n+1} - J^{n}) \vert \tilde{\mathbf{u}}^{n+1} \vert^2
\end{aligned}
\end{equation}
Bounding the first divergence term using integration by parts and Cauchy-Schwartz inequality, it follows
\begin{equation}
\begin{aligned}
\int_{\Omega_0} Div\left( J^{\circ\circ} F^{-1}_{\circ\circ} (\tilde{\mathbf{u}}^{n+1} - \tilde{\mathbf{u}}^{n}) \right) p^n & \leq \int_{\Omega_0} \frac{\rho}{2\tau} J^{\circ\circ} \vert \tilde{\mathbf{u}}^{n+1} - \tilde{\mathbf{u}}^{n} \vert^2 + \int_{\Omega_0} \frac{\tau}{2\rho} J^{\circ\circ} \vert F^{-T}_{\circ\circ} Grad(p^n) \vert^2
\end{aligned}
\end{equation}
Thus, the following energy estimate can be obtained:
\begin{equation}
\label{eq:energy_estimate_proof}
\begin{aligned}
\int_{\Omega_0} \rho \frac{J^{n+1}}{2\tau} \vert \tilde{\mathbf{u}}^{n+1} \vert^2 - \int_{\Omega_0} \rho \frac{J^{n}}{2\tau} \vert \tilde{\mathbf{u}}^{n} \vert^2 \leq & \int_{\Omega_0} \frac{\rho}{2\tau} (J^{n+1} - J^{\star\star}) \vert \tilde{\mathbf{u}}^{n+1} \vert^2 + \int_{\Omega_0} \frac{\rho}{2\tau} (J^{\star\star} - J^{n}) \vert \tilde{\mathbf{u}}^{n} \vert^2 \\
& - \int_{\Omega_0} \frac{\rho}{2\tau} J^{\star\star} \vert \tilde{\mathbf{u}}^{n+1} - \tilde{\mathbf{u}}^{n} \vert^2 - \int_{\Omega_0} J^{\star} 2\mu \vert \epsilon^{\star} (\tilde{\mathbf{u}}^{n+1}) \vert^2 \\
& + \int_{\Omega_0} \frac{\rho}{2\tau} J^{\circ\circ} \vert \tilde{\mathbf{u}}^{n+1} - \tilde{\mathbf{u}}^{n} \vert^2 + \int_{\Omega_0} \frac{\tau}{2\rho} J^{\circ\circ} \vert F^{-T}_{\circ\circ} Grad(p^n) \vert^2 \\
& + \int_{\Omega_0} Div\left( (J^{\circ\circ} F^{-1}_{\circ\circ} - J^{\circ}F^{-1}_{\circ}) \tilde{\mathbf{u}}^{n}\right) p^n - \int_{\Omega_0} \frac{\tau}{\rho} J^{\circ} \vert F^{-T}_{\circ} Grad(p^n) \vert^2 \\
& - \int_{\Omega_0} \frac{\rho}{2\tau} (J^{n+1} - J^{n}) \vert \tilde{\mathbf{u}}^{n+1} \vert^2
\end{aligned}
\end{equation}
From estimate \eqref{eq:energy_estimate_proof} it follows that whenever $\circ = \circ \circ = \star \star = n$ unconditional energy stability is attained, where $\star$ remains free of choice.
\qed
\end{proof}
%\begin{remark}
%In the general case where $\Gamma_N \cup \Gamma_D = \partial \Omega_0$ with $\Gamma_N \neq \emptyset$ a similar bound is attained with the addition of backflow stabilization term, where the same procedure as above holds for all terms except the one on $\Gamma_N$.
%\end{remark}
\section{Numerical examples}
\label{sec:numerical_examples}
We consider a rectangular domain with opposite vertices $\{ (0, -1), (6, 1) \}$ where the iNSE-ALE formulation \eqref{eq:continuous_formulation} will be simulated over the interval $[0, 2] \, [s]$ with non-zero initial condition of the form $ \mathbf{u}(0) := \gamma (1 - \mathbf{X}_1^2) \mathbf{X}_0 (6 - \mathbf{X}_0), \, \gamma = 0.001$.
The domain is deformed using $\mathcal{X}(\mathbf{X}, t) := \big( (1 + 0.9 sin(8 \pi t)) \mathbf{X}_0,\, \mathbf{X}_1 \big)$. %, i.e. an oscillation with initial expansion and frequency of $4 \pi$.
Discretization setup for Formulation \eqref{eq:discretized_monolithic_formulation} and \eqref{eq:chorin_temam_pps}-\eqref{eq:chorin_temam_fvs} is done choosing a time step $\tau = 0.01$ and space triangulation with elements diameter $h \approx 0.01 $, implemented through FEniCS \cite{FEniCS2015} using Python for interface and postprocessing.
To exemplify the theoretical results from previous sections, four schemes are taken into account. Monolitic (M) Formulation \eqref{eq:discretized_monolithic_formulation} is taken with linearized convective term and implicit treatment, i.e., $(\star, \ast, \ast\ast) = (n+1, n, n+1)$ where for $\star\star$ we consider two choices, denoted $\text{M}\, \star\star = n$ and $\text{M}\, \star\star = n+1$.
For both cases the space discretization is carried out with $ \mathbf{V}/Q = [\mathbb{P}_2]^d/\mathbb{P}_1$ Lagrange finite elements. Similarly, Chorin-Temam (CT) scheme \eqref{eq:chorin_temam_fvs}-\eqref{eq:chorin_temam_pps} is taken with linearized convective term and implicit treatment, i.e. $(\star, \ast\ast, \circ, \circ\circ) = (n+1, n+1, n, n)$ with $\star\star$ as before, denoting each scheme by $\text{CT}\, \star\star = n$ and $\text{CT}\, \star\star = n+1$ with space discretization done through $\mathbf{\widetilde{V}}/\widetilde{Q} = [\mathbb{P}_1]^d/\mathbb{P}_1$ elements. In all cases $\alpha=\beta=1$.
The results are assessed using time-dependent normalized parameters $ \hat{\delta}_{\text{M}}:= \delta_{\text{M}}/E_{st}^{\star}, \hat{\delta}_{\text{CT}}:= \delta_{\text{CT}}/E_{st}^{\star}$ defined as: %, where $E_{st}^{\star}$ denote the \CB{\mod{strain energy}{viscous dissipation}} and $\delta_{M}, \delta_{CT}$ residual errors defined as:
\begin{equation}
\label{eq:energy_error}
\begin{aligned}
\delta_{M}^{n+1} &:= D^{n+1} + E_{st}^{\star} + \int_{\Omega_0} \frac{\rho J^{\star\star}}{2\tau} \vert \mathbf{u}^{n+1} - \mathbf{u}^{n} \vert^2, \quad \delta_{CT}^{n+1} := D^{n+1} + E^{\star}_{st} + \int_{\Omega_0} \frac{\tau J^{\circ}}{2 \rho} \vert F^{-T}_{\circ} Grad(p^n) \vert^2 \\
D^{n+1} &:= \int_{\Omega_0} \frac{\rho}{2\tau} \left( J^{n+1} \vert \mathbf{u}^{n+1} \vert^2 - J^n \vert \mathbf{u}^n \vert^2 \right), \quad E^{\star}_{st} = \int_{\Omega_0} 2 \mu J^{\star} \vert \epsilon^{\star} (\mathbf{u}^{n+1}) \vert^2
\end{aligned}
\end{equation}
Figure \ref{fig:delta_hat_oscs} shows $\hat{\delta}_{\text{M}}, \hat{\delta}_{\text{CT}}$ values for each tested scheme. Propositions \ref{prop:energy_estimate_monolithic} and \ref{prop:energy_estimate_chorin_temam} are confirmed since $\hat{\delta}_{\text{M}}=0$ and $\hat{\delta}_{\text{CT}} \leq 0$ if $\star\star = n$. For $\star \star = n+1$, peaks appearing throughout the simulation are defined by the sign change of domain velocity, i.e. in the change from expansion to contraction.
Importantly, the spurious numerical energy rate related to discretization of the GCL condition appear to be positive in expansion, therefore being a potential source of numerical instabilities.
%Moreover the expected energy decay from Propositions \ref{prop:energy_estimate_monolithic}, \ref{prop:energy_estimate_chorin_temam} for the choice $\star \star = n$ in $\text{M, CT}$ schemes is obtained.
%\CB{Recall that $\hat{\delta}_{\text{M}}, \hat{\delta}_{\text{CT}}$ in \eqref{eq:energy_error} are taken in such form since on scheme $\text{M}\, \star\star = n$ it follows $\hat{\delta}_{MT} = 0$ and similarly on $\text{CT}\, \star\star = n$ it follows $\hat{\delta}_{\text{CT}} \sim 0$, defining the base case.}
\begin{figure}[!hbtp]
\centering
\includegraphics[width=\textwidth]{figs/Comparison_Delta_Hat_Value_GCL_True_solver_LU.png}
\caption{Summary of the numerical experiment in terms of energy balance. Left: Monolithic residual error values $\hat{\delta}_{\text{M}}$;
Right: Chorin-Temam residual error values $\hat{\delta}_{\text{CT}}$. %Both cases are simulated with deformation mapping $\mathcal{A}$ defined previously and shown on the interval $[0, 1] \, [s]$.}
}
\label{fig:delta_hat_oscs}
\end{figure}
\section{Conclusion}
Reported first order time discretization schemes for the iNSE-ALE have been reviewed, theoretically analyzed and numerically assessed. For the monolithic case, two schemes lead to well-posed energy-stable problems whenever $\alpha=\beta=1$ with $\star \star = n$ as studied by \cite{LeTallec2001, Lozovskiy2018, smaldone2014, Wang2020}.
To the best of the authors knowledge, the unconditionally stable Chorin-Temam scheme derived in this work has not been reported yet. %, and moreover the numerical experiments studied here validate the energy stable propositions.
%\appendix
%\printcredits
%% Loading bibliography style file
%\bibliographystyle{model1-num-names}
\bibliographystyle{cas-model2-names}
% Loading bibliography database
\bibliography{cas-refs}
%\vskip3pt
%\bio{}
%Author biography without author photo.
%Author biography. Author biography. Author biography.
%Author biography. Author biography. Author biography.
%Author biography. Author biography. Author biography.
%Author biography. Author biography. Author biography.
%Author biography. Author biography. Author biography.
%Author biography. Author biography. Author biography.
%Author biography. Author biography. Author biography.
%Author biography. Author biography. Author biography.
%Author biography. Author biography. Author biography.
%\endbio
%
%%\bio{figs/pic1}
%Author biography with author photo.
%Author biography. Author biography. Author biography.
%Author biography. Author biography. Author biography.
%Author biography. Author biography. Author biography.
%Author biography. Author biography. Author biography.
%Author biography. Author biography. Author biography.
%Author biography. Author biography. Author biography.
%Author biography. Author biography. Author biography.
%Author biography. Author biography. Author biography.
%Author biography. Author biography. Author biography.
%\endbio
%\bio{figs/pic1}
%Author biography with author photo.
%Author biography. Author biography. Author biography.
%Author biography. Author biography. Author biography.
%Author biography. Author biography. Author biography.
%Author biography. Author biography. Author biography.
%\endbio
\end{document}
Reference in New Issue View Git Blame Copy Permalink