documentation.md 107 KB
<body>
<div>
\(
\definecolor{florange}{rgb}{0.832 ,0.367 ,0}
\definecolor{flblue}{rgb}{0 ,0.445 ,0.695}
\definecolor{flgreen}{rgb}{0 ,0.617 ,0.449}
\definecolor{flred}{rgb}{1 ,0 ,0 }
\)
</div>
Documentation
This document contains extracts from the theses and articles listed in References.
Abbreviations
- BC: boundary conditions
- LS: levelset (distance to a fluid-fluid interface or fluid-solid interface i.e. a wall)
- p-grid: scalar grid, see staggered grid
- u-grid: staggered grid (with regard to p-grid, in x-axis), used to discretize the x-component of the velocity
- v-grid: staggered grid (with regard to p-grid, in y-axis), used to discretize the y-component of the velocity
- ghost cells: cells defined outside of the computational domain so that the same stencils may be used near the wall (see
allocate_ghost_matrices_2,init_ghost_neumann_2,IIOE_normal_indices_2!; for MPI parallelization, ghost cells may also refer to cells outside the computational domain of the current processor, not implemented here)
Algorithm, main function
run_forward!
!!! todo "change order scalar transport electrical potential"
Algorithm
- Initialize simulation parameters
- Allocations
- Initialisation of bulk and interfacial values
- Set matrices/operators
- Solve heat equation
- Poisson equation (electrical potential)
- Reevaluate the current i_current for boundary conditions
- Scalar transport: update boundary conditions from electrical current
- Scalar transport: solve
- Phase change
- Advection of the levelset
- Pressure-velocity coupling (Navier-Stokes)
!!! info "Viewing the main parts of the algorithm in VScode"
Comments starting with
#region and ending with
#endregion
You can visualise the regions correponding to this algorithm more easily with tools like "Region marker" in VScode.
!!! danger "update LS"
See [Updating the operator from Levelset](@ref)
Parameters for Flower.jl
Numerical
Definitions
Convention and orientation of the bubble interface
The signed distance is positive in the liquid and negative in the bubble.
The normal defined by LS.\alpha is oriented towards the liquid, so we have to take the opposite to define the outward normal for manufactured solutions for example (so that it points towards the interior of the bubble when solving in the liquid).
Multiple levelsets
In Flower, multiple levelsets may be defined. The last levelset is the combination of all levelsets and is used to define the cell and its centroid.
!!! todo "TODO"
Logic center refers to indices in the matrices.
"
In this work, a levelset function \phi \cite{Sethian1999} is defined on the computational domain \Omega to map the locus of one of its iso-levels (\left \{ \left . x \in \Omega \right | \phi \left ( x, t \right ) = \phi _ 0 \right \}) to an interface \Gamma \left ( t \right ) that separates two non-overlapping domains, \Omega _ 1 \left ( t \right ) and \Omega _ 2 \left ( t \right ), each occupied by a different phase. The value \phi is defined as the signed distance to the interface,
\begin{equation}
\phi \left ( x, t \right ) = \left \{ \begin{aligned}
-d \left ( x, \Gamma \left ( t \right ) \right ),& \ x \in \Omega _ 1 \left ( t \right ) \\
0,& \ x \in \Gamma \left ( t \right ) \\
d \left ( x, \Gamma \left ( t \right ) \right ),& \ x \in \Omega _ 2 \left ( t \right )
\end{aligned} \right .,
\end{equation}
where d \left ( x, \Gamma \left ( t \right ) \right ) denotes the minimal distance between the point x and the interface \Gamma \left ( t \right ),
"
<figure>
<a name="levelset_doc"></a>
<img src="./assets/levelset_doc.svg" alt="Levelset" title="Levelset defining an interface $\Gamma$ separating two phases $\Omega _ 1$ and $\Omega _ 2$">
<figcaption>Levelset defining an interface $\Gamma$ separating two phases $\Omega _ 1$ and $\Omega _ 2$ </figcaption>
</figure>
Geometric moments
The mesh is assumed to be rectilinear with n _ x points along x and n _ y along y. For wet and partially wet cells numbered i along x and j along y, the coordinates of the fluid center of mass are denoted x ^ \omega _ {i, j} \quad \mathrm{and} \quad y ^ \omega _ {i, j} and represent the components of two vectors of size n _ x n _ y (denoted x ^ \omega and y ^ \omega). Likewise, the vectors x ^ \gamma and y ^ \gamma, of size n _ x n _ y, store the coordinates of the center of mass of the boundary of partially wet cells (referred to as boundary cells).
In the following, discrete values can be cell-, face- or node-centered, therefore a numbering convention must be adopted for the indices \left ( i, j \right ). Fig.~\ref{fig:numbering} showcases the convention used in this work for both cell and face-centered quantities (the \left ( i, j \right ) node-centered quantities is located at the bottom-left of the cell).
<figure>
<a name="numbering_doc"></a>
<img src="./assets/numbering_doc.svg" alt="Numbering convention" title="Numbering convention">
<figcaption> Numbering convention </figcaption>
</figure>
"
Cut-cell methods are firmly grounded in the Finite Volume Method, which defines the primary discrete variables as cell-wise averages over mesh elements. The design of the Finite Volume operators is then based on the application of the Divergence theorem. For example, given a scalar field p, this theorem states that in a Cartesian coordinate system, the x component of the gradient q \equiv \nabla p averaged over a cell \Omega may be computed as
\begin{equation}
\Omega q _ x = \int _ \Omega \frac{\partial p}{\partial x} V = \oint _ {\partial \Omega} p e _ x \cdot S
\label{eq:Stokes}
\end{equation}
where S denotes the outward-pointing surface element and e_x the unit vector along the x direction.
For the sake of presentation, the case displayed in Fig.~\ref{fig:moments} is considered where $\Omega \equiv \mathcal{V} _ {i,j}$ consists of the intersection of a phase domain and a computational cell (a right hexahedron). The contour $\partial \Omega$ then consists of the union of the four planar faces $\mathcal{A} _ {x, i, j}$, $\mathcal{A} _ {x, i+1, j}$, $\mathcal{A} _ {y, i, j}$ and $\mathcal{A} _ {y, i, j+1}$ as well as the boundary surface $\Gamma _ {i, j}$. A piece-wise linear approximation of $\Gamma _ {i, j}$, denoted $\widetilde{\Gamma} _ {i, j}$ and of unit normal $\left ( n _ {x, i, j}, n _ {y, i, j} \right )$, can be defined by applying Eq.~\eqref{eq:Stokes} to $\mathcal{V} _ {i, j}$ with $p = 1$, yielding
\int _ {\widetilde{\Omega}} \frac{\partial 1}{\partial x} \mathrm{d} V = \mathcal{A} _ {x, i+1, j} - \mathcal{A} _ {x, i, j} + n _ {x, i, j} \widetilde{\Gamma} _ {i, j} = 0
which highlights the existence of a fundamental relation \begin{equation}
\mathcal{A} _ {x, i+1, j} - \mathcal{A} _ {x, i, j} = -n _ {x, i, j} \widetilde{\Gamma} _ {i, j}
\label{eq:approximate}
\end{equation} sometimes referred to as a Geometric Surface Conservation Law. The same can be done in the $y$-direction in order to obtain a relation between $\mathcal{A} _ {y, i, j+1}, \mathcal{A} _ {y, i, j}, n _ {y, i, j}$ and $\widetilde{\Gamma} _ {i, j}$. Note that to simplify the exposition, the notations ($\mathcal{A} _ {x, i, j}$, $\mathcal{V} _ {i, j}$...) are used to denote both a region in space and its measure (length or area, in 2D). "
Domains
\Gamma is the interface and \omega is the wet volume of the cell without the interface: \omega = \Omega \setminus \Gamma.
Centroids
In the code, the cell centroid for the p-grid (scalar grid, here noted gp) corresponding to the first levelset (LS[1]) is x^\omega = (x_{\text{centroid}},y_{\text{centroid}}):
x_centroid = gp.x .+ getproperty.(gp.LS[1].geoL.centroid, :x) .* gp.dx
y_centroid = gp.y .+ getproperty.(gp.LS[1].geoL.centroid, :y) .* gp.dy
For the v-grid:
x_centroid = gv.x .+ getproperty.(gv.LS[1].geoL.centroid, :x) .* gv.dx
y_centroid = gv.y .+ getproperty.(gv.LS[1].geoL.centroid, :y) .* gv.dy
<figure>
<a name="numbering_doc"></a>
<img src="./assets/staggered_coupled.svg" alt="Staggered grids" title="Staggered grids">
<figcaption> Staggered grids </figcaption>
</figure>
julia +1.10.5 --project=../Flower.jl --threads=1 ../Flower.jl/examples/convergence.jl ../Flower.jl/examples/convergence_diffusion_constant_conductivity_bubble_wall_perio.yml
The interface centroid x^\gamma is defined as "the mid point of the segment crossing the cell which will be used in the computation of the Stefan condition." (Fullana 2022).
x_bc = gp.x .+ getproperty.(gp.LS[1].mid_point, :x) .* gp.dx
y_bc = gp.y .+ getproperty.(gp.LS[1].mid_point, :y) .* gp.dy
!!! info
``(gp.x, gp.y)`` is the cell center, not the centroid. In the current post-processings the cell center is used if not said otherwise.
!!! info
The divergence ``\mathrm{div}(0,\mathrm{grad})`` is a bulk field, it is located at the cell centroids.
H
The boundary field is always located at the borders of the p-grid. So for scalars, the distance is always half a cell in both directions because the cell centroid is in the cell center. For the u-grid, since it's shifted half a cell in x, then from the wet cell centroid to the boundary there's a quarter of a cell at the left and right borders. And the same for the v-grid, it's displaced half a cell in y, so from the wet cell centroid to the boundary there's a quarter of a cell of distance.
!!! todo "TODO"
Without interface, with a constant mesh spacing in every direction, non-dimensionalized by the mesh spacing H should be equal to:
| Border | p-grid | u-grid | v-grid |
|---|---|---|---|
| left | 0.5 | 0.25 | 0.5 |
| bottom | 0.5 | 0.5 | 0.25 |
| right | 0.5 | 0.25 | 0.5 |
| top | 0.5 | 0.5 | 0.25 |
Capacities
The face capacities for the liquid phase of levelset n° iLS for grid are stored in grid.LS[iLS].geoL.cap[II,1:4]:
grid.LS[iLS].geoL.cap[II,1] # left, also known as A1 in the code
grid.LS[iLS].geoL.cap[II,2] # bottom, also known as A2 in the code
grid.LS[iLS].geoL.cap[II,3] # right, also known as A3 in the code
grid.LS[iLS].geoL.cap[II,4] # top, also known as A4 in the code
The cell-centered volume is:
grid.LS[iLS].geoL.cap[II,5] # cell-centered volume
The cell-centered heights are:
grid.LS[iLS].geoL.cap[II,6] # cell-centered, also known as B1
grid.LS[iLS].geoL.cap[II,7] # B2
Others:
grid.LS[iLS].geoL.cap[II,8] # W1
grid.LS[iLS].geoL.cap[II,9] # W2
grid.LS[iLS].geoL.cap[II,10] # W3
grid.LS[iLS].geoL.cap[II,10] # W4
!!! danger "half cells for u and v cap = 1/2 ?"
!!! info
In the code, B1 B2 are placed at the cell centroid, not the cell center.
iMx is the inverse of a volume. The border cells for u and v have half the normal cell size because the borders of the domain pass through the cell centers, not the cell faces like they do for p.
Example for 1D (equations)
This section explains the equations and their implementation for a simple geometry with one levelset describing a square bubble and a wall at the limits of the computational domain.
Operators in a cell of the mesh
Let us consider a square bubble in contact with a wall. For simplicity, the faces of a regular cartesian p-grid are displayed in black. The wall (in red) and bubble (in green) interfaces are shifted so that two boundary conditions can be imposed simultaneously ( see Simultaneous boundary conditions).
<figure>
<img src="./assets/mesh_square_bubble.svg" alt="Square bubble in contact with wall" title="Square bubble in contact with wall">
<figcaption>Square bubble in contact with wall</figcaption>
</figure>
Each boundary condition is solved in a separate cell.
<figure>
<img src="./assets/mesh_square_bubble_2.svg" alt="Square bubble in contact with wall" title="Square bubble in contact with wall">
<figcaption>Square bubble in contact with wall</figcaption>
</figure>
The scalar control volume $\mathcal{V}$, $\mathcal{B}$ capacities (height of wetted part in x and y directions) are associated with the centroid of the current phase (e.g. liquid): $x^\omega$ <span class="hover_img"> <a href="#centroidcontrolvolume"> (see figure) <span><img src="./assets/Vt_mesh_square_bubble_wall_1.svg" alt="image" height="800" /></span></a></span>.
The centroid is required even if it does not intervene in the discretization.
The LibGEOS.jl library, based on GEOS, is used for the geometrical computations, the VOFI library could also be used: (bna et al., 2016).
The scalar control volume \mathcal{W}, \mathcal{A} capacities (height of wetted part in x and y directions) are associated with the center of the faces of the cells.
The first kind capacities are x^\omega (cell), V (cell), A^x (face), the second kind capacities are B^x (cell) and \mathcal{W}^x (face).
The interface of the left wall is shifted at o^- and the bubble interface is also shifted so that two different boundary conditions (BC) may be imposed simultaneously.
When a volume is null, the associated capacity \mathcal{B} must be null as well.
At wall
The scalar control volume of cell (i=1,j) and its associated centroid are represented below:
<figure>
<a name="centroidcontrolvolume"></a>
<img src="./assets/V_mesh_square_bubble_wall_1.svg" alt="Centroid control volumes at i = 1" title="Centroid control volumes at i = 1">
<figcaption>Centroid control volumes at i = 1</figcaption>
</figure>
The control volume \mathcal{W}_{x,i=1} associated with the face is represented below:
<figure>
<a name="facecontrolvolume"></a>
<img src="./assets/Wx_mesh_square_bubble_wall_2.svg" alt="Face control volume " title="Face control volume">
<figcaption>Face control volume</figcaption>
</figure>
The values of the capacities are detailed here:
<figure>
<a name="capacities"></a>
<img src="./assets/Vt_mesh_square_bubble_wall_3.svg" alt="Centroid control volumes at i=0 and i = 1" title="Centroid control volumes at i=0 and i = 1">
<figcaption>Centroid control volumes at i=0 and i = 1</figcaption>
</figure>
\begin{array}{cccccccc}
\hline
Variable & Face ~-\frac 12 & Bulk ~ i = 0 & \color{flred}{0^-} & Face ~\frac 12 & Bulk ~ i = 1 & Face ~\frac 32 & Bulk ~ i = 2 \\
\hline
V & & 0 & \color{flred}{|} & & h_xh_y & & h_xh_y\\
A_x & 0 & &\color{flred}{|} &h_y & & h_y & \\
B_x & & 0 &\color{flred}{|} & & h_y & & h_y \\
W_x & 0 & &\color{flred}{|} &\frac{h_xh_y}{2} & & h_xh_y &\\
\hline
\end{array}
The location of the centroid of \mathcal{W}_x is not used in the current discretization.
Gradient
x-component
At cell (i,j), with \Gamma the interface and c the scalar, the x-component of the cell-averaged gradient is:
\begin{aligned}
{\mathrm{grad}(c^ω, c^{γ})}_{i +\frac 12} &= \mathcal{W}_{i +\frac 12}^\dagger (B_{i+1}c_{i+1}^\omega-B_{i}c_{i}^\omega \\
&+ (A_{i+\frac 12} - B_{i+1})c_{i+1}^\gamma + (B_{i}-A_{i+\frac 12})c_{i}^\gamma)
\end{aligned}
In the case of
<span class="hover_img"> <a href="#capacities">this figure <span><img src="./assets/Vt_mesh_square_bubble_wall_3.svg" alt="image" height="800" /></span></a></span>, the above expression is the same as with finite differences:
\begin{aligned}
{\mathrm{grad}(c^ω, c^{γ})}_{0 +\frac 12} &= \mathcal{W}_{\frac 12}^\dagger (B_{1}c_{1}^\omega-B_{0}c_{0}^\omega \\
&+ (A_{\frac 12} - B_{1})c_{1}^\gamma + (B_{0}-A_{\frac 12})c_{0}^\gamma) \\
&= \frac{1}{\frac{h_x}{2}h_y} (h_yc_{1}^\omega- 0c_{0}^\omega \\
&+ (h_y - h_y)c_{1}^\gamma + (0-h_y)c_{0}^\gamma) \\
&=\frac{c_1^\omega-c_0^\gamma}{\frac{h_x}{2}}
\end{aligned}
Away from the wall and without interface, the formula is:
\begin{aligned}
{\mathrm{grad}(c^ω, c^{γ})}_{i +\frac 12} =\frac{c_{i+1}^\omega-c_i^\omega}{h_x}
\end{aligned}
y-component
The y-component is given below:
\begin{aligned}
{\mathrm{grad}(c^ω, c^{γ})}_{i +\frac 12} &= \mathcal{W}_{i +\frac 12}^\dagger (B_{i+1}^yc_{j+1}^\omega-B_{i}^yc_{j}^\omega \\
&+ (A_{j+\frac 12} - B_{j+1})c_{i+1}^\gamma + (B_{j}-A_{j+\frac 12})c_{j}^\gamma)
\end{aligned}
Since \Gamma is shifted upwards in the y-direction A_{\frac 12}=0:
\begin{aligned}
{\mathrm{grad}(c^ω, c^{γ})}_{0 +\frac 12} &= \mathcal{W}_{\frac 12}^\dagger (B_{1}c_{1}^\omega-B_{0}c_{0}^\omega \\
&+ (A_{\frac 12} - B_{1})c_{1}^\gamma + (B_{0}-A_{\frac 12})c_{0}^\gamma) \\
&= \frac{1}{\frac{h_x}{2}h_y} (h_x c_{1}^\omega- 0c_{0}^\omega \\
&+ (0 - h_x)c_{1}^\gamma + (0-0)c_{0}^\gamma) \\
&=\frac{c_i^\omega-c_0^\gamma}{\frac{h_y}{2}}
\end{aligned}
Right wall
For the right wall, with nx the number of cells in the x direction:
!!! todo "TODO"
finish
\begin{aligned}
{\mathrm{grad}(c^ω, c^{γ})}_{nx +\frac 12} &= \mathcal{W}_{\frac 12}^\dagger (B_{nx+1}c_{nx+1}^\omega-B_{nx}c_{nx}^\omega \\
&+ (A_{nx+\frac 12} - B_{nx+1})c_{1}^\gamma + (B_{nx}-A_{nx+\frac 12})c_{nx}^\gamma) \\
&= \frac{1}{\frac{h_x}{2}h_y} (0 c_{nx+1}^\omega- h_yc_{nx}^\omega \\
&+ (h_y - 0)c_{nx+1}^\gamma + (h_y-h_y)c_{nx}^\gamma) \\
&=\frac{-c_1^\omega+c_{nx+1}^\gamma}{\frac{h_x}{2}}
\end{aligned}
Interface below
!!! todo "TODO"
finish and check
Since \Gamma is shifted downwards in the y-direction A_{\frac 12}=0:
\begin{aligned}
{\mathrm{grad}(c^ω, c^{γ})}_{} &= \mathcal{W}_{\frac 12}^\dagger (B_{}c_{1}^\omega-B_{0}c_{0}^\omega \\
&+ (A_{\frac 12} - B_{1})c_{1}^\gamma + (B_{0}-A_{\frac 12})c_{0}^\gamma) \\
&= \frac{1}{\frac{h_x}{2}h_y} (0 c_{1}^\omega- h_x c_{0}^\omega \\
&+ (0 - 0)c_{1}^\gamma + (h_x-0)c_{0}^\gamma) \\
&=\frac{-c_i^\omega+c_0^\gamma}{\frac{h_y}{2}}
\end{aligned}
Divergence
At cell (i,j), with \Gamma the interface
<a name="divergencetag"></a>
\begin{equation}
\label{divergence}
\begin{aligned}
{\mathrm{div}(q^ω, q^{γ})}_i &= \color{florange}{A^x_{i+\frac{1}{2}} q_{i+\frac{1}{2}}^\omega-A^x_{i-\frac{1}{2}} q_{i-\frac{1}{2}}^\omega} \\
&+\color{flblue}{(B_i^x -A^x_{i+\frac{1}{2}}) q_{i+\frac{1}{2}}^\gamma +(A^x_{i-\frac{1}{2}}-B_i^x ) q_{i-\frac{1}{2}}^\gamma } \\
&+\color{flgreen}{A^y_{j+\frac{1}{2}} q_{j+\frac{1}{2}}^\omega -A^y_{j-\frac{1}{2}} q_{j-\frac{1}{2}}^\omega } \\
&+\color{flblue}{(B_j^y -A^y_{j+\frac{1}{2}}) q_{j+\frac{1}{2}}^\gamma +(A^y_{j-\frac{1}{2}}-B_j^y ) q_{j-\frac{1}{2}}^\gamma } \\
\end{aligned}
\end{equation}
\begin{aligned}
\int_\Omega \nabla \cdot q dV &= \oint q \cdot n dS \\
& = \color{florange}{\oint_{i+\frac 12} q \cdot (+e_x) dS} + \color{florange}{\oint_{i-\frac 12} q \cdot (-e_x) dS} \\
& + \color{flgreen}{\oint_{j+\frac 12} q \cdot (+e_y) dS} + \color{flgreen}{\oint_{j-\frac 12} q \cdot (-e_y) dS} \\
& + \color{flblue}{\oint_\Gamma q \cdot n dS}
\end{aligned}
In the case of <span class="hover_img">
<a href="#capacities">this figure <span><img src="./assets/Vt_mesh_square_bubble_wall_3.svg" alt="image" height="800" /></span></a>
</span>, the x-contribution from the wall (from the above expression) is:
<a name="tagdivleftwall"></a>
\color{flblue}{(0-h_y) (\frac{c_i^\omega-c_0^\gamma}{\frac{h_x}{2}}) + (0-0) q_{i-\frac{1}{2}}^\gamma }
Since
<a name="tagdivbubbleinterface"></a>
The y-contribution for the green bubble is (by denoting j=1 the index of the cell cut by the interface, in the phase under consideration, i.e. above the bubble here):
\color{flblue}{(h_x-h_x) q_{\frac 32}^y + (0-h_x) q_{\frac 12}^y=-h_x (\frac{c_1^\omega-c_1^\gamma}{\frac{h_y}{2}}) }
The two former equations give an example of discretization of the boundary conditions at the left wall and at the interface.
Right wall
At nx+1
\color{flblue}{(0-0) (\frac{c_i^\omega-c_0^\gamma}{\frac{h_x}{2}}) + (h_y-0) q_{i-\frac{1}{2}}^\gamma }
\color{flblue}{ h_y q_{i-\frac{1}{2}}^\gamma \frac{-c_1^\omega+c_{nx+1}^\gamma}{\frac{h_x}{2}}}
Since
<a name="tagdivbubbleinterface"></a>
The y-contribution for the green bubble (under the bubble) is:
\color{flblue}{(h_x-0) q_{i+\frac 12}^y + (h_x-h_x) q_{i-\frac 12}^y=h_x (\frac{-c_1^\omega+c_1^\gamma}{\frac{h_y}{2}}) }
If h_x=h_y, we have the following coefficients on the rows of the interfacial values: +2 for interfacial, -2 for bulk
!!! danger "Check sign"
Coefficients in a simple configuration
<a name="Coefficients-in-a-simple-configuration"></a>
<table class="styled-table">
<thead>
<tr>
<th> BC </th>
<th>Border </th>
<th>Bulk </th>
</tr>
</thead>
<tbody>
<tr>
<th>Left</th>
<td> +2</td>
<td> -2</td>
</tr>
<tr>
<th>Right</th>
<td> +2</td>
<td> -2</td>
</tr>
<tr>
<th>Bottom</th>
<td> +2</td>
<td> -2</td>
</tr>
<tr>
<th>Top</th>
<td> +2</td>
<td> -2</td>
</tr>
</tbody>
</table>
Example away from interfaces
Away from the wall, without an interface, with a regular mesh of constant spacings h_x and h_y is:
\begin{aligned}
{\mathrm{div}(q^ω, q^{γ})}_i &= \color{florange}{A^x_{i+\frac{1}{2}} q_{i+\frac{1}{2}}^\omega-A^x_{i-\frac{1}{2}} q_{i-\frac{1}{2}}^\omega} \\
&+\color{flblue}{(B_i^x -A^x_{i+\frac{1}{2}}) q_{i+\frac{1}{2}}^\gamma +(A^x_{i-\frac{1}{2}}-B_i^x ) q_{i-\frac{1}{2}}^\gamma } \\
&+\color{flgreen}{A^y_{j+\frac{1}{2}} q_{j+\frac{1}{2}}^\omega -A^y_{j-\frac{1}{2}} q_{j-\frac{1}{2}}^\omega } \\
&+\color{flblue}{(B_j^y -A^y_{j+\frac{1}{2}}) q_{j+\frac{1}{2}}^\gamma +(A^y_{j-\frac{1}{2}}-B_j^y ) q_{j-\frac{1}{2}}^\gamma } \\
&= \color{flblue}{h_y (\frac{c_{i+1}^\omega-c_i^\omega}{h_x}) - h_y (\frac{c_{i}^\omega-c_{i-1}^\omega}{h_x})} \\
&+ \color{flgreen}{h_x (\frac{c_{j+1}^\omega-c_j^\omega}{h_y}) - h_x (\frac{c_{j}^\omega-c_{j-1}^\omega}{h_y})}
\end{aligned}
With h_x = h_y, this simplifies to:
{\mathrm{div}(q^ω, q^{γ})}_{i,j} = -4 c_{i,j} + c_{i-1,j} + c_{i+1,j} + c_{i,j-1} + c_{i,j+1}
Approximations inside a cell
The following hypothesis is made, q^\omega=q^\gamma=q, which means that the cell-averaged gradient (corresponding to the control volume ..., associated with the centroid ...), is associated at the middle of the part of interface, as illustrated by the black arrows below.
<figure>
<img src="./assets/Wx_mesh_square_bubble_wall_4.svg" alt="Face control volume " title="Face control volume">
<figcaption>Face control volume</figcaption>
</figure>
Poisson equation
!!! todo "Old implementation"
In [`set_poisson`](@ref) the system is `` + \nabla \cdot \nabla p = f`` :
In solve_poisson system is + \nabla \cdot \nabla p = f :
\begin{cases}
+ \mathrm{div} (q^\omega, q^\gamma ) &= V f^\omega\\
I_a I p^\gamma + I_b \mathrm{div} (0, q^\omega) &= I g^\omega\\
q^\omega &= \mathrm{grad} ( p^\omega, p^\gamma ) \\
q^\gamma &= q^\omega
\end{cases}
Second line: rechecking
!!! danger "Question"
since div in bulk and BC, and BC part of div is of opposite sign, why not opposite signs in matrix ?
full divergence
<a name="tagfulldivergence"></a>
\begin{equation}
\mathrm{div}(q^ω, q^{γ} ) = − (G^⊤ + H^{Γ1,⊤} + H^{Γ2,⊤}) q^ω + H^{Γ1,⊤}q^{γ1} + H^{Γ2,⊤}q^{γ2}
\end{equation}
<a name="grad_concise"></a>
\mathrm{grad}(q^ω, q^{γ} ) = W ^ \dagger \left ( G p ^ \omega + H p ^ \gamma \right )
\begin{equation}
\left [ \begin{array}{>{\centering\arraybackslash$} p{2.0cm} <{$} >{\centering\arraybackslash$} p{3.2cm} <{$}}
-G ^ \top W ^ \dagger G & -G ^ \top W ^ \dagger H \\
I _ b (H ^ \top W ^ \dagger G) & I _ b (H ^ \top W ^ \dagger H) + I _ a I _ \Gamma
\end{array} \right ] \left [ \begin{array}{c}
p ^ \omega \\
p ^ \gamma
\end{array} \right ] \simeq \left [ \begin{array}{c}
V f ^ \omega \\
I _ \Gamma g ^ \gamma
\end{array} \right ],
\end{equation}
Interface length
The interface length is defined as:
\chi = \sqrt{{(A_{\frac 12}^x-B_0^x)}^2+{(B_0^x-A_{-\frac 12}^x)}^2 + {(A_{\frac 12}^y-B_0^y)}^2+{(B_0^y-A_{-\frac 12}^y)}^2}
In Rodriguez 2024, it is defined as:
Y_\Gamma = \mathrm{diag}(|H^⊤1|)
With G=\begin{bmatrix} D^{x}_- B^x \\ D^{y}_- B^y \end{bmatrix} and H=\begin{bmatrix}A^x D^x_- -D^x_- B^x_-\\A^y D^y_- -D^y_- B^y_- \end{bmatrix}
In set_cutcell_matrices! and scalar_transport!, the interface length \chi is thus computed:
χx = (grid.LS[iLS].geoL.dcap[:,:,3] .- grid.LS[iLS].geoL.dcap[:,:,1]) .^ 2
χy = (grid.LS[iLS].geoL.dcap[:,:,4] .- grid.LS[iLS].geoL.dcap[:,:,2]) .^ 2
op.χ[iLS].diag .= sqrt.(vec(χx .+ χy))
System
For the bulk: 1:ni
With the pad function:
pad
A[1:ni,1:ni] = pad(L, -4.0)
A[1:ni,end-nb+1:end] = bc_L_b
The function laplacian returns
L = BxT * iMx * Bx .+ ByT * iMy * By
bc_L[iLS]= BxT * iMx * Hx[iLS] .+ ByT * iMy * Hy[iLS]
bc_L_b = (BxT * iMx_b * Hx_b .+ ByT * iMy_b * Hy_b)
!!! info "Correpondence article-implementation"
cf. [`divergence_B!`](@ref)
ni: bulk so BxT →-G^T , Bx→G in the article
Hx → H
Hx_b →H
So A[1:ni,1:ni] = pad(L, -4.0) corresponds to \mathrm{div}: -G ^ \top W ^ \dagger G
A[1:ni,end-nb+1:end] = bc_L_b
corresponds to -G ^ \top W ^ \dagger H (Hx_b is for the wall)
For the wall: A[end-nb+1:end,:]
A[end-nb+1:end,1:ni] = b_b * (HxT_b * iMx_b' * Bx .+ HyT_b * iMy_b' * By)
corresponds to
I _ b (H ^ \top W ^ \dagger G)
in cf. set_border_matrices!:
cf. bc_matrix_borders! : - sign to the left (cf the outward normal), +sigh to the right
and with
mat_assign_T!(opC_p.HxT_b, sparse(opC_p.Hx_b'))
Hx_b is the transpose of HxT_b
cf. :
A[end-nb+1:end,end-nb+1:end] = pad(b_b * (HxT_b * iMx_bd * Hx_b .+ HyT_b * iMy_bd * Hy_b) .+ χ_b * a1_b, -4.0)
corresponds to
I _ b (H ^ \top W ^ \dagger H) + I _ a I _ \Gamma
!!! todo "-4 or 4"
-4 (convention, value not important ?)
A[1:ni,sb] = bc_L[iLS]
corresponds -G ^ \top W ^ \dagger H (Hx[iLS] is for the wall in laplacian)
For the interface (other than wall)
# Boundary conditions for inner boundaries
A[sb,1:ni] = b * (HxT[iLS] * iMx * Bx .+ HyT[iLS] * iMy * By)
# Contribution to Neumann BC from other boundaries
for i in 1:num.nLS
if i != iLS
A[sb,i*ni+1:(i+1)*ni] = b * (HxT[iLS] * iMx * Hx[i] .+ HyT[iLS] * iMy * Hy[i])
end
end
A[sb,sb] = pad(
b * (HxT[iLS] * iMx * Hx[iLS] .+ HyT[iLS] * iMy * Hy[iLS]) .+ χ[iLS] * a1 .-
a2 * Diagonal(diag(fs_mat)), -4.0
)
#TODO was +b... here in set_poisson
A[sb,end-nb+1:end] = b * (HxT[iLS] * iMx_b * Hx_b .+ HyT[iLS] * iMy_b * Hy_b) #why was it different in set_poisson ?
# Boundary conditions for outer boundaries
A[end-nb+1:end,sb] = b_b * (HxT_b * iMx_b' * Hx[iLS] .+ HyT_b * iMy_b' * Hy[iLS])
end
veci(rhs,grid,iLS+1) .= +χ[iLS] * vec(a0[iLS]) #was - in set_poisson
end
vecb(rhs,grid) .= +χ_b * vec(a0_b) #was - in set_poisson
A[sb,i*ni+1:(i+1)*ni] = b * (HxT[iLS] * iMx * Hx[i] .+ HyT[iLS] * iMy * Hy[i])
!!! todo "Sign free-surface, check"
In solve_poisson:
function solve_poisson(
bc_type, num, grid, a0, opC, opC_u, opC_v,
A, L, bc_L, bc_L_b, BC,
ls_advection)
@unpack Bx, By, Hx, Hy, HxT, HyT, χ, M, iMx, iMy, Hx_b, Hy_b, HxT_b, HyT_b, iMx_b, iMy_b, iMx_bd, iMy_bd, χ_b = opC
ni = grid.nx * grid.ny
nb = 2 * grid.nx + 2 * grid.ny
rhs = fnzeros(grid, num)
a0_b = zeros(nb)
_a1_b = zeros(nb)
_b_b = zeros(nb)
for iLS in 1:num.nLS
set_borders_poisson!(grid, grid.LS[iLS].cl, grid.LS[iLS].u, a0_b, _a1_b, _b_b, BC, num.n_ext_cl)
end
a1_b = Diagonal(vec(_a1_b))
b_b = Diagonal(vec(_b_b))
if ls_advection
# Poisson equation
A[1:ni,1:ni] = pad(L, -4.0)
A[1:ni,end-nb+1:end] = bc_L_b
# Boundary conditions for outer boundaries
A[end-nb+1:end,1:ni] = b_b * ( (HxT_b * iMx_b' * Bx .+ HyT_b * iMy_b' * By) )
A[end-nb+1:end,end-nb+1:end] = pad(b_b * (HxT_b * iMx_bd * Hx_b .+ HyT_b * iMy_bd * Hy_b) .+ χ_b * a1_b,- 4.0)
end
for iLS in 1:num.nLS
if ls_advection
if is_dirichlet(bc_type[iLS])
__a1 = 1.0
__a2 = 0.0
__b = 0.0
elseif is_neumann(bc_type[iLS])
__a1 = 0.0
__a2 = 0.0
__b = 1.0
elseif is_robin(bc_type[iLS])
__a1 = 1.0
__a2 = 0.0
__b = 1.0
elseif is_fs(bc_type[iLS])
__a1 = 0.0
__a2 = 1.0
__b = 0.0
elseif is_wall_no_slip(bc_type[iLS])
__a1 = 0.0
__a2 = 0.0
__b = 1.0
elseif is_navier(bc_type[iLS])
__a1 = 0.0
__a2 = 0.0
__b = 1.0
elseif is_navier_cl(bc_type[iLS])
__a1 = 0.0
__a2 = 0.0
__b = 1.0
else
__a1 = 0.0
__a2 = 0.0
__b = 1.0
end
_a1 = ones(grid) .* __a1
a1 = Diagonal(vec(_a1))
_a2 = ones(grid) .* __a2
a2 = Diagonal(vec(_a2))
_b = ones(grid) .* __b
b = Diagonal(vec(_b))
fs_mat = HxT[iLS] * Hx[iLS] .+ HyT[iLS] * Hy[iLS]
sb = iLS*ni+1:(iLS+1)*ni
# Poisson equation
A[1:ni,sb] = bc_L[iLS]
# Boundary conditions for inner boundaries
A[sb,1:ni] = b * (-1) * (HxT[iLS] * iMx * Bx .+ HyT[iLS] * iMy * By)
# Contribution to Neumann BC from other boundaries
for i in 1:num.nLS
if i != iLS
A[sb,i*ni+1:(i+1)*ni] = b * (-1) * (HxT[iLS] * iMx * Hx[i] .+ HyT[iLS] * iMy * Hy[i])
end
end
A[sb,sb] = pad(
b * (-1) *(HxT[iLS] * iMx * Hx[iLS] .+ HyT[iLS] * iMy * Hy[iLS]) .+ χ[iLS] * a1 .-
a2 * Diagonal(diag(fs_mat)), -4.0
)
#TODO was +b... here in set_poisson
A[sb,end-nb+1:end] = b * (-1) * (HxT[iLS] * iMx_b * Hx_b .+ HyT[iLS] * iMy_b * Hy_b) #why was it different in set_poisson ?
# Boundary conditions for outer boundaries
A[end-nb+1:end,sb] = b_b * (-1) * (HxT_b * iMx_b' * Hx[iLS] .+ HyT_b * iMy_b' * Hy[iLS])
end
veci(rhs,grid,iLS+1) .= +χ[iLS] * vec(a0[iLS]) #was - in set_poisson
end
vecb(rhs,grid) .= +χ_b * vec(a0_b) #was - in set_poisson
return rhs
end
Simultaneous boundary conditions
Corners
!!! todo "TODO"
corners
Poisson equation with variable coefficient
The gradient is not collocated with the scalar, so they don't have the same shapes for each direction. Five diagonal matrices with coeffD inside are required:
mat_coeffDx = Diagonal(vec(coeffDx_bulk)) # coeffDx_bulk is a 2d matrix with shape (grid_u.ny, grid_u.nx), multiplies Bx
mat_coeffDy = Diagonal(vec(coeffDy_bulk)) # coeffDx_bulk is a 2d matrix with shape (grid_v.ny, grid_v.nx), multiplies By
mat_coeffDx_i = Diagonal(vec(coeffDx_interface)) # coeffDx_interface is a 2d matrix with shape (grid_u.ny, grid_u.nx), multiplies Hx
mat_coeffDy_i = Diagonal(vec(coeffDy_interface)) # coeffDx_interface is a 2d matrix with shape (grid_v.ny, grid_v.nx), multiplies Hy
mat_coeffDx_b = Diagonal(vec(coeffD_borders)) # coeffDx_interface is a 1d vector with shape (2grid.ny + 2grid.nx), multiplies Hx_b and Hy_b
\begin{equation}
\grad\left( p^\omega,p^\gamma \right)=W^\dagger \left( G p^\omega + H p^\gamma \right)
\end{equation}
With variable coefficient:
\begin{equation}
\kappa \grad\left( p^\omega,p^\gamma \right)=\kappa^\omega W^\dagger \left( G p^\omega + H p^\gamma \right)
\end{equation}
\begin{equation}
\mathrm{div}\left( \mathbf{q}^\omega,\mathbf{q}^\gamma \right)=-\left( G^T + H^T \right) \mathbf{q}^\omega + H^T \mathbf{q}^\gamma
\end{equation}
\begin{cases}
-\mathrm{div}\left( \mathbf{q}^\omega, \mathbf{q}^\gamma \right) =& Vf^\omega \\
I_a I_\Gamma p^\gamma - I_b \mathrm{div}\left( 0,\kappa^\gamma \mathbf{q}^\gamma \right) =& I_\Gamma g^\gamma \\
\mathbf{q}^\omega =& \kappa \grad\left( p^\omega,p^\gamma \right)\\
\mathbf{q}^\gamma =& \mathbf{q}^\omega
\end{cases}
\begin{cases}
-\mathrm{div}\left( \mathbf{q}^\omega,\mathbf{q}^\gamma \right) =& Vf^\omega \\
I_a I_\Gamma p^\gamma - I_b \mathrm{div}\left( 0,\kappa^\gamma \mathbf{q}^\gamma \right) =& I_\Gamma g^\gamma \\
\mathbf{q}^\omega =& W^\dagger \left( \kappa^\omega G p^\omega + \kappa^\gamma H p^\gamma \right)\\
\mathbf{q}^\gamma =& \mathbf{q}^\omega
\end{cases}
\begin{bmatrix}
G^T W^\dagger \kappa^\omega G & G^T W^\dagger \kappa^\gamma H \\
I_b H^T W^\dagger \kappa^\omega G & I_b H^T W^\dagger \kappa^\gamma H + I_a I_\Gamma
\end{bmatrix}
\begin{bmatrix}
p^\omega \\
p^\gamma
\end{bmatrix} \simeq
\begin{bmatrix}
V f^\omega \\
I_\Gamma g^\gamma
\end{bmatrix}
math
\begin{equation}
\mathrm{div}(UA)=U(divA)+A⋅\grad(U)
\end{equation}
\begin{equation}
\mathrm{div}(\kappa \mathbf{q} )=\kappa (\mathrm{div}(\mathbf{q}))+ \mathbf{q}⋅\grad(\kappa)
\end{equation}
% div(kappa grad(A))=U(divgrad(A))+grad(A)⋅gradU
ni = grid.nx * grid.ny
nb = 2 * grid.nx + 2 * grid.ny
nt = (num.nLS + 1) * ni + nb
Aphi_eleL = spzeros(nt, nt)
phL.phi_ele .= reshape(veci(phL.phi_eleD,grid,1), grid)
A[end-nb+1:end,1:ni] = -b_b * (HxT_b * iMx_b' * coeff[end-nb+1:end,1:ni] * Bx .+ HyT_b * iMy_b' * coeff[end-nb+1:end,1:ni] * By)
@inline zeros(g::Grid) = zeros(g.ny, g.nx)
@inline fnzeros(g::Grid, n::Numerical) = zeros((n.nLS + 1) * g.ny * g.nx + 2 * g.nx + 2 * g.ny)
set_border_matrices!
# Implicit part of heat equation
A[1:ni,1:ni] = pad_crank_nicolson(M .- 0.5 .* timestep_n .* diffusion_coeff_scal .* LT, grid, timestep_n)
A[1:ni,ni+1:2*ni] = - 0.5 .* timestep_n .* diffusion_coeff_scal .* LD
A[1:ni,end-nb+1:end] = - 0.5 .* timestep_n .* diffusion_coeff_scal .* LD_b
# Interior BC
A[ni+1:2*ni,1:ni] = b * (HxT[1] * iMx * Bx .+ HyT[1] * iMy * By)
A[ni+1:2*ni,ni+1:2*ni] = pad(b * (HxT[1] * iMx * Hx[1] .+ HyT[1] * iMy * Hy[1]) .- χ[1] * a1)
A[ni+1:2*ni,end-nb+1:end] = b * (HxT[1] * op.iMx_b * op.Hx_b .+ HyT[1] * op.iMy_b * op.Hy_b)
# Border BCs
A[end-nb+1:end,1:ni] = b_b * (op.HxT_b * op.iMx_b' * Bx .+ op.HyT_b * op.iMy_b' * By)
A[end-nb+1:end,ni+1:2*ni] = b_b * (op.HxT_b * op.iMx_b' * Hx[1] .+ op.HyT_b * op.iMy_b' * op.Hy[1])
A[end-nb+1:end,end-nb+1:end] = pad(b_b * (op.HxT_b * op.iMx_bd * op.Hx_b .+ op.HyT_b * op.iMy_bd * op.Hy_b) .- op.χ_b * a1_b, 4.0)
# Explicit part of heat equation
B[1:ni,1:ni] = M .+ 0.5 .* timestep_n .* diffusion_coeff_scal .* LT .- timestep_n .* CT
B[1:ni,ni+1:2*ni] = 0.5 .* timestep_n .* diffusion_coeff_scal .* LD
B[1:ni,end-nb+1:end] = 0.5 .* timestep_n .* diffusion_coeff_scal .* LD_b
Poisson equation (old implementation)
We consider the Poisson equation -\nabla \cdot \nabla p = f~\text{in}~\Omega and the Robin boundary condition \frac{\partial p}{\partial n} = q = g where a and b are two scalar coefficients
is written:
\color{flblue}{\oint_\Gamma q \cdot n dS} = g \chi
<a name="tagRobin"></a>
Robin boundary condition ap + b \frac{\partial p}{\partial n} = g
The hypothesis q^\omega=q^\gamma=q can be made for system resolution, which simplifies the above expression:
\begin{aligned}
\int_\Omega \nabla \cdot q dV = B_i^x (q_{i+ \frac 12}^x-q_{i - \frac 12}^x) + B_i^y (q_{j+ \frac 12}^y-q_{i - \frac 12}^y)
\end{aligned}
"The scalar value along the immersed boundary ( p^\gamma ) is not explicitly known, and has to be inferred from the values of the normal component of the gradient along the immersed boundary and g. The discretization of this equation requires the definition of the boundary vectors g^\gamma , a^\gamma and b^\gamma for boundary cells, whose components are set to,
g^\gamma = g(x_{i,j}^\gamma, y_{i,j}^\gamma)
a^\gamma = a(x_{i,j}^\gamma, y_{i,j}^\gamma)
b^\gamma = b(x_{i,j}^\gamma, y_{i,j}^\gamma)
where a^\gamma and b^\gamma define the diagonal matrices I_a and I_b, I_a = \mathrm{diag} (a^\gamma) and I_b = \mathrm{diag} (b^\gamma). Likewise, the components of the discrete load f^\omega are defined as f^\omega_{i,j} = f ( x_{i,j}^\gamma, y_{i,j}^\gamma ) in fluid cells and zero elsewhere.
"
See this equation for the
<a href="documentation.html#divergencetag">divergence</a> and:
" In order to discretize the Poisson and Robin
See this equation for the
<a href="test.html#tagPoisson">Poisson</a> and <a href="test.html#tagRobin">Robin</a> equations
, the term with the derivative in the normal direction (\partial_n p) is set as the boundary contribution to the divergence. The resulting system is given by:
\begin{cases}
− \mathrm{div} (q^\omega, q^\gamma ) &= V f^\omega\\
I_a I p^\gamma − I_b \mathrm{div} (0, q^\omega) &= I g^\omega\\
q^\omega &= \mathrm{grad} ( p^\omega, p^\gamma ) \\
q^\gamma &= q^\omega
\end{cases}
where I measures the boundary within a given cell.
"
"
The intermediate variables q^\omega and q^\gamma can be eliminated from the above linear system, which can then be expressed in terms of the previously defined matrices as
\begin{equation}
\left [ \begin{array}{>{\centering\arraybackslash$} p{2.0cm} <{$} >{\centering\arraybackslash$} p{3.2cm} <{$}}
G ^ \top W ^ \dagger G & G ^ \top W ^ \dagger H \\
I _ b H ^ \top W ^ \dagger G & I _ b H ^ \top W ^ \dagger H + I _ a I _ \Gamma
\end{array} \right ] \left [ \begin{array}{c}
p ^ \omega \\
p ^ \gamma
\end{array} \right ] \simeq \left [ \begin{array}{c}
V f ^ \omega \\
I _ \Gamma g ^ \gamma
\end{array} \right ],
\label{eq:roblapmat}
\end{equation}
"
by eliminating all intermediate variables and using Eqs.
<a href="documentation.html#grad_concise">(gradient)</a> and <a href="documentation.html#div_concise">(divergence)</a>.
"The matrix on the left-hand side is symmetric if the coefficients b are equal to one and it is straightforward to obtain a symmetric matrix if they are not. Appendix C presents the resulting expression for a given cell of the one-dimensional Poisson's equation, for the sake of clarity."
From (Rodriguez et al. 2024)
"
$I^\Gamma$ is needed in order to have a dimensionally consistent equation for the
See this equation for the
<a href="test.html#tagRobin">Robin</a> equations.
The intermediate variables q^\omega and q^\gamma can be eliminated from the above linear system, which can then be expressed in terms of the previously defined matrices as
\begin{bmatrix}
G W^\dagger G & G W^\dagger H\\
I_b H W ^\dagger G & I_b H W ^\dagger H + I_a I_\Gamma
\end{bmatrix}
\begin{bmatrix}
p^\omega\\
p^\gamma
\end{bmatrix} =
\begin{bmatrix}
V f^\omega\\
I g^\gamma
\end{bmatrix}
"
At the interface
Dirichlet __a1 = -1.0
Neumann __b = 1.0
Robin __a1 = -1.0
__a2 = 0.0 __b = 1.0
In we have
In set_borders!, we have:
Dirichlet: a1 = -1 Neumann: b =1 Robin: a1 = -1 b = 1 a0 : BC value
!!! todo "Signs"
* why a1 = -1
* why -a0
A[sb,sb] = -pad(
b * (HxT[iLS] * iMx * Hx[iLS] .+ HyT[iLS] * iMy * Hy[iLS]) .- χ[iLS] * a1 .+
a2 * Diagonal(diag(fs_mat)), 4.0
)
A[sb,sb] = pad(
b * (-1) * (HxT[iLS] * iMx * Hx[iLS] .+ HyT[iLS] * iMy * Hy[iLS]) .+ χ[iLS] * a1 .-
a2 * Diagonal(diag(fs_mat)), -4.0
)
!!! todo
``a1 = -1``
A[sb,end-nb+1:end] = b * (HxT[iLS] * iMx_b * Hx_b .+ HyT[iLS] * iMy_b * Hy_b)
Why +b? and everywhere else -b ?
A[sb,end-nb+1:end] = b * (HxT[iLS] * iMx_b * Hx_b .+ HyT[iLS] * iMy_b * Hy_b)
# Boundary conditions for outer boundaries
A[end-nb+1:end,sb] = -b_b * (HxT_b * iMx_b' * Hx[iLS] .+ HyT_b * iMy_b' * Hy[iLS])
end
veci(rhs,grid,iLS+1) .= -χ[iLS] * vec(a0[iLS])
end
vecb(rhs,grid) .= -χ_b * vec(a0_b)
set_poisson
cf. (Rodriguez et al. 2024) for a 1D expression in the i-th cell, x component
\begin{aligned}
-\mathcal{B}_{x,i} [&\mathcal{W}_{x,i+1} (\mathcal{B}_{x,i+1} p^\omega_{i+1} - \mathcal{B}_{x,i} p^\omega _i ) - \mathcal{W}_{x,i} (\mathcal{B}_{x,i} p^\omega_i - \mathcal{B}_{x,i-1} p^\omega_{i-1} )] \\
-\mathcal{B}_{x,i} \{&\mathcal{W}_{x,i+1} [(\mathcal{B}_{x,i+1} - \mathcal{A}_{x,i+1} )p^\gamma_{i+1} - (\mathcal{A}_{x,i+1} - \mathcal{B}_{x,i} ) p^\gamma_i ] \\
&-\mathcal{W}_{x,i} [(\mathcal{B}_{x,i} - A_{x,i} ) p^\gamma_i + (A_{x,i} - \mathcal{B}_{x,i-1}) p^\gamma_{i-1} ] \} \\
= V_i f^\omega_i&
\end{aligned}
function set_poisson(
bc_type, num, grid, a0, opC, opC_u, opC_v,
A, L, bc_L, bc_L_b, BC,
ls_advection)
@unpack Bx, By, Hx, Hy, HxT, HyT, χ, M, iMx, iMy, Hx_b, Hy_b, HxT_b, HyT_b, iMx_b, iMy_b, iMx_bd, iMy_bd, χ_b = opC
ni = grid.nx * grid.ny
nb = 2 * grid.nx + 2 * grid.ny
rhs = fnzeros(grid, num)
a0_b = zeros(nb)
_a1_b = zeros(nb)
_b_b = zeros(nb)
# Dirichlet: a1 = -1
# Neumann: b =1
# Robin: a1 = -1 b = 1
# a0 : BC value
for iLS in 1:num.nLS
set_borders!(grid, grid.LS[iLS].cl, grid.LS[iLS].u, a0_b, _a1_b, _b_b, BC, num.n_ext_cl)
end
a1_b = Diagonal(vec(_a1_b))
b_b = Diagonal(vec(_b_b))
if ls_advection
# Poisson equation
A[1:ni,1:ni] = pad(L, -4.0)
A[1:ni,end-nb+1:end] = bc_L_b
# Boundary conditions for outer boundaries
A[end-nb+1:end,1:ni] = -b_b * (HxT_b * iMx_b' * Bx .+ HyT_b * iMy_b' * By)
A[end-nb+1:end,end-nb+1:end] = -pad(b_b * (HxT_b * iMx_bd * Hx_b .+ HyT_b * iMy_bd * Hy_b) .- χ_b * a1_b, 4.0)
end
for iLS in 1:num.nLS
if ls_advection
if is_dirichlet(bc_type[iLS])
__a1 = -1.0
__a2 = 0.0
__b = 0.0
elseif is_neumann(bc_type[iLS])
__a1 = 0.0
__a2 = 0.0
__b = 1.0
elseif is_robin(bc_type[iLS])
__a1 = -1.0
__a2 = 0.0
__b = 1.0
elseif is_fs(bc_type[iLS])
__a1 = 0.0
__a2 = 1.0
__b = 0.0
elseif is_wall_no_slip(bc_type[iLS])
__a1 = 0.0
__a2 = 0.0
__b = 1.0
elseif is_navier(bc_type[iLS])
__a1 = 0.0
__a2 = 0.0
__b = 1.0
elseif is_navier_cl(bc_type[iLS])
__a1 = 0.0
__a2 = 0.0
__b = 1.0
else
__a1 = 0.0
__a2 = 0.0
__b = 1.0
end
_a1 = ones(grid) .* __a1
a1 = Diagonal(vec(_a1))
_a2 = ones(grid) .* __a2
a2 = Diagonal(vec(_a2))
_b = ones(grid) .* __b
b = Diagonal(vec(_b))
fs_mat = HxT[iLS] * Hx[iLS] .+ HyT[iLS] * Hy[iLS]
sb = iLS*ni+1:(iLS+1)*ni
# Poisson equation
A[1:ni,sb] = bc_L[iLS]
# Boundary conditions for inner boundaries
A[sb,1:ni] = -b * (HxT[iLS] * iMx * Bx .+ HyT[iLS] * iMy * By)
# Contribution to Neumann BC from other boundaries
for i in 1:num.nLS
if i != iLS
A[sb,i*ni+1:(i+1)*ni] = -b * (HxT[iLS] * iMx * Hx[i] .+ HyT[iLS] * iMy * Hy[i])
end
end
A[sb,sb] = -pad(
b * (HxT[iLS] * iMx * Hx[iLS] .+ HyT[iLS] * iMy * Hy[iLS]) .- χ[iLS] * a1 .+
a2 * Diagonal(diag(fs_mat)), 4.0
)
A[sb,end-nb+1:end] = b * (HxT[iLS] * iMx_b * Hx_b .+ HyT[iLS] * iMy_b * Hy_b)
# Boundary conditions for outer boundaries
A[end-nb+1:end,sb] = -b_b * (HxT_b * iMx_b' * Hx[iLS] .+ HyT_b * iMy_b' * Hy[iLS])
end
veci(rhs,grid,iLS+1) .= -χ[iLS] * vec(a0[iLS])
end
vecb(rhs,grid) .= -χ_b * vec(a0_b)
return rhs
end
Boundary conditions
Neumann boundary condition
The interface length is used to impose boundary conditions, for instance, a Neumann boundary condition \frac{\partial p}{\partial n} = q = g is written:
\color{flblue}{\oint_\Gamma q \cdot n dS} = g \chi
Robin boundary condition
ap + b \frac{\partial p}{\partial n} = g
Assumption
The hypothesis q^\omega=q^\gamma=q can be made for system resolution, which simplifies the above expression:
\begin{aligned}
\int_\Omega \nabla \cdot q dV = B_i^x (q_{i+ \frac 12}^x-q_{i - \frac 12}^x) + B_i^y (q_{j+ \frac 12}^y-q_{i - \frac 12}^y)
\end{aligned}
Other formalisms
The aforementioned operators may be rewritten in a more concise form.
<a name="grad_concise"></a>
\mathrm{grad}(q^ω, q^{γ} ) = W ^ \dagger \left ( G p ^ \omega + H p ^ \gamma \right )
\mathrm{div}(q^ω, q^{γ} ) = − G^⊤ + H^{Γ1,⊤} + H^{Γ2,⊤} q^ω + H^{Γ1,⊤}q^{γ1} + H^{Γ2,⊤}q^{γ2}
where
−G^⊤ − H^{Γ1,⊤} − H^{Γ2,⊤} = D_x+A_x^\Omega D_y+A_y^\Omega
and
H^{Γi,⊤} = B_x^{Γi} D_x+ − D_x+A_x^{Γi} B_y^{Γi} D_y+ − D_y+A_y^{Γi}
" It should be noted that if the vector field represents a gradient, which is a higher order quantity, no distinction is made between its values in the fluid and at the boundary, i.e. $q^\gamma$ is simply set to $q^ \omega$, simplifying the divergence operator to
<a name="div_concise"></a>
\mathrm{div}(q^ω, q^{γ} ) = -G ^T q^ \omega
cf fig 2.3 in Rodriguez 2024
Divergence of a face-centered vector field
<figure>
<a name="geometric_moments"></a>
<img src="./assets/moments_doc.svg" alt="Geometric moments" title="Geometric moments">
<figcaption>"Geometric moments. $\mathcal{A} _ x$ in red dashed lines, $\mathcal{A} _ y$ in blue dashed lines, $\mathcal{B} _ x$ in green dashed lines. $\mathcal{B} _ y$ are not shown."</figcaption>
</figure>
!!! todo "Staggered capacities and gradients"
for Poiseuille
- u grid: xmin, xmin+dx, ...xmax-dx,xmax
v grid: ymin, ymin+dx, ...ymax-dy,ymax
set_bc_bndsif is_neumann(BC_v.bottom) @inbounds Dy[1,:] .= v[1,:] .+ Hv[1,:] .* BC_v.bottom.val elseif is_dirichlet(BC_v.bottom) @inbounds Dy[1,:] .= BC_v.bottom.val elseif is_periodic(BC_v.bottom) @inbounds Dy[1,:] .= v[end,:] end
but gamma for u v grids ?
!!! todo "add param special init BC"
" As explained in Chapter 1, a staggered grid arrangement has to be employed to discretize the incompressible Navier-Stokes equations if a conservative method is to be developed. The use of staggered quantities means that Eqs. (2.2) have to be applied to staggered finite volumes. This leads to the need of additional geometrical information as follows. For each domain $\Omegai$, two staggered grids are defined in the $x$ and $y$-directions. These staggered grids are defined by tracing vertical and horizontal lines passing through each corresponding cell centroid $(x{i,j}^c, y_{i,j}^c)$ (and through the cell centers if the domain $\Omega_i$ is not defined in the cell). The intersection of these lines with the domain $\Omega_i$ define two additional surface moments $B_x$ and $B_y$. The planar faces $B_x$, $A_y$ and the boundary surface $\Gamma$ enclose the staggered volumes $W_x$ and likewise, the planar faces $B_y$, $A_x$ and the boundary surface $\Gamma$ define the staggered volumes $Wy$. For the sake of presentation, in Fig. 2.3, the $x$-staggered grid and the staggered volume $W{x,(i+1,j-1)}$ are also shown. These additional geometric moments define an additional set of diagonal matrices.
"
" To avoid any confusions, the $\alpha$ component of a gradient will be denoted $\left[\text{grad}(p^{\omega}, p^{\eta})\right]{\alpha}$ for a scalar field $p$ and $\left[\text{grad}{\beta} (u{\beta}^{\omega}, u{\beta}^{\eta})\right]_{\alpha}$ for the $\beta$ component of a vector field $u$.
"
cf A levelset based cut-cell method for two-phase flows. Part 2: Free-surface flows and dynamic contact angle treatment:
"Just like the gradient operator, the (volume-integrated) divergence operator is multilinear and takes one more argument than the number of interfaces:
\mathrm{div}(q^ω, q^{γ1} , q^{γ2} ) = − G^⊤ + H^{Γ1,⊤} + H^{Γ2,⊤} q^ω + H^{Γ1,⊤}q^{γ1} + H^{Γ2,⊤}q^{γ2}
where
−G^⊤ − H^{Γ1,⊤} − H^{Γ2,⊤} = D_x+A_x^\Omega D_y+A_y^\Omega
and
H^{Γi,⊤} = B_x^{Γi} D_x+ − D_x+A_x^{Γi} B_y^{Γi} D_y+ − D_y+A_y^{Γi}
The reader may find details in Rodriguez 2024 and Fullana (2022)
for Morinishi's notation for proofs.
Capacities
See section 2.2.1 Geometric moments in "Numerical methods for modeling and optimization of interfacial flows"
opC_p.Hx_b: left: -dy (cell height along y), bottom: 0, right: +dy, top: 0
opC_p.Hy_b: left: 0 , bottom: -dx, right: 0, top: +dx
Geometry
Example:
x = LinRange(0, L0, n+1)
y = LinRange(0, L0, n+1)
Or
x = collect(LinRange(0, L0, n+1))
y = collect(LinRange(0, L0, n+1))
That sets the location of the cell faces, not the cell center. The borders of the domain will be at 0 and at L0.
Along x-axis, without interface, constant mesh spacing
The first cell centroid (when there is no two-fluid interface and the mesh spacing is constant) is at:
x[1]+dx[1]/2
i.e.
xmin+dx[1]/2
The cell centroids are: x[1]+dx[1]/2, x[1]+dx[1]*3/2,..., x[end]-dx[1]3/2,x[end]-dx[1]/2
The meshes are initialized by init_meshes which in turn calls Mesh
init_meshes
Mesh
The borders of the domain are defined by x and y, now stored ax x\_node and y\_node in the Mesh structure, you can also compute the coordinates of the borders of the domain with:
x_bc_left = gp.x[:,1] .- gp.dx[:,1] ./ 2.0
y_bc_bottom = gp.y[1,:] .- gp.dy[1,:] ./ 2.0
y_bc_top = gp.y[end,:] .+ gp.dy[end,:] ./ 2.0
x_bc_right = gp.x[:,end] .+ gp.dx[:,end] ./ 2.0
Storage of bulk and interfacial variables
Variables are stored in the following order: bulk, interfacial (based on levelsets) and interfacial (borders, in the order left bottom right top, cf grid.ind.b_left[1], grid.ind.b_bottom[1], grid.ind.b_right[1] and grid.ind.b_top[1]).
The system size is nt \times nt:
ni = gp.nx * gp.ny #size of bulk/interfacial (LS) data
nb = 2 * gp.nx + 2 * gp.ny #size of border
nt = (num.nLS + 1) * ni + nb
A = spzeros(nt, nt)
Example at (j=1,i=1), index of corner values
i_corner_vecb_left = (num.nLS + 1) * ni + 1
i_corner_vecb_bottom = (num.nLS + 1) * ni + gp.ny + 1
You can use the following functions to access bulk and interfacial fields:
veci
vecb(a, g::G) where {G<:Grid} = @view a[end-2*g.ny-2*g.nx+1:end]
vecb
vecb_L(a,g::G) where {G<:Grid} = @view vecb(a, g)[1:g.ny]
vecb_L
vecb_B(a,g::G) where {G<:Grid} = @view vecb(a, g)[g.ny+1:g.ny+g.nx]
vecb_B
vecb_R(a,g::G) where {G<:Grid} = @view vecb(a, g)[g.ny+g.nx+1:2*g.ny+g.nx]
vecb_R
vecb_T(a,g::G) where {G<:Grid} = @view vecb(a, g)[2*g.ny+g.nx+1:2*g.ny+2*g.nx]
vecb_T
Indices
In Flower.jl, the dimensions are stored in this order (y,x) i.e. to access the data of cell (i,j), one may use data[j,i].
Indices
Mesh and how to access data at (i,j)
For a 1D vector, grid p
II = CartesianIndex(jplot, iplot) #(id_y, id_x)
pII = lexicographic(II, grid.ny)
!!! todo "TODO CHECK"
For a 1D vector, grid v
```julia
II = CartesianIndex(jplot, iplot) #(id_y, id_x)
pII = lexicographic(II, grid.ny +1)
```
diag[pII] is faster than [pII,pII]
Accessing iMx (inverse of weight)
The function lexicographic can be used to access iMx (inverse of weight) in a cell.
II = CartesianIndex(j,i)
tmp = lexicographic(II,gp.ny)
II is a two-dimensional index, and in iMx 1D index is required.
It is recommended to use iMx.diag[id], it's faster than accessing iMx[id,id].
!!! todo "TODO
Do I need to add 1 to n for the border terms ? ``pII = lexicographic(δy⁺(II), grid.ny)``
In set_cutcell_matrices!, grid.ny already has the good size no matter what grid you're using
But for all the grids, when you set iMy, you have to do
pII = lexicographic(δy⁺(II), grid.ny+1) #for the top border
pII = lexicographic(II, grid.ny+1) #for the rest
Because iMy has one extra layer on y (true for all the grids).
!!! todo "TODO"
v mesh left wall : in the x-direction the distances are h and h/2
Initialisation
When a Neumann Boundary Condition (BC) is imposed at the interface, an extrapolation can be used to get an approximation of the value at the interface.
get_height!
init_fields_multiple_levelsets!
get height for u grid: solid then liquid
H,LS.mid_point, dx, dy and geo.centroid are of size (ny, nx+1)
get height for v grid: solid then liquid
H,LS.mid_point, dx, dy and geo.centroid are of size (ny+1, nx)
get height for p grid: solid then liquid
H, LS.mid_point, dx, dy and geo.centroid are of size (ny, nx)
See Definitions
!!! todo "TODO"
finish doc on centroids
```julia
x_centroid = gp.x .+ getproperty.(gp.LS[1].geoL.centroid, :x) .* gp.dx
y_centroid = gp.y .+ getproperty.(gp.LS[1].geoL.centroid, :y) .* gp.dy
x_bc = gp.x .+ getproperty.(gp.LS[1].mid_point, :x) .* gp.dx
y_bc = gp.y .+ getproperty.(gp.LS[1].mid_point, :y) .* gp.dy
#Initialize bulk value
vec1(phL.TD,gp) .= vec(ftest.(x_centroid,y_centroid))
#Initialize interfacial value
vec2(phL.TD,gp) .= vec(ftest.(x_bc,y_bc))
#Initialize border value
vecb_L(phL.TD, gp) .= ftest.(gp.x[:,1] .- gp.x[:,1] ./ 2.0, gp.y[:,1])
vecb_B(phL.TD, gp) .= ftest.(gp.x[1,:], gp.y[1,:] .- gp.y[1,:] ./ 2.0)
vecb_R(phL.TD, gp) .= ftest.(gp.x[:,end] .+ gp.x[:,1] ./ 2.0, gp.y[:,1])
vecb_T(phL.TD, gp) .= ftest.(gp.x[1,:], gp.y[end,:] .+ gp.y[1,:] ./ 2.0)
```
Electrical potential
Poisson equation for electrical potential
Variable coefficients are interpolated. The Laplacian matrices L, bc_L_b and bc_L are defined in the functions laplacian and laplacian_bc. coeffD is a coefficient involving the electrical conductivity. Then, you have to put the coefficient coeffD before the matrices Bx, By, Hx, Hy, Hx_b and Hy_b. The thing is that the gradient is not collocated with the scalar, so they don't have the same shapes for each direction. Five diagonal matrices with coeffD inside are required.
In set_borders!, we have:
- Dirichlet: a1 = -1
- Neumann: b =1
- Robin: a1 = -1 b = 1
a0 : BC value
interpolate_scalar!(grid, grid_u, grid_v, u, uu, uv)aux_interpolate_scalar!(II_0, II, u, x, y, dx, dy, u_faces) static_stencil(a, II::CartesianIndex) faces_scalar(itp, II_0, II, x, y, dx[II], dy[II])\phi(\bar{x}, \bar{y}) = \begin{bmatrix} \bar{x}^2 & \bar{x} & 1 \end{bmatrix} \begin{bmatrix} m_{0,0} & m_{0,1} & m_{0,2} \\ m_{1,0} & m_{1,1} & m_{1,2} \\ m_{2,0} & m_{2,1} & m_{2,2} \end{bmatrix} \begin{bmatrix} \bar{y}^2 \\ \bar{y} \\ 1 \end{bmatrix} = X^T MY
This matrix is computed in B_BT.
B_BT(II_0, x, y)
where $\bar{x} = \frac{x - x{i,j}}{x{i+1,j} - x{i-1,j}}$ and $\bar{y} = \frac{y - y{i,j}}{y{i,j+1} - y{i,j-1}}$ are normalized Cartesian coordinates, and $M$ is a matrix with the unknown interpolation coefficients. The values of the level-set at the cells of the $3 \times 3$ stencil can be used to set the following system of equations
\Phi = AMB.
where
\Phi = \begin{bmatrix}
\phi_{i-1,j-1} & \phi_{i-1,j} & \phi_{i-1,j+1} \\
\phi_{i,j-1} & \phi_{i,j} & \phi_{i,j+1} \\
\phi_{i+1,j-1} & \phi_{i+1,j} & \phi_{i+1,j+1}
\end{bmatrix},
A = \begin{bmatrix}
(\bar{x}_{i-1,j})^2 & \bar{x}_{i-1,j} & 1 \\
0 & 0 & 1 \\
(\bar{x}_{i+1,j})^2 & \bar{x}_{i+1,j} & 1
\end{bmatrix},
B = \begin{bmatrix}
(\bar{y}_{i,j-1})^2 & 0 & (\bar{y}_{i,j+1})^2 \\
\bar{y}_{i,j-1} & 0 & \bar{y}_{i,j+1} \\
1 & 1 & 1
\end{bmatrix}.
The interpolation matrix $M$ can be computed after inverting matrices $A$ and $B$
M = A^{-1}\Phi B^{-1}.
Once the interpolation matrix is obtained, the following equation can be used to interpolate the value of the level-set at a normalized position $(\bar{x}, \bar{y})$:
\phi(\bar{x}, \bar{y}) = \begin{bmatrix} \bar{x}^2 & \bar{x} & 1 \end{bmatrix}
\begin{bmatrix}
m_{0,0} & m_{0,1} & m_{0,2} \\
m_{1,0} & m_{1,1} & m_{1,2} \\
m_{2,0} & m_{2,1} & m_{2,2}
\end{bmatrix}
\begin{bmatrix} \bar{y}^2 \\ \bar{y} \\ 1 \end{bmatrix} = X^T MY
"
The level-set has to be interpolated to compute the required geometric moments.
To that end, a biquadratic interpolation defined on the 3 \times 3 stencil centered at a cell (i, j) will be used.
The interpolant is defined as
"
biquadratic(m, x, y)
Example
To determine the values of the corners, we perform a bi-quadratic interpolation on the 3 x 3 stencil centered on the cell of interest. For this exercise, we assume a constant spacing of the Cartesian grid in all dimensions. We want to calculate the
The interpolation matrix ( M ) is given by:
[ \phi(x, y) = \sum{i=0}^{2} \sum{j=0}^{2} m_{ij} x^i y^j ]
This can be expressed as:
[ \phi(x, y) = \begin{bmatrix} x^2 & x & 1 \end{bmatrix} \begin{bmatrix} m{2,2} & m{2,1} & m{2,0} \ m{1,2} & m{1,1} & m{1,0} \ m{0,2} & m{0,1} & m_{0,0} \end{bmatrix} \begin{bmatrix} y^2 \ y \ 1 \end{bmatrix} = X M Y^T ]
With (\phi(x, y)) being the set of known data points of the level set function, yielding:
[ \Phi = G M G^T ]
Where:
[ \Phi = \begin{bmatrix} \phi{-1,-1} & \phi{-1,0} & \phi{-1,1} \ \phi{0,-1} & \phi{0,0} & \phi{0,1} \ \phi{1,-1} & \phi{1,0} & \phi{1,1} \end{bmatrix} = \begin{bmatrix} (-1)^2 & -1 & 1 \ 0^2 & 0 & 1 \ 1^2 & 1 & 1 \end{bmatrix} G \begin{bmatrix} m{2,2} & m{2,1} & m{2,0} \ m{1,2} & m{1,1} & m{1,0} \ m{0,2} & m{0,1} & m{0,0} \end{bmatrix} M \begin{bmatrix} (-1)^2 & 0 & (1)^2 \ -1 & 0 & 1 \ 1 & 1 & 1 \end{bmatrix} G^T ]
!!! todo "Interpolation for Poisson with variable coefficient"
does not take border and interfacial values into account
For the wall term on the left part of the domain, corresponding face control volume <a href="#facecontrolvolume"> (see figure), we should use: $\kappa = \frac{\kappa^\gamma+\kappa^\omega_{1,j}}{2}$.<span><img src="./assets/Wx_mesh_square_bubble_wall_2.svg" alt="image" height="800" /></span></a></span>
interpolate_scalar_to_staggered_u_v_grids_at_border!
solve_poisson_variable_coeff!
Laplacian: bulk
L = BxT * iMx * Bx
L is of size nx*ny tmp_x ((nx+1)*ny , nx*ny) BxT (nx*ny, (nx+1)*ny ) Bx ((nx+1)*ny, nx*ny) iMx ((nx+1)*ny, (nx+1)*ny ) iMx_b (nx+1)*ny,2*nx +2*ny Hx_b 2*nx +2*ny,2*nx +2*ny
The main loop for the resolution of the Poisson equation is in:
solve_poisson_loop!
Electrical current
compute_grad_phi_ele!
Scalar transport
scalar_transport!
The concentrations are checked during the scalar transport:
compute_bulk_or_interface_average
Phase change
integrate_mass_transfer_rate_over_interface
Electrical conductivity
update_electrical_conductivity!
Fresh and dead cells
" The last step of the method presented in this study is the treatment of fresh and dead cells. As the interface moves through the Cartesian grid, the temperature field needs to be initialized according to the front position. Two cases can be distinguished:
- A full or partial cell becomes an empty cell.
- An empty cell becomes a partial cell.
In the first case, the previous temperature value is simply set to 0. In the second case, the temperature previously non-existing needs to be initialized. The algorithm is similar to that of the Stefan condition (Section 4.2.1)
A shifted $3 \times 3$ stencil is chosen, as shown in Figure 4.7.
A line from the interface centroid is cast in the opposite normal direction $-\mathbf{n}$.
The crossing points $A$ and $B$ of this line and the vertical (or horizontal depending on the normal orientation) segments of the neighboring 3 points, are identified. "
Fullana (2022)kill_dead_cells!
Butler-Volmer equation
update_BC_electrical_potential_left!(num,grid,BC_phi_ele,elec_cond,elec_condD,i_butler)
update_electrical_current_from_Butler_Volmer!(num,grid,heat,phi_eleD,i_butler;T=nothing)
butler_volmer_concentration(alpha_a,alpha_c,c_H2,c0_H2,c_H2O,c0_H2O,c_KOH,c0_KOH,Faraday,i0,phi_ele,phi_ele1,Ru,temperature0)
butler_volmer_no_concentration(alpha_a,alpha_c,Faraday,i0,phi_ele,phi_ele1,Ru,temperature0)
butler_volmer_no_concentration!(alpha_a,alpha_c,Faraday,i0,phi_ele,phi_ele1,Ru,temperature0,i_current)
butler_volmer_no_concentration_concentration_Neumann(alpha_a,alpha_c,Faraday,i0,phi_ele,phi_ele1,Ru,temperature0,diffusion_coeff,inv_stoechiometric_coeff)
Heat equation
set_heat!
Navier-Stokes
!!! info "Operators"
"I do this:
```julia
opC_u.AxT * opC_u.Rx
```
so that the same staggered volume is used for all the terms. Otherwise, there were some velocity cells that were not seeing any pressure gradient"
FE_set_momentum
Navier-slip boundary conditions
" If the Navier-slip boundary condition is to be applied to a straight surface, the Navier-slip model
\begin{cases}
u_t &= \lambda \frac{\partial u_t}{\partial n} \\
u_n &= 0 \\
\end{cases}
can be discretized with Robin boundary conditions: ap + b \frac{\partial p}{\partial n} = g
" For a general Robin boundary condition as in Eq. (2.44) applied on the boundary $\Gamma_i$, the discretized version reads
\begin{equation} \label{eq:robin_discretized}
Y^{\Gamma_i}_a Y^{\Gamma_i} p^{\gamma_i} + Y^{\Gamma_i}_b H^{\Gamma_i, \top} W^{\Gamma_i} \left[ G p^{\omega} + H^{\Gamma_i} p^{\gamma_i} + H^{\Gamma_{3-i}} p^{\gamma_{3-i}} \right] = Y^{\Gamma_i} g^{\gamma_i},
\end{equation}
where the diagonal matrices $Y^{\Gamma_i}_a = \text{diag}(a^{\gamma_i})$ and $Y^{\Gamma_i}_b = \text{diag}(b^{\gamma_i})$ are coefficients that can be used to impose either Dirichlet, Neumann or Robin boundary conditions, the diagonal matrix $Y^{\Gamma_i} = \text{diag}(|H^{\Gamma_i, \top} \mathbf{1}|)$ measures the boundary within a given cell.
It should be noted that the last term on the left-hand-side of Eq. (\ref{eq:robindiscretized}) represents the contribution from the second boundary: $\Gamma{3-i}$
\Gamma_{3-i}
to the projection of the gradient onto the direction normal to $\Gamma_i$.
"
. However, if the boundary condition is applied to an arbitrarily oriented boundary surface, some modifications are required.
For an arbitrary orientation, the tangential component of the velocity $u_t$ depends on both velocity components $u_x$ and $u_y$. As presented in Sec. \ref{sec:2.4}, the proposed methodology (and operators previously defined) leverages on bulk and boundary values, the later being defined implicitly as roots of the discretized boundary conditions. This results in a system of algebraic equation that consists of the discretized bulk (momentum) and boundary conditions, that is solved to determine both bulk and boundary fields.
In Part I, both boundary velocity components ($u_x^i$ and $u_y^i$) were solved. However, for the Navier-slip condition, a different approach is used: in the absence of blowing or suction, the velocity component normal to the boundary is identically zero. As far as the velocity component tangential to the boundary is concerned, it is implicitly defined by a discretized form of Eq. \eqref{eq:3.7} (see below). For simplicity, the boundary field $u_t^i$ is collocated with the cells where $p$ is defined. However, in order to assemble the viscous and convective transport terms and the velocity boundary conditions, this information (together with the vanishing normal velocity) must be interpolated to the mesh faces before being included in the viscous transport term.
To illustrate this, let us consider the case where the Navier-slip condition (Eq. \eqref{eq:3.7}) is applied on $\Gammai$, and a different boundary condition (e.g., a free-surface or Dirichlet boundary condition) on $\Gamma{3-i}$. The discrete version of Eq. \eqref{eq:3.7} then reads
\begin{equation}
Y^{\Gamma_i} u_t^{\gamma_i} - \lambda \left[ H^{\Gamma_i, \top} W^{\Gamma} H^{\Gamma_i} u_t^{\gamma_i} +
H^{\Gamma_i, \top} W^{\Gamma} H^{\Gamma_{3-i}} \left( \text{diag}(\sin \alpha) S_x^- u_x^{\gamma_{3-i}} + \text{diag}(-\cos \alpha) S_y^- u_y^{\gamma_{3-i}} \right) +
H^{\Gamma_i, \top} W^{\Gamma} G \left( \text{diag}(\sin \alpha) S_x^- u_x^{\gamma_i} + \text{diag}(-\cos \alpha) S_y^- u_y^{\gamma_i} \right) \right] = Y^{\Gamma_i} U_b. \tag{3.60}
\end{equation}
The projection and interpolation of the tangential velocity onto the mesh faces can be performed following
\begin{align}
u_x^{\gamma_i} &= \text{diag}(\sin \alpha) S_x^+ u_t^{\gamma_i}, \tag{3.61} \\
u_y^{\gamma_i} &= \text{diag}(-\cos \alpha) S_y^+ u_t^{\gamma_i},
\end{align}
which is used in the viscous and convective transport terms and in boundary conditions where the Cartesian components of the velocity appear, as mentioned earlier.
Regarding the velocity divergence which appears in the right-hand-side of the Poisson's equation for the pressure, there is no contribution from the Navier-slip boundary since there is no velocity in the normal direction.
"
set_velocity_boundary_conditions
The resolution is coupled if a Navier BC is activated.
FE_set_momentum_coupled
!!! todo "change bc_type to BCu and BCv"
New version
!!! todo "capacities for divergence"
```julia
Duv = opC_p.AxT * vec1(u_predictionD,grid_u) .+ opC_p.Gx_b * vecb(u_predictionD,grid_u) .+
opC_p.AyT * vec1(v_predictionD,grid_v) .+ opC_p.Gy_b * vecb(v_predictionD,grid_v)
for iLS in 1:nLS
if !is_navier(bc_int[iLS]) && !is_navier_cl(bc_int[iLS])
Duv .+= opC_p.Gx[iLS] * veci(u_predictionD,grid_u,iLS+1) .+
opC_p.Gy[iLS] * veci(v_predictionD,grid_v,iLS+1)
end
end
```
\begin{equation}
\begin{pmatrix}
F_u & B_p^T \\
B_u & 0
\end{pmatrix}
\begin{pmatrix}
\vec{u}^{n+1} \\
p^{n+1}
\end{pmatrix}
=
\begin{pmatrix}
\vec{f} \\
0
\end{pmatrix}
\end{equation}
where ( F_u = M_u^{(\rho)} + N_u^{(\rho)} + L_u^{(\mu)} ), and finally:
- ( Bu \vec{u}^{n+1} \approx \int{\Omega_{i,j,k}} \nabla \cdot \vec{u}^{n+1} dV ) is the discrete part of the divergence,
- ( Bp^T p^{n+1} \approx \int{\Omega_{i,j,k}} \nabla p^{n+1} dV ) is the discrete gradient pressure operator,
- ( Mu^{(\rho)} \vec{u}^{n+1} \approx \int{\Omega{i,j,k}} \left( \frac{\rho^{n+1} \vec{u}^{n+1} + \overline{B}{\text{penu}} f(\vec{u}^{n+1})}{\Delta t} \right) dV ) is the mass matrix,
- ( Lu^{(\mu)} \vec{u}^{n+1} \approx -\int{\Sigma_{i,j,k}} \overline{\mathbf{T}^{n+1}} \cdot \vec{n} dS ) is the viscous stresses operator,
- ( Nu^{(\rho)} \vec{u}^{n+1} \approx \int{\Sigma_{i,j,k}} \rho \vec{u}^{n+1} \vec{u}^n \cdot \vec{n} dS ) is the discrete part of the convective term,
- ( \vec{f} \approx \int{\Omega{i,j,k}} \left( \frac{\rho^n \vec{u}^n + \overline{B}{\text{penu}} \vec{u}\infty + \vec{F}_s}{\Delta t} \right) dV ) is the second member associated with the speed.
!!! todo "Should nNavier be set to 1? Or only for 2 LS?"
!!! todo "explain interpolation uv to grid p with averaging coefficients"
why volume used for interpolating for iLS and difference in liquid heights for i ?
interpolating_coefficient_Navier_uv_grids_to_p_grid_volume
interpolating_coefficient_Navier_uv_grids_to_p_grid_height
FE_set_momentum_coupled2
pressure_projection!
In "An important issue is that the boundary condition for u§ must be consistent with Eq. (10), although at the time the boundary conditions are applied the function phi is not yet known and hence must be approximated"
set_Forward_Euler!
set_Crank_Nicolson!
One-fluid model
solve_one_fluid_NS!
Matrix
The matrix for the one-fluid model is defined here:
FE_set_momentum_coupled2_one_fluid
Convection
compute_fluxes(u, v, rho, dx_p, dy_p)
The interpolations for the viscosity are done with:
bilinear_interpolation
<figure>
<a name="levelset_doc"></a>
<img src="./assets/staggered_viscosity.svg" alt="Staggered grids" title="Staggered grids">
<figcaption>Staggered grids </figcaption>
</figure>
Surface tension
If num.surface_tension = 1, the surface tension is computed based on the levelset.
Based on the levelset
When the surface tension is computed based on the levelset, the levelset is firt smoothed with:
levelset_heavyside
compute_surface_tension_LS!
Based on the VOF
compute_surface_tension_VOF!
get_curvature
Velocity-pressure coupling
"
\begin{equation*}
\frac{\mathbf{u}^{*} - \mathbf{u}^{n}}{\Delta t} + \nabla q = \frac{\nu}{2} \nabla^2 (\mathbf{u}^{n} + \mathbf{u}^{*})
\end{equation*}
where $\nabla q$ is related to the pressure. The velocity satisfies $\nabla \cdot \mathbf{u}^{n+1} = 0$ and is given by
\begin{equation*}
\mathbf{u}^{n+1} = \mathbf{u}^{*} - \Delta t \nabla \phi^{n+1}
\end{equation*}
and the pressure is updated with
\begin{equation*}
p^{n+1/2} = q + L \phi^{n+1}
\end{equation*}
where $L$ is a linear differential operator.
Referring to Eqs. (8), (10), and (12), three combinations of $q$ and $L$ will be considered:
- a projection method similar to that of Bell, Colella, and Glaz, described by $q = p^{n-1/2}$ and $L = I$. This combination will be referred to as projection method I (PmI).
- a similar projection method that uses the improved pressure-update formula Eq. (13). This combination corresponds to $q = p^{n-1/2}$ and $L = I - \frac{\nu \Delta t}{2} \nabla^2$ and will be referred to as PmII.
- a projection method similar to that of Kim and Moin's, which corresponds to $q = 0$ and $L = I - \frac{\nu \Delta t}{2} \nabla^2$. This method will be referred to as PmIII.
"
"
\subsection*{Projection Methods with a Lagged Pressure Term}
Projection method I (PmI)
The method first described in this section is referred to as PmI. It is similar to the method developed by Bell \textit{et al.} \cite{bell1989second}, except in the treatment of the advective derivatives which are computed as in \cite{adams1959method} with a second-order Adams--Bashforth formula.
The first step of the projection method is found by solving
\begin{equation*}
\frac{\mathbf{u}^{*} - \mathbf{u}^{n}}{\Delta t} + \nabla p^{n-1/2} = -[(\mathbf{u} \cdot \nabla_h)\mathbf{u}]^{n+1/2} + \frac{\nu}{2} \nabla_h^2 (\mathbf{u}^{n} + \mathbf{u}^{*})
\end{equation*}
for the intermediate field $\mathbf{u}^{*}$ with boundary conditions
\begin{equation*}
\mathbf{u}^{*} = \mathbf{u}_b^{n+1}.
\end{equation*}
Next, $\mathbf{u}^{n+1}$ is recovered from the projection of $\mathbf{u}^{*}$ by solving
\begin{equation}
\Delta t \nabla_h^2 \phi^{n+1} = \nabla_h \cdot \mathbf{u}^{*} \quad \text{in } \Omega
\end{equation}
\begin{equation}
\hat{\mathbf{n}} \cdot \nabla_h \phi^{n+1} = 0 \quad \text{on } \partial \Omega
\end{equation}
and setting
\begin{equation}
\mathbf{u}^{n+1} = \mathbf{u}^{*} - \Delta t \nabla_h \phi^{n+1}.
\end{equation}
The new pressure is computed as in \cite{bell1989second, kim1985application} by
\begin{equation}
\nabla_h p^{n+1/2} = \nabla_h p^{n-1/2} + \nabla_h \phi^{n+1}.
\end{equation}
As discussed before, this formula is not consistent with a second-order discretization of the Navier--Stokes equations since, due to Eq. (72), the normal component of the pressure gradient will remain constant in time at the boundary.
Projection method II (PmII)
A second implementation of the method just described can be made by utilizing the correct pressure update given by Eq. (13). Specifically, Eqs. (70)--(72) are used in combination with
\begin{equation}
\nabla_h p^{n+1/2} = \nabla_h p^{n-1/2} + \nabla_h \phi^{n+1} - \frac{\nu \Delta t}{2} \nabla_h \nabla_h^2 \phi^{n+1}.
\end{equation}
<section class="algorithm">
<div class="counter"></div>
<div class="procedure procedure__start">
<span class="sc">PmII Algorithm</span>(<i>Initial velocity</i> <b>u</b><sup>n</sup>, <i>pressure</i> <b>p</b><sup>n-1/2</sup>, <i>viscosity</i> ν, <i>time step</i> Δt)
</div>
<div class="comment">
Update velocity <b>u</b><sup>n+1</sup> and pressure <b>p</b><sup>n+1/2</sup>
</div>
<div class="counter"></div>
<div class="procedure procedure__start">
<span class="sc">Step 1:</span> Solve for the intermediate velocity field <b>u</b><sup>*</sup>
</div>
<div class="counter"></div>
<div class="state ml1">
\[
\frac{\mathbf{u}^* - \mathbf{u}^n}{\Delta t} + \nabla p^{n-1/2} = -[(\mathbf{u} \cdot \nabla) \mathbf{u}]^{n+1/2} + \frac{\nu}{2} \nabla^2 (\mathbf{u}^* + \mathbf{u}^n)
\]
</div>
<div class="comment">
with boundary conditions <b>u</b><sup>*,γ</sup> = <b>u</b><sup>γ</sup>
</div>
<div class="counter"></div>
<div class="procedure procedure__start">
<span class="sc">Step 2:</span> Perform the projection
</div>
<div class="counter"></div>
<div class="state ml1">
\[
\Delta \phi = \frac{1}{\Delta t} \nabla \cdot \mathbf{u}
\]
</div>
<div class="comment">
With <span>\(\frac{\partial \phi}{\partial n}=0\)</span>
</div>
<div class="counter"></div>
<div class="state ml1">
\[
\mathbf{u}^{n+1} = \mathbf{u}^* - \Delta t \nabla \phi^{n+1}
\]
</div>
<div class="state ml1">
\[
\nabla \cdot \mathbf{u}^{n+1} = 0
\]
</div>
<!-- <div class="comment">
with boundary conditions consistent with <span>\(B(\mathbf{u}^*) = 0\)</span> and <span>\(\mathbf{u}^{n+1}|_{\partial \Omega} = \mathbf{u}_b^{n+1}\)</span>
</div> -->
<div class="counter"></div>
<div class="procedure procedure__start">
<span class="sc">Step 3:</span> Update the pressure
</div>
<div class="counter"></div>
<div class="state ml1">
\[
p^{n+1/2} = p^{n-1/2} + \phi^{n+1} - \frac{\nu \Delta t}{2} \nabla^2 \phi^{n+1}
\]
</div>
<div class="comment">
Re-impose BC or naturally satisfied if <span>\(p^0\)</span> satisfies pressure BC ? <span>\(\phi\)</span>: OK but <span>\(\Delta \phi\)</span> : <span>\(\nabla \nabla^2 \phi = \nabla \nabla \cdot \mathbf{u}=0\)</span>?
But with accumulation of numerical errors ?
</div>
<div class="counter"></div>
<div class="procedure procedure__end"></div>
</section>
!!! todo "Example algo html"
<section class="algorithm">
<div class="counter"></div>
<div class="procedure procedure__start">
<span class="sc">Euclid</span>(<i>a, b</i>)
</div>
<div class="comment">
The g.c.d. of <i>a</i> and <i>b</i>
</div>
<div class="counter"></div>
<div class="state ml1"><i>r ← a </i>mod<i> b</i></div>
<div class="counter"></div>
<div class="while while__start ml1">
<i>r ≠ </i>0
</div>
<div class="comment">
We have the answer if <i>r</i> is 0
</div>
<div class="counter"></div>
<div class="ml2"><i>a ← b</i></div>
<div class="counter"></div>
<div class="state ml2"><i>b ← r</i></div>
<div class="counter"></div>
<div class="ml2"><i>r ← a </i>mod<i> b</i></div>
<div class="counter"></div>
<div class="while while__end ml1"></div>
<div class="counter"></div>
<div class="ml1"><strong>return</strong> <i>b</i></div>
<div role="region" aria-label="comment" class="comment">
The gcd is <i>b</i>
</div>
<div role="region" aria-label="line" class="counter"></div>
<div class="procedure procedure__end"></div>
</section>
Projection method III (PmIII), a Projection Method without Pressure Gradient
It is similar to the method of Kim and Moin \cite{KimMoin}, but uses a different spatial discretization and a slightly different treatment of the boundary conditions. The momentum equation is discretized by
\begin{equation}
\frac{\mathbf{u}^* - \mathbf{u}^n}{\Delta t} = -[(\mathbf{u} \cdot \nabla_h) \mathbf{u}]^{n+1/2} + \frac{\nu}{2} \nabla_h^2 (\mathbf{u}^n + \mathbf{u}^*)
\end{equation}
and we first consider boundary conditions
\begin{align}
\hat{\mathbf{n}} \cdot \mathbf{u}^*|_{\partial \Omega} &= \hat{\mathbf{n}} \cdot \mathbf{u}_b^{n+1} \\
\hat{\mathbf{t}} \cdot \mathbf{u}^*|_{\partial \Omega} &= \hat{\mathbf{t}} \cdot (\mathbf{u}_b^{n+1} + \Delta t \nabla_h \phi^n)|_{\partial \Omega}.
\end{align}
As before, $\mathbf{u}^{n+1} = \mathbf{P}(\mathbf{u}^)$, i.e., $\mathbf{u}^{n+1} = \mathbf{u}^ - \Delta t \nabla_h \phi^{n+1}$, where $\phi^{n+1}$ satisfies Eqs. (71) and (72).
"
Either pressure gradient in prediction or remove component in velocity boundary conditions
Flower original projection
Problem: BC for velocity used as is but pressure gradient not in prediction, also Flower test case "channel.jl" has homogeneous Neumann BC at the injection whereas there is a pressure gradient in Poiseuille (in the direction of the flow)
With the parameter num.prediction, the following methods can be called:
- "Flower": original pressure-velocity coupling in Flower, the pressure gradient is not in the prediction and the boundary conditions are not modified accordingly
- "PmI": pressure gradient in prediction
- "PmII":
"PmIII": TODO, the pressure gradient is not in the prediction " In the present work, the method referred to as projection method II (PmII) by Brown et al. 2001, which ensures a second order discretization of the equations, is employed. The first step consists in updating an intermediate velocity field $\mathbf{u}^*$ using the momentum equation (referred to as prediction step), where the diffusive transport term is solved using a Crank--Nicolson scheme and the convective transport term using a second-order Adams--Bashforth integrator, "
\frac{\mathbf{u}^* - \mathbf{u}^n}{\tau} + [\mathbf{u} \cdot \nabla \mathbf{u}]^{n+1/2} = -\nabla p^{n-1/2} + \frac{\mu}{2} \Delta (\mathbf{u}^n + \mathbf{u}^*).
Here, $\tau$ refers to the time step and the superscript $n$ the time iteration. In general, the velocity field $\mathbf{u}^*$ is not divergence-free, and the next step enforces this condition by first solving the Poisson's equation
\tau \Delta \psi^{n+1} = \nabla \cdot \mathbf{u}^*,
for the intermediate pressure field $\psi$, prior to correcting the intermediate velocity field according to
\mathbf{u}^{n+1} = \mathbf{u}^* - \tau \nabla \psi^{n+1}.
Finally, the pressure is updated for the velocity prediction at the next time-step:
p^{n+1/2} = p^{n-1/2} + \psi^{n+1} - \frac{\tau \mu}{2} \Delta \psi^{n+1}.
"
The discrete cut-cell operators are used to discretize the pressure-projection method. The prediction step reads
\frac{u_{\alpha}^{\omega_{i, *}} - u_{\alpha}^{\omega_{i, n}}}{\tau} + \frac{3}{2} \text{conv}_{\alpha} \left( \boldsymbol{u}^{\omega_{i, n}}, \boldsymbol{u}^{\gamma, n} \right) - \frac{1}{2} \text{conv}_{\alpha} \left( \boldsymbol{u}^{\omega_{i, n-1}}, \boldsymbol{u}^{\gamma, n-1} \right) = -V_{\alpha} \left[ \text{grad} \left( p^{\omega_{i, n-1/2}}, p^{\gamma, n-1/2} \right) \right]_{\alpha} + \frac{1}{2 \text{Re}} \left[ \text{div}_{\alpha} \left( \text{grad}_{\alpha} \left( u_{\alpha}^{\omega_{i, *}}, u_{\alpha}^{\gamma, *} \right) \right) + \text{div}_{\alpha} \left( \text{grad}_{\alpha} \left( u_{\alpha}^{\omega_{i, n}}, u_{\alpha}^{\gamma, n} \right) \right) \right] \quad \forall \ \alpha \in \{x, y\}
The boundary conditions applicable to $\boldsymbol{u}^{*}$ are those of the velocity field at the next time step $\boldsymbol{u}^{n+1}$. The discrete Poisson’s equation results in
\text{div} \left( \text{grad} \left( \psi^{\omega_{i, n+1}}, \psi^{\gamma, n+1} \right) \right) = \frac{1}{\tau} \text{div} \left( \boldsymbol{u}^{\omega_{i, *}}, \boldsymbol{u}^{\gamma, *} \right)
The velocity field is ultimately corrected as
u_{\alpha}^{\omega_{i, n+1}} = u_{\alpha}^{\omega_{i, *}} - \tau V_{\alpha} \left[ \text{grad} \left( \psi^{\omega_{i, n+1}}, \psi^{\gamma, n+1} \right) \right]_{\alpha} \quad \forall \ \alpha \in \{x, y\}
and the pressure is finally updated as
p^{\omega_{i, n+1/2}} = p^{\omega_{i, n-1/2}} + \psi^{n+1} - \frac{\tau}{2 \text{Re}} \text{div} \left( \text{grad} \left( \psi^{\omega_{i, n+1}}, \psi^{\gamma, n+1} \right) \right)
where the last term ensures the second order accuracy of the pressure field.
As explained in Sec. 1.1.3, in order to obtain an energy preserving discretization of the incompressible Navier–Stokes equations, the pressure gradient must be equal to the negative transpose of the velocity divergence. This means that in the term $\left[ \text{grad} \left( p^{\omega{i, n-1/2}}, p^{\gamma, n-1/2} \right) \right]{\alpha}$, instead of using the definition given by Eq. (2.37), the transpose of the operators used in Eq. (2.39) must be employed. Regarding the boundary conditions for the pressure, to ensure the skew-symmetricity of the gradient and the divergence, opposite boundary conditions must be used. This means that whenever a Dirichlet boundary conditions is used in the velocity, a Neumann boundary condition is employed in the pressure and vice-versa.
"
Either pressure gradient in prediction or remove component in velocity boundary conditions
!!! todo "TODO"
BC document and test
"For moving boundaries the gradient of pressure was removed from the prediction because the simulations weren't very stable."
!!! info ""
To print information about the matrix in the velocity-pressure coupling: ``gp.ny +1`` in lexicographic for v and ``gp.ny`` in lexicographic for u
All the operators in the code are volume integrated,
so if you want to obtain the actual gradient or divergence you always have to divide by the volume
(In the finite volume method you have the volume that appears dividing)
For the viscous term, you cannot divide by opC_v.iMx_bd, the viscous term is collocated with the velocity component, and you're dividing by a staggered volume. You have to divide by the collocated volume. You have to divide by opC_p.iMy .
The parameter
num.prediction
is used to choose the prediction method.
!!! todo "Adams"
set_convection!
vector_convection!
" For moving boundaries the gradient of pressure was removed from the prediction because the simulations weren't very stable."
To print information about the matrix in the velocity-pressure coupling:
gp.ny +1 in lexicographic for v and gp.ny in lexicographic for u
All the operators in the code are volume integrated, so if you want to obtain the actual gradient or divergence you always have to divide by the volume (In the finite volume method you have the volume that appears dividing)
For the viscous term, you cannot divide by opC_v.iMx_bd, the viscous term is collocated with the velocity component, and you're dividing by a staggered volume. You have to divide by the collocated volume. You have to divide by opC_p.iMy .
p not v
Coupled
set_first_cells!
Convection
The convective transport of a vector field $\boldsymbol{q}$ by a flow velocity field $\boldsymbol{u}$, which can be rewritten in conservative form if $\boldsymbol{u}$ is divergence-free as
\begin{equation}
\left ( \boldsymbol{u} \cdot \nabla \right ) \boldsymbol{q} = \nabla \cdot \left ( \boldsymbol{q} \otimes \boldsymbol{u} \right ) - \left ( \nabla \cdot \boldsymbol{u} \right ) \boldsymbol{q} = \nabla \cdot \left ( \boldsymbol{q} \otimes \boldsymbol{u} \right ),
\end{equation}
param num.Cdivu
"
which in discrete form can be defined as
\begin{equation} \label{eq:conv1}
\begin{split}
\operatorname{conv} _ x \left ( \boldsymbol{u} ^ {\omega _ i}, \boldsymbol{u} ^ \gamma, \boldsymbol{q} ^ {\omega _ i}, \boldsymbol{q} ^ \gamma \right ) &= C _ x\left ( \boldsymbol{u} ^ {\omega _ i} \right ) q _ x ^ {\omega _ i} - \dfrac{K _ x \left ( \boldsymbol{u} ^ \gamma \right ) q _ x ^ {\omega _ i} + K _ x \left ( \boldsymbol{q} ^ \gamma \right ) u _ x ^ {\omega _ i}}{2}, \\
\operatorname{conv} _ y \left ( \boldsymbol{u} ^ {\omega _ i}, \boldsymbol{u} ^ \gamma, \boldsymbol{q} ^ {\omega _ i}, \boldsymbol{q} ^ \gamma \right ) &= C _ y \left ( \boldsymbol{u} ^ {\omega _ i} \right ) q _ y ^ {\omega _ i} - \dfrac{K _ y \left ( \boldsymbol{u} ^ \gamma \right ) q _ y ^ {\omega _ i} + K _ y \left ( \boldsymbol{q} ^ \gamma \right ) u _ y ^ {\omega _ i}}{2},
\end{split}
\end{equation}
where $\boldsymbol{q} ^ \omega$, $\boldsymbol{q} ^ \gamma$, $\boldsymbol{u} ^ \omega$ and $\boldsymbol{u} ^ \gamma$ denote face-centered discrete vector fields and
\begin{equation}
\begin{split}
C _ x \left ( \boldsymbol{u} ^ {\omega _ i} \right ) &= D _ {x, x} ^ + \operatorname{diag} \left ( S _ {x, x} ^ - A _ x u ^ {\omega _ i} _ x \right ) S _ {x, x} ^ - + D _ {x, y} ^ + \operatorname{diag} \left ( S _ {x, x} ^ - A _ y u ^ {\omega _ i} _ y \right ) S _ {x, y} ^ -, \\
K _ x \left ( \boldsymbol{u} ^ \gamma \right ) &= \operatorname{diag} \left ( S ^ + _ {x, x} H ^ \top \boldsymbol{u} ^ \gamma \right ), \\
C _ y \left ( \boldsymbol{u} ^ {\omega _ i} \right ) &= D _ {y, x} ^ + \operatorname{diag} \left ( S _ {y, y} ^ - A _ x u ^ {\omega _ i} _ x \right ) S _ {y, x} ^ - + D _ {y, y} ^ + \operatorname{diag} \left ( S _ {y, y} ^ - A _ y u ^ {\omega _ i} _ y \right ) S _ {y, y} ^ -, \\
K _ y \left ( \boldsymbol{u} ^ \gamma \right ) &= \operatorname{diag} \left ( S ^ + _ {y, y} H ^ \top \boldsymbol{u} ^ \gamma \right ).
\end{split}
\end{equation}
The definition of the convective term is based on the skew-symmetric convective operator given in \cite{Morinishi1998} for finite differences, which has been modified to work with the cut-cell method. The operators acting on the bulk field are skew-symmetric, as required to have a conservative discretization of the Navier-Stokes equation. Eq. \eqref{eq:conv1} can be rearranged to show the skew-symmetric operator
\begin{equation}
\begin{split}
\operatorname{conv} _ x \left ( \boldsymbol{u} ^ {\omega _ i}, \boldsymbol{u} ^ \gamma, \boldsymbol{q} ^ {\omega _ i}, \boldsymbol{q} ^ \gamma \right ) &= \underbrace{\left [ C _ x \left ( \boldsymbol{u} ^ {\omega _ i} \right ) - \dfrac{1}{2} K _ x \left ( \boldsymbol{u} ^ \gamma \right ) \right ]} _ \text{Skew-symmetric} v _ x ^ {\omega _ i} - \dfrac{1}{2} K _ x \left ( \boldsymbol{q} ^ \gamma \right ) u _ x ^ {\omega _ i}, \\
\operatorname{conv} _ y \left ( \boldsymbol{u} ^ {\omega _ i}, \boldsymbol{u} ^ \gamma, \boldsymbol{q} ^ {\omega _ i}, \boldsymbol{q} ^ \gamma \right ) &= \underbrace{\left [ C _ y \left ( \boldsymbol{u} ^ {\omega _ i} \right )- \dfrac{1}{2} K _ y \left ( \boldsymbol{u} ^ \gamma \right ) \right ]} _ \text{Skew-symmetric} v _ y ^ {\omega _ i} - \dfrac{1}{2} K _ y \left ( \boldsymbol{q} ^ \gamma \right ) u _ y ^ {\omega _ i}.
\end{split}
\end{equation}
In the following, whenever $\boldsymbol{q} = \boldsymbol{u}$, the convective term will be referred as $\operatorname{conv} _ \alpha \left ( \boldsymbol{u} ^ \omega, \boldsymbol{u} ^ \gamma \right )\ \forall\ \alpha \in { x, y }$. The discrete convective operator also exhibits the free-streaming condition, hence for $\boldsymbol{u} ^ \omega = \boldsymbol{u} ^ \gamma = \boldsymbol{c}$, $\operatorname{conv} _ x \left ( \boldsymbol{c}, \boldsymbol{c} \right ) = 0$ and $\operatorname{conv} _ y \left ( \boldsymbol{c}, \boldsymbol{c} \right ) = 0$. "
!!! todo "multiple LS"
\begin{equation}
\text{conv}_\alpha \left( u^\omega, u^\gamma, v^\omega, v^\gamma \right)= C_x \left( \mathbf{u}^\omega \right) v_x^\omega - \frac{1}{2} K_x\left( \right)
\end{equation}
\begin{equation}
\text{conv}_\alpha \left( \mathbf{u}^\omega, \mathbf{u}^\gamma, T^\omega, T^\gamma \right)= \sum_\alpha \left\{ \frac{\delta A_\alpha u_\alpha^\omega \overline{T^\omega}^\alpha }{\delta \xi_\alpha} + \left[ \frac{\delta\left( \overline{B_\alpha}^\alpha-A_\alpha\right) \mathbf{u}_\alpha^\gamma}{\delta \xi_\alpha} -\overline{\frac{\delta B_\alpha \mathbf{u}_\alpha^\gamma}{\delta \xi_\alpha}}^\alpha\right]\frac{T^\omega + T^\gamma}{2}\right\}
\end{equation}
In set_convection! which uses vector_convection!
Du_y[1,:] .= vecb_B(uD,grid_u) + dt* grd_x[1,:]
Extract from (Rodriguez et al. 2024):
" This section presents the proposed discretization of the incompressible Navier-Stokes equations for an isotropic Newtonian fluid
\begin{cases}
\frac{\partial u}{\partial t} + \left ( u \cdot \nabla \right ) u & = -\dfrac{\nabla p}{\rho} + \dfrac{1}{\mathrm{Re}} \nabla \cdot \nabla u \\
\nabla \cdot u & = 0
\label{eq:ns}
\end{cases}
where $u$ and $p$ respectively denote the fluid's velocity and pressure fields, $\rho$ its constant density and $\mathrm{Re}$ denotes the Reynolds number. This equation can be rewritten in semi-discrete form using the linear operators previously defined, yielding the momentum equation in direction $\alpha$
\begin{equation} \label{eq:nsSemi} V _ \alpha \frac{\mathrm{d} u _ \alpha ^ \omega}{\mathrm{d} t} + \mathrm{conv} _ \alpha \left ( u ^ \omega, u ^ \gamma \right ) = - \dfrac{1}{\rho} \left [ \mathrm{grad} \left ( p ^ \omega, p ^ \gamma \right ) \right ] _ \alpha + \dfrac{1}{\mathrm{Re}} \mathrm{div} _ \alpha \left ( \mathrm{grad} _ \alpha \left ( u _ \alpha ^ \omega, u _ \alpha ^ \gamma \right ) \right ), \end{equation}
and the divergence-free condition \begin{equation} \label{eq:cont}
\mathrm{div} \left ( u ^ \omega, u ^ \gamma \right ) = 0,
\end{equation} where $u ^ \omega$ and $u ^ \gamma$ are two face-centered vector fields and $p ^ \omega$ and $p ^ \gamma$ are two cell-centered scalar fields. [...] The method referred to as projection method II (PmII) by Brown et al. 2001 is employed. It ensures a second order discretization of the equations.
In this projection method, the convective term is discretized using the explicit second-order Adams-Bashforth scheme:
!!! todo "boundary conditions"
and the viscous term is discretized using the implicit Crank-Nicolson scheme. The first step of the method consists of obtaining an intermediate velocity field $u ^ \star$ by solving:
\begin{multline} \label{eq:pmPred}
V _ \alpha \frac{u _ \alpha ^ {\omega, \star} - u _ \alpha ^ {\omega, n}}{\tau} + \frac{3}{2} \mathrm{conv} _ \alpha \left ( u ^ {\omega, n}, u ^ {\gamma, n} \right ) \\
- \frac{1}{2} \mathrm{conv} _ \alpha \left ( u ^ {\omega, n - 1}, u ^ {\gamma, n - 1} \right ) = - \dfrac{1}{\rho} \left [ \mathrm{grad} \left ( p ^ {\omega, n - 1 / 2}, p ^ {\gamma, n - 1 / 2} \right ) \right ] _ \alpha \\
+ \dfrac{1}{2 \mathrm{Re}} \left [ \mathrm{div} _ \alpha \left ( \mathrm{grad} _ \alpha \left ( u _ \alpha ^ {\omega, \star}, u _ \alpha ^ {\gamma, \star} \right ) \right ) + \mathrm{div} _ \alpha \left ( \mathrm{grad} _ \alpha \left ( u _ \alpha ^ {\omega, n}, u _ \alpha ^ {\gamma, n} \right ) \right ) \right ]\ \forall\ \alpha \in \{ x, y \},
\end{multline}
where $\tau$ denotes the time step and the superscript $n$ the iteration number. The boundary conditions applicable to $u ^ \star$ (the predicted velocity field) are those of the velocity field at the next time step.
The Dirichlet boundary conditions applied for example to impose the no-slip condition or the inflow conditions read \begin{equation}
u _ \alpha ^ {\gamma, \star} = u _ {\alpha, \mathrm{ref}}, \quad \forall \alpha \in \{ x, y \},
\end{equation} where $u _ {\alpha, \mathrm{ref}}$ is the value of the velocity to be imposed at the boundary in the $\alpha$-direction. The homogeneous Neumann boundary condition applied to impose outflow conditions is given by \begin{equation}
-\mathrm{div} _ \alpha \left ( 0, \, \mathrm{grad} _ \alpha \left ( u _ \alpha ^ {\omega, \star}, u _ \alpha ^ {\gamma, \star} \right ) \right ) = 0, \quad \forall \alpha \in \{ x, y \},
\end{equation}
In the projection step, the velocity field is updated by projecting $u ^ \star$ using the intermediate pressure field $\psi ^ {n + 1}$, which is obtained by solving the following Poisson equation \begin{equation} \label{eq:pmPoisson} \mathrm{div} \left ( \mathrm{grad} \left ( \psi ^ {\omega, n + 1}, \psi ^ {\gamma, n + 1} \right ) \right ) = \dfrac{\rho}{\tau}\mathrm{div} \left ( u ^ {\omega, \star}, u ^ {\gamma, \star} \right ). \end{equation} Whenever no-slip or inflow conditions are to be imposed, a homogeneous Neumann boundary condition is applied to the intermediate pressure field following \begin{equation}
-\mathrm{div} \left ( 0, \, \mathrm{grad} \left ( \psi ^ {\omega, n + 1}, \psi ^ {\gamma, n + 1} \right ) \right ) = 0.
\end{equation} With outflow boundary conditions, a Dirichlet boundary condition is imposed following \begin{equation}
\psi ^ {\gamma, n + 1} = \psi _ {\mathrm{ref}},
\end{equation}
The velocity field is ultimately corrected as \begin{equation} u_\alpha ^ {\omega, n + 1} = u ^ {\omega, \star} _\alpha - \dfrac{\tau}{\rho} \left [ \mathrm{grad} \left ( \psi ^ {\omega, n + 1}, \psi ^ {\gamma, n + 1} \right ) \right ] _ \alpha\ \forall\ \alpha \in { x, y }. \label{eq:pm_projection} \end{equation} and the pressure is finally updated as \begin{equation} p ^ {\omega, n + 1 / 2} = p ^ {\omega, n - 1 / 2} + \psi ^ {n + 1} - \frac{\tau}{2 \mathrm{Re}} \mathrm{div} \left ( \mathrm{grad} \left ( \psi ^ {\omega, n + 1}, \psi ^ {\gamma, n + 1} \right ) \right ), \label{eq:pm_update} \end{equation} where the last term ensures the second order accuracy of the pressure field.
"
Boundary conditions
Boundary conditions for border
The idea there is to have at the corners the contributions from both borders. bcTx is only used for derivatives in x and bcTy for derivatives in y, so it doesn't matter that bcTx is not filled at the top and bottom and the same for bcTy. Otherwise we would have to select one of the contributions.
set_borders!
Interface definition
Liquid cells have an isovalue of 0, solid cells fave an isovalue of 15.0 and mixed cells have an other value (cf marching squares).
Convention, orientation of the surface
The signed distance is positive in the liquid and negative in the bubble. The normal is oriented towards the liquid.
update_all_ls_data
marching_squares!
get_cells_indices
Interface transport (with Levelset)
IIOE_normal!
IIOE_normal_indices_2!
indices_extension
inflow_outflow
field_extension!
S2IIOE!
Operators
Gradient of a scalar field
!!! todo "u v nodes at wall"
At the wall, the u and v nodes are associated with the midpoints of the border but are in fact in the centroid of the half control volume.
!!! todo "difference grad p cap and u v cap"
difference: half control volumes ...
compute_grad_T_x_T_y_array_u_v_capacities!
compute_grad_T_x_T_y_array!
compute_grad_T_x_array!
compute_grad_T_y_array!
compute_grad_p!
normal_gradient
Laplacian
set_matrices!
laplacian
laplacian_bc
There are several ways to present the method: G and H, Morinishi's notations and a geometric formulation: (Rodriguez et al. 2024).
diamond
Functions for divergence
-G^T= [ B_xD_x^+ B_yD_y^+ ]
with
\begin{equation}
( D ^ - _ x ) ^ \top = -D ^ + _ x \quad \mathrm{and} \quad ( S ^ - _ x ) ^ \top = S ^ + _ x,
\label{eq:multitrans}
\end{equation}
\Delta ^ + = \left [ \begin{array}{rrrrr}
-1 & 1 & & & \\
& -1 & 1 & & \\
& & \ddots & \ddots & \\
& & & -1 & 1 \\
& & & & 0
\end{array} \right ]
divergence_B!
Then, for the gradient, B_x corresponds to -BxT' i.e. -(-G^T)^T=G
Updating the operator from Levelset
!!! todo "Order of updates"
At every iteration, [`update_all_ls_data`](@ref) is called twice, once inside run.jl and another one
(if there's advection of the levelset) inside [`set_heat!`](@ref).
The difference between both is a flag as last argument, inside run.jl is implicitly defined as true
and inside [`set_heat!`](@ref) is false. If you're calling your version of [`set_heat!`](@ref) several times,
then you're calling the version with the flag set to false, but for the convective term it has to be set to true.
The flag=true, the capacities are set for the convection, the flag=false they are set for the other operators
In Flower.jl:
DiscreteOperators
In scalar_transport, this corresponds to:
# Interior BC
A[sb,1:ni] = b * (HxT[iLS] * iMx * Bx .+ HyT[iLS] * iMy * By)
# Contribution to Neumann BC from other boundaries
for i in 1:num.nLS
if i != iLS
A[sb,i*ni+1:(i+1)*ni] = b * (HxT[iLS] * iMx * Hx[i] .+ HyT[iLS] * iMy * Hy[i]) #-b in Poisson
end
end
A[sb,sb] = pad(b * (HxT[iLS] * iMx * Hx[iLS] .+ HyT[iLS] * iMy * Hy[iLS]) .- op.χ[iLS] * a1)
#TODO no need for fs_mat
A[sb,end-nb+1:end] = b * (HxT[iLS] * op.iMx_b * op.Hx_b .+ HyT[iLS] * op.iMy_b * op.Hy_b)
Interpolation
Acpm
interpolate_grid_liquid!
interpolated_temperature
static_stencil
Interpolating to display solid and liquid
@views fwd.trans_scal[1,:,:,iscal] .= phL.trans_scal[:,:,iscal].*LS[end].geoL.cap[:,:,5] .+ phS.trans_scal[:,:,iscal].*LS[end].geoS.cap[:,:,5]
Convection
Capacities
The mass matrix:
M: face-centered mass matrices, diagonal with coefficients V_\alpha (the volume of the staggered control volumes)
is stored in the struct Operators and defined in set_matrices! in set_cutcell_matrices!:
M.diag .= vec(geo[end].dcap[:,:,5])
Example:
M_vL = Diagonal(fzeros(grid_v))
Example: V(1,1) ⋅ ⋅
⋅ V(...,...) ⋅
⋅ ⋅ V
New capacities from x^2
bc_matrix_borders!
capacities
mat_assign_T!
mass_matrix_borders!
periodic_bcs_borders!
set_boundary_indicator!
!!! todo "Rx Ry"
periodic_bcs_R!
B vs H
\chi
for iLS in 1:num.nLS
χx = (geo.dcap[:,:,3] .- geo.dcap[:,:,1]) .^ 2
χy = (geo.dcap[:,:,4] .- geo.dcap[:,:,2]) .^ 2
χ[iLS].diag .= sqrt.(vec(χx .+ χy))
end
Small cells
Small cells are cut with num.\epsilon and num.\epsilon_{wall}. This is used to prevent abnormally high values
The mid points are not exactly on the interface when the epsilon parameter is not null:
num.\epsilon = 0
num.\epsilon = 5%
Solver
For LS: gmres
grid.LS[iLS].u .= reshape(gmres(grid.LS[iLS].A, grid.LS[iLS].B * vec(grid.LS[iLS].u) .+ rhs_LS), grid)
Elsewhere: Julia's solver: /
- Diagonal scaling does not help (direct resolution)
- MUMPS similar
Changing
- diag + eps
- eps small cell
with diagonal scaling (not useful for direct resolution?): norm2(Ax-b)/norm2(b)=21.5550725612-12
System
A[end-nb+testn,end-nb+1:end] is the contribution to the gradient in the normal direction from the border, and
A[end-nb+testn,ni+1:2*ni] from the interface
If the interface is nearly perpendicular to the border, there's no much contribution from the interface to the gradient in the x-direction.
Cutcell
ilp2cap
set_cutcell_matrices!
set_other_cutcell_matrices!
Sets Gx_b for divergence of velocity
bc_matrix_borders!(grid, Hx_u, Hy_v, Hx_p, Hy_p, dcap)
called in
set_border_matrices!
Interface transport
select_advection!
Levelset
!!! todo "TODO"
The levelset is not extended for the evaluation of capacities. At the borders, the interpolation is shifted by one cell. You can see this in [marching_squares!](@ref) (the ``II\_0`` variable).
Ghost cells
allocate_ghost_matrices_2
init_ghost_neumann_2
Temporal scheme
!!! todo "Time-step
change coefficients for extrapolation in case of adaptive time-step
Crank-Nicolson
Forward-Euler
Initialization
Interface position for storage
"Moreover, we define the interface centroid as the mid point of the segment crossing the cell which will be used in the computation of the Stefan condition." Fullana (2022)
Store segments in 2D to save up space
#Initialize liquid phase
x_centroid = gp.x .+ getproperty.(gp.LS[1].geoL.centroid, :x) .* gp.dx
y_centroid = gp.y .+ getproperty.(gp.LS[1].geoL.centroid, :y) .* gp.dy
x_bc = gp.x .+ getproperty.(gp.LS[1].mid_point, :x) .* gp.dx
y_bc = gp.y .+ getproperty.(gp.LS[1].mid_point, :y) .* gp.dy
#Initialize bulk value
vec1(phL.TD,gp) .= vec(ftest.(x_centroid,y_centroid))
#Initialize interfacial value
vec2(phL.TD,gp) .= vec(ftest.(x_bc,y_bc))
Convection
!!! todo "Advection equation of a vector-valued field 2.5.1 Discrete operators for staggered quantities" Rodriguez (2024)
Boundaries
Boundaries
BoundariesInt
Masks
compute_mask_1D!
Small cells
Small cells introduce errors, the parameter \epsilon is used to delete small cells. Merging the cells could help with this problem, though it would involve modifying the capacities and the structure of the cut-cell operators.
Need concentration near wall for conductivity for Poisson equation
take the value of the concentration in the bulk field instead
BC_phi_ele.left.val .= -i_butler./elec_cond[:,1]
for the conductivity used in the Poisson equation, instead of doing :
BC_phi_ele.left.val .= -i_butler./vecb_L(elec_condD, grid)
The conductivity is computed from the scalar. Taking the value in the bulk field is recommended as long as cell merging is not implemented. Due to small cells, we may have slivers/small cells at the left wall, then the divergence term is small, which produces a higher concentration in front of the contact line.
Solver
\section{Solvers}
In Flower for LS: gmres
grid.LS[iLS].u .= reshape(gmres(grid.LS[iLS].A, grid.LS[iLS].B * vec(grid.LS[iLS].u) .+ rhs_LS), grid)
Elsewhere: Julia's solver: /
- Diagonal scaling does not help (direct resolution)
- MUMPS similar
Changing
- diag + eps
- eps small cell
with diagonal scaling (not useful for direct resolution?): norm2(Ax-b)/norm2(b)=21.5550725612-12
Timestep
adapt_timestep!
Possible improvements
Restart (TODO restart with PDI) cf hello_access.jl
if sim.restart == 1: cf hello_access.jl endInterface length
Remove
replace!(N, NaN=>0.0)imposed constant velocity field : after 100 it : error after scalar transport 3.88e-14
pressure BC during prediction
small cells: linear algebra, time step ...
Epsilon to prevent NaN, with eps : coefficients not exact
sometimes Julia reports out of memory
Compilation time: upcoming development: \url{https://jbytecode.github.io/juliac/}
Store segments in 2D to save up space
Memory, allocations
check conservation of variables after n iterations
stencils
parallelisation
check BC wall pressure (wall no slip : Neumann it seems)
List of variables which can be removed:
- emptied
- cap (already dcap)
- ...
List all variables
- bulk+ interface \rightarrow *2
- phL.u phL.uD duplicate for bulk
- u, v, T, LS, moments (heights,...), phi elec, i current, scalars ...
- grids
- liquid/solid
- successive substitution reuse LU decomposition: use "factorize(A)"
- BitArray to mask values (small cells or solid) to compute min, mask, ...
- integral quantities
- test radial + BC wall
- manufactured solution for diffusion
- sparse array for mask 1D
!!! todo "Document capacities p, u, v"
```julia
# opC_p = op.opC_p
# ∇ϕ_x = opC_p.iMx * opC_p.Bx * vec1(pD,grid) .+ opC_p.iMx_b * opC_p.Hx_b * vecb(pD,grid)
# ∇ϕ_y = opC_p.iMy * opC_p.By * vec1(pD,grid) .+ opC_p.iMy_b * opC_p.Hy_b * vecb(pD,grid)
# for iLS in 1:num.nLS
# ∇ϕ_x .+= opC_p.iMx * opC_p.Hx[iLS] * veci(pD,grid,iLS+1)
# ∇ϕ_y .+= opC_p.iMy * opC_p.Hy[iLS] * veci(pD,grid,iLS+1)
# end
∇ϕ_x = opC_u.AxT * opC_u.Rx * vec1(pD,grid) .+ opC_u.Gx_b * vecb(pD,grid)
∇ϕ_y = opC_v.AyT * opC_v.Ry * vec1(pD,grid) .+ opC_v.Gy_b * vecb(pD,grid)
for iLS in 1:num.nLS
∇ϕ_x .+= opC_u.Gx[iLS] * veci(pD,grid,iLS+1)
∇ϕ_y .+= opC_v.Gy[iLS] * veci(pD,grid,iLS+1)
end
```
Epsilon/handling small or cancelled cells
num.\epsilonis uded for small cellsweights for cancelled cells
inv_weight_clip
Sources of errors
inv_weight_clip- eps added to diagonal (
A[i,i] += 1e-10) - epsilon for small cells (to prevent NaN) (
num.\epsilon), the coefficients are no longer exact in simple cases (like Laplacian: 1 1 -4 1 1) - Approximation of gradient inside a cell (
q^\omega \approx q^\gamma) - approximation of interface
Tools for grids
create_2D_grid_x
create_2D_grid_y
IO/post-processing
convert_interfacial_D_to_segments
The output files are generated with PDI. The on-line PDI documentation is available at PDI documentation.
TODO
Need details on dimensions: construction not explained (sparsity)
Scripts
- For diagrams mesh_square_bubble_wall.py