Commit 87d5c1f6 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts

pffff

parent 895e75b8
This diff is collapsed.
\section{Finite Element Formulation}
We present a finite element formulation for the simulation of a viscous incompressible flow that is an intermediate between a free flow (Stokes flow) and porous media flow (Darcy flow).
Consider our problem in a domain $\Omega$ with homogeneous Dirichlet boundary conditions
on its boundary $\Gamma$:
\begin{eqnarray}
\nabla p + {\mu \over K} \mvu - \bar{\mu} \nabla^2 \mvu &=& {\mathbf 0}~~\mbox{in}~~ \Omega\nonumber \\
-\nabla \cdot (\mvu + \mvu_s) &=& 0~~\mbox{in}~~ \Omega\nonumber \\
\mvu + \mvu_s& =& {\mathbf 0} ~~\mbox{on}~~ \Gamma\nonumber
\end{eqnarray}
We use the following notations:
The variational formulation of the problem is the following. Find $(\mvu,p) \in [H^1_0(\Omega)]^3 \times L^2_0(\Omega)$ such that
\begin{eqnarray}
\label{eq:var}
&&\left(\nabla p, \mvv\right)_{\Omega} +
\left(\mvu, \nabla q\right) _{\Omega} +
\left(\bar{\mu} \nabla \mvu, \nabla \mvv \right) _{\Omega} +
\left({\mu \over K} \mvu, \mvv \right) _{\Omega}
= \left(\nabla \cdot \mvu_s, q\right) _{\Omega}\nonumber \\
&&\forall (\mvv,q) \in [H^1_0(\Omega)]^3 \times L^2_0(\Omega).
\end{eqnarray}
Equation \eqref{eq:var} can be written in the usual abstract form: find
$(\mvu,p) \in [H^1_0(\Omega)]^3 \times L^2_0(\Omega)$ such that
$$B(\mvu,p;\mvv,q) = L(\mvv,q),~~~~\forall (\mvv,q) \in [H^1_0(\Omega)]^3 \times L^2_0(\Omega).$$
This continuous problem has a unique solution.
We assume a partitioning ${\mathcal T}$ of the domain $\Omega$. We denote $T$ an element of the partitioning,
and its size of $h_T$. Let us first define standard continuous finite element spaces
$${\mathbf V}_h =\left\{\mvv \in [H^1_0(\Omega) \cap C(\Omega)]^3 ~|~\mvv|_T \in [P_k(K)]^3,~~\forall T \in {\mathcal T} \right\},$$
$${ Q}_h =\left\{p \in L^2(\Omega) \cap C(\Omega)~|~p|_T \in P_k(K),~~\forall T \in {\mathcal T} \right\}.$$
We use a stabilized method for ensuring the stability of the discrete formulation: find $(\mvu_h, p_h) \in {\mathbf V}_h \times Q_h$
such that
$$B_h(\mvu_h,p_h;\mvv,q)=L_h(\mvv,q)~~ \forall (\mvv_h, q_h) \in {\mathbf V}_h \times Q_h.$$
For that, we define a the \emph{stabilized bilinear form} \cite {Juntunen}
\begin{eqnarray}
&&B_h(\mvu_h,p_h;\mvv,q) = B(\mvu_h,p_h;\mvv,q) - \nonumber \\
&&\alpha \sum_{T \in {\mathcal T}} {h_T^2 \over K+h_T^2}
\left (
\nabla p_h + {\mu \over K} \mvu_h - \mu \nabla^2 \mvu_h,
{K\over\mu} \nabla q + \mvv - K \nabla^2 \mvv
\right)_T
\end{eqnarray}
and
\begin{eqnarray}
L_h(\mvv,q) = L(\mvv,q) -
\alpha \sum_{T \in {\mathcal T}} {h_T^2 \over K+h_T^2}
\left (\nabla \cdot \mvu_s,q
\right)_T
\end{eqnarray}
Coefficient $\alpha$ has to be positive.
If linear finite elements are used, this simplifies to
\begin{eqnarray}
&&B_h(\mvu_h,p_h;\mvv,q) = B(\mvu_h,p_h;\mvv,q) - \nonumber \\
&&\alpha \sum_{T \in {\mathcal T}} {h_T^2 \over K+h_T^2}
\left[
\left ({K\over\mu} \nabla p_h,\nabla q\right)_T +
\left ({\mu\over K} \mvu_h,\mvv\right)_T +
\left (\nabla p_h,\mvv\right)_T +
\left (\mvu_h, \nabla q\right)_T
\right].\nonumber
\end{eqnarray}
\subsection {FE formulation ($v^β = 0$)}
\begin{eqnarray*}
«φv_\alpha ·ṉ» - ‹\nabla φ·v_\alpha&=& 0 \\
«\mu ψn · \nabla v_\alpha » -‹\mu \big(\nabla ψ·\nabla v_\alpha + ψ\nabla v^\alpha ·\nabla c_\alpha \big)› - ‹ψc_\alpha \nabla p^\alpha › -‹\frac{\mu (1-c_\alpha}{kc²_\alpha}v_\alpha&=& 0
\end{eqnarray*}
where
\[
\nabla v^\alpha = \frac{\nabla v_\alpha}{c_\alpha} - \frac{v_\alpha \nabla c_\alpha}{_\alpha}
\]
\section{Avancement des particules}
\section{Couplage temporel}
\subsection{Avancement des particules}
Pour le moment, on résout en déplacement (si j'ai bien compris ce qu'on entendait par là). C'est à dire que:
\begin{align*}
a^i &= a(x^i, v^i)\\
......@@ -22,7 +23,7 @@ x^{i+1} &= x^i + v^{i+1}\, Δt
}
\section{Accélération partiellement implicite}
\subsection{Accélération partiellement implicite}
\begin{align*}
\acc(x, v) &= \af(x, v) + \aext\\
\af(x, v) &= -β∇p -α(v-\vf)\,(v-\vf)
......@@ -46,7 +47,7 @@ telle que $\tia(x, v, v) = a(x, v)$ et
$v$ est traité implicitement dans l'approximation $\aim$ et explicitement dans la correction $\aex$. Comme $\tiv$, on prend la vitesse au début du pas de temps.
\nb{\item si on applique RK à $x$, il faut aussi introduire un $\tilde x$ pour la partie implicite}
\section{Formalisme IMEX}
\subsection{Formalisme IMEX}
Je suis ici les notations d'\cite{ascher_implicit-explicit_1995}: $u$ est le vecteur inconnu, $g$ la partie implicite et $f$ la partie explicite.
Pour un pas de temps, on a:
\begin{description}
......@@ -64,10 +65,10 @@ $ u^n_i = u^n + Δt \sum_{j=1}^i\left(a_{i,j}\,K_j+\hat a_{i+1,j}\,\hat K_j\righ
NB: Lorsque $\hat b_{s+1} = 0$, l'évaluation de $\hat K_{s+1}$ n'est pas nécessaire.
\section{Application}
\subsection{Application}
Dans notre cas, vu que la linéarisation se fait autour de $v^n$ et que $x$ ne rentre pas dans le schéma RK, on a $\hat K_1 = 0$, ce qui réduit de 1 le nombre d'évaluation de la partie explicite.
Le $u^{n+1}$ désigne en fait $\vfree^{n+1}$ et la seule résolution des contacts a lieu, à la fin de l'itération.
\subsection{Forward-backward Euler (1,1,1)}
\subsubsection{Forward-backward Euler (1,1,1)}
\[
\text{implicite: }\begin{array}{c|cc}
0 & 0 & 0 \\
......@@ -96,7 +97,7 @@ On résout le problème de Stokes en utilisant:
= \frac{-β∇p-α(v^n-\vf)(v^n - \vf + Δt\,\aext)}{1 + Δt\,α(v^n-\vf)}\]
Et on obtient la vitesse libre:
\[\vfree^{n+1} = v^{n} + Δt\, (\af + \aext)\]
\subsection{Forward-backward Euler (1,2,1)}
\subsubsection{Forward-backward Euler (1,2,1)}
\[
\text{implicite: }\begin{array}{c|cc}
0 & 0 & 0 \\
......@@ -130,7 +131,7 @@ où $\haf{2}$ est obtenu en résolvant le problème de Stokes en utilisant $v =
\[\haf{2} = -β∇p-α(v_2-\vf)(v_2 - \vf)\]
On trouve alors $\vfree^{n+1}$ à partir du {$\haf{2}$} explicite :
\[\vfree^{n+1} = v^n + Δt\,({\af}_2 + \aext + \haf{2} - {\af}_2) = v^n + Δt\,(\haf{2} + \aex)\]
\subsection{Implicit-Explicit midpoint (1,2,2)}
\subsubsection{Implicit-Explicit midpoint (1,2,2)}
\[
\text{implicite: }\begin{array}{c|cc}
0 & 0 & 0 \\
......@@ -152,5 +153,5 @@ Même coût que le précédent mais second ordre de précision.
\end{description}
Les calculs sont identiques, il faut juste utiliser $\frac{Δt}{2}$ au lieu de $Δt$ pour le calcul de ${\af}_2$ et $v_2$.
\subsection{plus haut ordre de précision}
\subsubsection{plus haut ordre de précision}
On va déjà essayer les précédents ...
......@@ -56,10 +56,12 @@
\newpage
\input{equations.tex}
\input{particle-fluid-interaction.tex}
\input{fe.tex}
\input{stability.tex}
\input{stability2.tex}
\input{imex.tex}
\input{scontact.tex}
\input{validation.tex}
\newpage
\bibliography{zotero}
\end{document}
\section{Particle-fluid interaction forces}
Assume a porous medium formed of particles of average diameter $d_p$ and a fluid
characterized by its dynamic viscosity $\mu$ [$\mbox{Pa} \cdot \mbox{sec}$].
Assume that the representative elementary volume (REV) is caracterized by a porosity $\phi$ defined as
$$
\phi = {\mbox{Volume of fluid in the REV} \over \mbox{Volume of the REV}}.
$$
Let us call $\mvu = (u_1,u_2,u_3)$ and $p$ the average fluid velocity and the average fluid pressure over a REV.
Assume the pressure gradient to be proportional both to the velocity and its Laplacian (Brinkman \cite{brinkman}, Forchheimer \cite{forchheimer}):
$$\nabla p = -{\mu \over K} \mvu -{C_E \over {\sqrt{K}}} \|\mvu\| \mvu + \bar{\mu} \nabla^2 \mvu.$$
where $K$ [$\mbox{m}^2$] is the permeability of the REV,
where $C_E$ [$\mbox{m}^{-1}$] is the Ergun coefficient that accounts for inertial (kinetic) effects
and where $\bar{\mu}$ is the effective viscosity.
The Ergun coefficient is strongly dependent on the flow regime.
For slow flows, $C_E$ is very small and the quadratic terme can be neglected.
The effective viscosity can be computed using Einstein’s formula \cite{einstein} for the effective viscosity
of a dilute suspension of rigid spherical particles in an ambient fluid:
$$\bar{\mu} \approx \mu \left(1+{5\over 2} (1-\phi)\right).$$
The permeability is interpreted in terms of spatial parameters $\phi$ and $d_p$
as
$$K = {d^2_p \phi^3 \over 150 (1-\phi)^2}.$$
Let us first neglect inertial effects ($C_E = 0$). Assume a domain of charachetistic length $L \gg d_p$. The following
adimensional number (called the Fidel Castro Number $\mbox{Fc}$) accounts for the relative importance of the large scale
viscous effects (Stokes term) and the small scale effects (Darcy term)
$$
\mbox{Fc} = {\mbox{Stokes} \over \mbox{Darcy}} = {K \over L^2} = {d^2_p \phi^3\over 150 L^2 (1-\phi)^2}.
$$
If a significative amount of particles are present in the REV ($\phi$ is not very close to $1$) , $\mbox{Fc} \ll 1$ and large scale effects viscous effects can be neglected. Yet, for
$\phi \simeq 1$, stokes effects can become proeminant and this is clearly the case when no particles are present in the REV. More precisely,
if $(1-\phi) L\gg d_p$, then large scale viscous effects are negligible.
\subsection*{New}
NB : all those papers deal with gaz-particle mixture in general and fluidized bed in particular. I'm not sure this can be applied to subsurface water.
While other forces can be seen in the literature \citep{zhu_discrete_2007}, the principal forces present in the gas-particle interactions are the drag force $G$ and the pressure gradient force $F$.
\subsubsection*{Pressure gradient force}
......
\section*{TODO}
\subsection*{For the Hydro-Quebec report}
\subsection*{1 Hydro-Quebec report {\small (Hopefully before July 3)}}
\begin{itemize}
\item coupling with lmgc
\item doc python
\item build with python instead of (or in addition of*) cmake
\item coupling with lmgc for hydroc-quebec 2d and 3d test cases
\item zone of prescribed porosity
\item doc python of hydroc-quebec 2d and 3d test cases
\item build with python instead of (or in addition of) cmake
\item doc scontactplot ?
\end{itemize}
\subsection*{For the model}
\subsection*{2 Model}
\begin{enumerate}
\item imex : include d/dt in the fluid equations, increase the order, find why there is a delay when we increase dt
\item friction in scontact
\item free surface/ale
\end{enumerate}
\subsection*{For the paper}
\subsection*{3 Paper}
\begin{itemize}
\item introduction
\item notation consistency
\item introduction + state of the art (partially done in my FNRS proposal)
\item conclusion
\item result section
\item result section (which testcases ?)
\item references (some can be found in my FNRS proposal)
\item polish/structure
\end{itemize}
\section{Validation}
\begin{figure}
\center
\includegraphics[width=0.7\textwidth]{fig/darcy/darcy_aligned_mesh.png}\\
\includegraphics[width=0.7\textwidth]{fig/darcy/darcy_aligned_pressure.png}\\
\includegraphics[width=0.7\textwidth]{fig/darcy/darcy_aligned_velocity.png}
\caption{aligned, $c_β = 0.4$, $d=0.0128$}
\end{figure}
\begin{figure}
\center
\includegraphics[width=0.7\textwidth]{fig/darcy/darcy_random_mesh.png}\\
\includegraphics[width=0.7\textwidth]{fig/darcy/darcy_random_pressure.png}\\
\includegraphics[width=0.7\textwidth]{fig/darcy/darcy_random_velocity.png}
\caption{random, $c_β = 0.4$, $d=0.0128$}
\end{figure}
\begin{figure}
\center
\begin{tabular}{cc}
\includegraphics[width=0.4\textwidth]{fig/darcy/pressure_aligned_01.pdf}&
\includegraphics[width=0.4\textwidth]{fig/darcy/pressure_aligned_02.pdf}\\
\includegraphics[width=0.4\textwidth]{fig/darcy/pressure_aligned_03.pdf}&
\includegraphics[width=0.4\textwidth]{fig/darcy/pressure_aligned_04.pdf}\\
\includegraphics[width=0.4\textwidth]{fig/darcy/pressure_aligned_05.pdf}&
\includegraphics[width=0.4\textwidth]{fig/darcy/pressure_aligned_06.pdf}\\
\includegraphics[width=0.4\textwidth]{fig/darcy/pressure_aligned_07.pdf}&
\includegraphics[width=0.4\textwidth]{fig/darcy/pressure_aligned_08.pdf}
\end{tabular}
\caption{Aligned particles}
\end{figure}
\begin{figure}
\center
\begin{tabular}{cc}
\includegraphics[width=0.4\textwidth]{fig/darcy/pressure_random_01.pdf}&
\includegraphics[width=0.4\textwidth]{fig/darcy/pressure_random_02.pdf}\\
\includegraphics[width=0.4\textwidth]{fig/darcy/pressure_random_03.pdf}&
\includegraphics[width=0.4\textwidth]{fig/darcy/pressure_random_04.pdf}\\
\includegraphics[width=0.4\textwidth]{fig/darcy/pressure_random_045.pdf}
\end{tabular}
\caption{Random particles}
\end{figure}
\begin{figure}
\center
\includegraphics[width=0.7\textwidth]{fig/darcy/K.pdf}\\
\includegraphics[width=0.7\textwidth]{fig/darcy/coeff.pdf}
\caption{Permeability}
\end{figure}
\section{Densité linéaire de particules}
$d'$ est une densité constante dans les coordonées $(x', y')$.
\begin{equation}
d(x, y) = \left(\frac{∂x}{∂x'} \frac{∂y}{∂y'} - \frac{∂x}{∂y'}\frac{∂y}{∂x'}\right)^{-1} d'
\end{equation}
une densité fonction de $x$ peut donc être obtenue en choisissant :
\begin{equation}
\frac{∂x}{∂x'} = \frac{∂y}{∂y'} = \left(\frac {d'}{d(x)}\right)^{\frac 1 2} \hspace*{1cm} \frac{∂x}{∂y'} \frac{∂y}{∂x'} = 0\hspace{1cm} \mbox{avec }L' = \frac{αL}{1-α}
\end{equation}
une densité constante en $y$ et linéaire (de 1 en 0 à $α$ en L) en $x$, s'obtient donc en posant :
\begin{equation}
\frac{∂x}{∂x'} = \frac{∂y}{∂y'} = \left(1 - \frac{x}{L'}\right)^{-\frac 1 2} \hspace*{1cm} \frac{∂x}{∂y'} \frac{∂y}{∂x'} = 0\hspace{1cm} \mbox{avec }L' = \frac{L}{1-α}
\end{equation}
d'où
\begin{align*}
x &= L - L\Big(1 ± \frac{3x'}{2L}\Big)^{\frac 2 3}
%x &= \Big[\Big(1 - \frac{3x'}{2L'}\Big)^\frac 2 3 - 1\Big] L' & y &= \Big(1 - \frac{3x'}{2L'}\Big)^{-\frac 1 3}y'
\end{align*}
\subsection{solution théorique}
Seulement Darcy :
\begin{align*}
∇p &= -\frac{μ}{K(x)}u \dx & K = \frac{_pφ(x)³}{150(1-φ(x))²}\\
Δp &= ∫_0^L -\frac{150(1-φ(x))²}{d_p²φ(x)³} μu\dx & φ(x)= φ_0 (1 + \frac x L (α - 1)) \\
&= -\frac{150 μuL}{d^2__0(α-1)}\left[-\frac{1}{2φ²} + \frac{2}{φ} + \log{φ}\right]_{φ_0}^{αφ_0}\\
&= -\frac{150 μuL}{d^2_p}\left(\frac{α+1}{2α²φ_} - \frac{2}{αφ_} + \frac{\log(α)}{(α - 1)φ_0}\right)
\end{align*}
Complete (avec $u = \text{cst}$) :
\begin{align*}
∇p &= \left(\frac{|∇φ|²}{φ³}-\frac{1}{K(x)}\right)μu\\
\end{align*}
Le rapport entre les deux est:
\[
\frac{1}{150}\left(\frac{|∇φ|d_p}{1-φ}\right)²
\]
Autrement dit, le terme supplémentaire est significatif quand $|∇φ| ≥ \frac{10}{d_p(1-φ)}$ ce qui est impossiple.
On a donc un terme supplémentaire :
\begin{align*}
_0^L\frac{φ_0² (α-1)²}{L²φ³}\dx
&=-\left[\frac{φ_0 (α-1)}{2Lφ^2}\right]_{φ_0}^{αφ_0}\\
&=-\frac{φ_0(α-1)(1-α^2)}{2Lα^_0^2}\\
&=\frac{(1-α)^2(1+α)}{2Lα^_0}
\end{align*}
%\subsection{flow around a porous cylinder}
%\begin{figure}
%\includegraphics[width=\textwidth]{fig/complete/porosity.png}\\
%\end{figure}
%\begin{figure}
%\includegraphics[width=\textwidth]{fig/complete/velocity.png}\\
%\includegraphics[width=\textwidth]{fig/complete/pressure.png}
%\caption {complete term with $-\tilde P^\alpha\nabla c_\alpha$}
%\end{figure}
%\begin{figure}
%\includegraphics[width=\textwidth]{fig/incomplete/velocity.png}\\
%\includegraphics[width=\textwidth]{fig/incomplete/pressure.png}
%\caption {incomplete term without $-\tilde P^\alpha\nabla c_\alpha$}
%\end{figure}
%\begin{figure}
%\includegraphics[width=\textwidth]{fig/orig/velocity.png}\\
%\includegraphics[width=\textwidth]{fig/orig/pressure.png}
%\caption {previous formulation}
%\end{figure}
\newpage
......@@ -4,3 +4,5 @@ from .dofManager import DofManager
from .meshJacobian import MeshJacobian
from .legendre import *
from .particleFluid import ParticleFluid
from .fluidSolver import FluidSolver
from .mesh import Mesh
import numpy as np
from .meshJacobian import MeshJacobian
from .particleFluid import ParticleFluid
from .Brinkman import Brinkman
from .dofManager import DofManager
from .Steady import Steady
from .mesh import Mesh
from .legendre import *
import time
import pyfefp
class Hydro :
def __init__(self, mesh, mesh_velocity_cb, dt, r, rhop = 1, rho = 1, mu = 1, g = [0, -1]) :
self._jac = pyfefp.MeshJacobian(mesh, 2)
self._pf = pyfefp.ParticleFluid(self._jac)
self._mesh_velocity = mesh_velocity_cb
class FluidSolver :
def __init__(self, mesh, dt, rhop, rho, mu, g) :
if isinstance(mesh, str):
mesh = Mesh(mesh)
self._jac = MeshJacobian(mesh)
self._pf = ParticleFluid(self._jac)
self._dt = dt
self.g = np.array(g)
self.rho = rho
#todo add vmesh in law
law = pyfefp.Brinkman(
law = Brinkman(
mu = mu,
pf = self._pf,
ns = False,
g = np.array([0, 0]),
rho = rho
)
d = pyfefp.DofManager(mesh)
d.add_field(pyfefp.TriangleP1Mini)
d.add_field(pyfefp.TriangleP1Mini)
d.add_field(pyfefp.TriangleP1)
d.add_boundary_condition("Pile", 0, cb = lambda x : mesh_velocity_cb(x)[:, 0])
d.add_boundary_condition("Pile", 1, cb = lambda x : mesh_velocity_cb(x)[:, 1])
d.add_boundary_condition("Box", 0, 0.)
d.add_boundary_condition("Box", 1, 0.)
d.add_boundary_condition("Top", 0, 0.)
d.add_boundary_condition("Top", 1, 0.)
d.add_boundary_condition("Top", 2, 0.)
d.complete()
self.rhop = rhop
g = np.array(g),
rho = rho,
rhop = rhop,
stab = dt)
d = DofManager(mesh)
if self._jac.dim == 2 :
d.add_field(TriangleP1Mini)
d.add_field(TriangleP1Mini)
d.add_field(TriangleP1)
elif self._jac.dim == 3:
d.add_field(TetrahedronP1Mini)
d.add_field(TetrahedronP1Mini)
d.add_field(TetrahedronP1Mini)
d.add_field(TetrahedronP1)
else :
raise ValueError("invalid mesh dimension : %i" % self._jac.dim)
self._dt = dt
self._law = law
self._dof = d
self._sol = d.new_vector()
self._solver = pyfefp.Steady(law, d, self._jac)
self._sol = None
def set_particles(self, pos, vol, vel, forces) :
self._pf.set_particles(pos, vol, vel, forces)
def add_boundary_condition(self, tag, field, *a, **b) :
self._dof.add_boundary_condition(tag, field, *a, **b)
def move_mesh(self) :
dx = self._mesh_velocity(self._jac.x.reshape([-1, 2])) * self._dt
self._jac.update(self._jac.x + dx.reshape([-1, 3, 2]))
def complete(self) :
self._dof.complete()
self._sol = self._dof.new_vector()
self._solver = Steady(self._law, self._dof, self._jac)
self._solverEx = Steady(self._law, self._dof, self._jac)
def solve(self) :
tic = time.clock()
v0 = self._pf.velocity.copy()
self._law.stab = self._dt *0.5
self._sol[:] = self._solver.solve()
self._pf.forces[:] += self._pf.volume * (self.rhop - self.rho) * self.g
self._pf.velocity += self._pf.forces * self._dt*0.5 / self._pf.mass
self._law.stab = 0
self._sol[:] = self._solverEx.solve()
self._pf.velocity[:] = v0
print("cpu : ", time.clock() - tic)
return self.forces
def write_position(self, odir, oiter, time, filetype="msh"):
if filetype == "msh" :
basename = "%s/f_%05i" % (odir, oiter)
X = np.zeros((9, self._jac.x.shape[0]))
X[[0,1,3,4,6,7], :] = self._jac.x.reshape([-1, 6]).T
self._jac.writeScalar(basename + "-xmesh.msh", "x", oiter, time, X)
else :
print("unknown export format : \"%s\"." % filetype)
def write_solution(self, odir, oiter, time, filetype="msh") :
print(oiter, time)
if filetype == "msh" :
basename = "%s/f_%05i" % (odir, oiter)
self._jac.writeScalar(basename + "-porosity.msh", "porosity", oiter, time, self._jac.dofManager.getField(0, self._pf.porosity))
self._jac.writeScalar(basename + "-pressure.msh", "pressure", oiter, time, self._dof.getField(2, self._sol))
X = np.zeros((9, self._jac.x.shape[0]))
X[[0,1,3,4,6,7], :] = self._jac.x.reshape([-1, 6]).T
self._jac.writeScalar(basename + "-xmesh.msh", "x", oiter, time, X)
X[0::3, :] = self._dof.getField(0, self._sol)[:-1,:]
X[1::3, :] = self._dof.getField(1, self._sol)[:-1,:]
self._jac.writeScalar(basename + "-velocity.msh", "velocity", oiter, time, X)
else :
print("unknown export format : \"%s\"." % filetype)
return
def mesh_position(self) :
return self._jac.x
def set_mesh_position(self, x) :
self._jac.update(x)
def set_particles(self, massp, volume, position, velocity) :
self.forces = np.zeros_like(velocity)
self._pf.set_particles(position, volume, velocity, self.forces, massp)
......@@ -71,7 +71,9 @@ class MeshJacobian :
return X
def __init__(self, mesh, dim) :
def __init__(self, mesh, dim = None) :
if not dim :
dim = max([e.dimension for e in mesh.entities])
d = dofManager.DofManager(mesh)
self.dofManager = d
self.dim = dim
......
......@@ -670,6 +670,22 @@ void particleProblemIterate(ParticleProblem *p, double alert, double dt, double
particleProblemSolve(p, alert, dt, tol, maxit);
for (size_t i = 0; i < vectorSize(p->position); ++i)
p->position[i] += p->velocity[i] * dt;
for (size_t i = 0; i < vectorSize(p->disks); ++i)
for (int j = 0; j < DIMENSION; ++j)
p->disks[i].x[j] += p->disks[i].v[j] * dt;
for (size_t i = 0; i < vectorSize(p->segments); ++i)
for (int j = 0; j < DIMENSION; ++j) {
p->segments[i].p[0][j] += p->segments[i].v[j] * dt;
p->segments[i].p[1][j] += p->segments[i].v[j] * dt;
}
#if DIMENSION == 3
for (size_t i = 0; i < vectorSize(p->triangles); ++i)
for (int j = 0; j < DIMENSION; ++j) {
p->triangles[i].p[0][j] += p->triangles[i].v[j] * dt;
p->triangles[i].p[1][j] += p->triangles[i].v[j] * dt;
p->triangles[i].p[2][j] += p->triangles[i].v[j] * dt;
}
#endif
}
size_t particleProblemAddBoundarySegment(ParticleProblem *p, const double x0[DIMENSION], const double x1[DIMENSION], int tag) {
......
......@@ -15,7 +15,7 @@ class Hydro :
mu = mu,
pf = self._pf,
ns = False,
g = self.g, #np.array([0, 0]),
g = self.g,
rho = rho,
rhop = rhop,
stab = dt)
......
from pyfefp.mesh import Mesh
from hydro import Hydro
import scontact2
import pyfefp
import numpy as np
import shutil
import time
import os
import sys
import lmgcInterface
import scontactInterface
class MeshVelocity :
def __init__(self) :
self.Y = 0
self.V = 0
self.position = 0
self.velocity = 0
def cb(self, x) :
def call(self, x) :
l = 0.1
L = 2
h = 0.5
H = 1.5
fx = np.minimum(1 - (np.abs(x[:, 0]) - l) / (L - l), 1.)
HX = self.Y * fx
HH = np.where(x[:, 1] > HX, H - HX , H + HX)
fy = np.minimum(1 - (np.abs(x[:, 1] - HX) - h) / (HH - h), 1.)
fx = np.minimum(1 - (np.abs(x[..., 0]) - l) / (L - l), 1.)
HX = self.position * fx
HH = np.where(x[..., 1] > HX, H - HX , H + HX)
fy = np.minimum(1 - (np.abs(x[..., 1] - HX) - h) / (HH - h), 1.)
v = np.zeros_like(x)
v[:, 0] = 0
v[:, 1] = fx * fy * self.V
v[..., 0] = 0
v[..., 1] = fx * fy * self.velocity
return v
def update(self, dt) :
self.Y += self.V * dt
self.position += self.velocity * dt
def writep(p, odir, i) :
outputname = "%s/part-%05i" % (odir, i)
tmpoutputname = "%s/part-%05i_tmp" % (odir, i)
p.write(tmpoutputname)
shutil.move(tmpoutputname, outputname)