Hi everyone!
I have been trying to implement a nonlinear viscosity (shear thickening) in a Navier-Stokes-like problem. for the reference this is the 2D equation I’m trying to solve:
\omega v_i =
\eta(\dot{\gamma})\, \nabla^2 v_i
+ \varepsilon_1\, \partial_i \!\left( \frac{\rho}{\rho + \rho_0} \right)
\\
\partial_t \rho
+ \nabla \cdot (\rho \mathbf{v})
=
- k \, (\rho - \rho_0)
+ D \, \nabla^2 \rho
v is the velocity, rho is density field, omega is friction and viscosity is dependent on the shear rate as follows:
\eta(\dot{\gamma})
= \dfrac{\eta_0 + \alpha \eta_1 \dot{\gamma}^2}
{1 + \alpha \dot{\gamma}^2}
I’m not sure whether I’m implementing the nonlinearity correctly. I have gone through the nonlinear examples in the documentation, but I found them difficult to follow since they use a different approach from mine. Here is my implementation so far:
I have define the equations in a class and most of the code is commented so it should be easy to follow, the nonlinear shear viscosity is in def setup nonlinear_form. Two points I’m not sure about:
-
how should I define the shear rate in my code? Currently, I use the velocity from previous time step just by calling gpu.components. is this the correct way?
-
I find that I need to add a constant viscosity term to the bilinear form in order to stabilize the equation, even though this term does not appear explicitly in my model equation. is it because the way I’m implementing the nonlinear viscosity?
Thanks for your help!
class Navierstokes:
def __init__(self,
width=1,
height=1,
maxh=0.05,
gamma=0.5,
eta_1=0.25,
k=0.01,
D=0.01,
chi0=0.15,
rho0=1,
alpha = 0.1,
eta0 = 1,
eta2 = 0.1,
):
self.width = width
self.height = height
self.maxh = maxh
# Physical parameters
self.gamma = gamma
self.eta_1 = eta_1
self.eta_2 = eta_2
self.k = k
self.D = D
self.chi0 = chi0
self.rho0 = rho0
self.alpha = alpha
self.eta0 = eta0
self.eta2 = eta2
# Initialize simulation components
self._create_mesh()
self._setup_function_spaces()
self._setup_initial_conditions()
def _create_mesh(self):
"""Create the rectangular mesh with named boundaries."""
shape = Rectangle(self.width, self.height).Face()
shape.edges.Min(X).name = "right"
shape.edges.Max(X).name = "left"
shape.edges.Min(Y).name = "up"
shape.edges.Max(Y).name = "down"
self.mesh = Mesh(OCCGeometry(shape, dim=2).GenerateMesh(maxh=self.maxh)).Curve(3)
def _setup_function_spaces(self):
self.V = VectorH1(self.mesh, order=3, dirichlet="right|left|up|down") # velocity field
self.R = H1(self.mesh, order=2) # density field
self.X = self.V * self.R
self.gfu = GridFunction(self.X)
self.velocity, self.density= self.gfu.components
self.time = Parameter(0)
def _setup_initial_conditions(self):
"""Set up initial conditions."""
self.density.Set(1)
v = CoefficientFunction((0, 0))
self.velocity.Set(v)
self.time.Set(0)
def _shear_dependent_viscosity(self, v):
gamma_dot = self._compute_shear_rate(v)
eta_apparent = (self.eta0 + self.eta2 * self.alpha * gamma_dot**2)/ (1+ self.alpha* gamma_dot**2)
return eta_apparent
def _compute_shear_rate(self, v):
D11 = grad(v)[0, 0]
D22 = grad(v)[1, 1]
D12 = 0.5 * (grad(v)[0, 1] + grad(v)[1, 0])
shear_rate_sq = 2 * (D11*D11 + D22*D22 + 2*D12*D12)
epsilon = 1e-10
shear_rate = sqrt(shear_rate_sq + epsilon)
return shear_rate
def _setup_bilinear_form(self, functions):
(v_trial, rho_trial), \
(v_test,rho_test) = functions
# velocity equation
force_bal = (self.gamma * InnerProduct(v_trial, v_test) * dx +
self.eta_1 * InnerProduct(grad(v_trial), grad(v_test)) * dx
)
# Density evolution equation
lhs_rho_time = (self.k * rho_trial * rho_test * dx +
self.D * InnerProduct(grad(rho_trial), grad(rho_test)) * dx)
return force_bal + lhs_rho_time
def _setup_nonlinear_form(self, functions):
(v_trial,rho_trial),\
(v_test,rho_test) = functions
# Advection term for density
advection_rho = (InnerProduct(grad(rho_trial), v_trial) * rho_test * dx +
rho_trial * div(v_trial) * rho_test * dx)
# Nonlinear term in velocity equation
force_bal = (self.chi0(self.time) * rho_trial * div(v_test) * 2 / (rho_trial + 1) * dx
)
# nonlinear shear viscosity
eta_eff = self._shear_dependent_viscosity(self.velocity)
viscosity_nonlinear = eta_eff * InnerProduct(grad(v_trial), grad(v_test)) * dx
return advection_rho + force_bal + viscosity_nonlinear
def _setup_linear_form(self, functions):
"""Set up linear form."""
_, (_, rho_test, _, _) = functions
return self.k(self.time) * self.rho0(self.time) * rho_test * dx
def _setup_inverse_form(self, functions, bilinear, tau):
"""Set up inverse. """
(_, rho_trial), \
(_, rho_test) = functions
inverse_form = (rho_trial * rho_test * dx +
tau * bilinear)
return inverse_form
def _setup_forms(self, tau):
"""Set up the bilinear and linear forms for the simulation."""
functions = self.X.TnT()
# Set up bilinear forms
bilinear = self._setup_bilinear_form(functions)
self.a = BilinearForm(self.X)
self.a += bilinear
self.a.Assemble()
# Nonlinear terms
self.nonlinear = BilinearForm(self.X, nonassemble=True)
self.nonlinear += self._setup_nonlinear_form(functions)
# Linear terms
self.f = LinearForm(self.X)
self.f += self._setup_linear_form(functions)
self.f.Assemble()
# Inverse form
mstar = BilinearForm(self.X)
mstar += self._setup_inverse_form(functions, bilinear, tau)
mstar.Assemble()
self.inv = mstar.mat.Inverse(freedofs=self.X.FreeDofs(), inverse="sparsecholesky")
def set_initial_density(self, density_function):
"""Set custom initial density distribution."""
self.density.Set(density_function)
def set_initial_velocity(self, velocity_function):
"""Set custom initial velocity field."""
self.velocity.Set(velocity_function)
def simulate(self, tend=10, tau=0.01, save_interval=1):
"""Run the simulation."""
t = 0
self._setup_forms(tau)
# Create multidimensional GridFunction for animation
self.gfut = GridFunction(self.gfu.space, multidim=0)
with TaskManager():
for i in tqdm(range(int(tend/tau))):
self.time.Set(t)
self.a.Assemble()
self.nonlinear.Assemble()
res = (self.a.mat * self.gfu.vec +
self.nonlinear.mat * self.gfu.vec -
self.f.vec)
self.gfu.vec.data -= tau * self.inv * res
if i % save_interval == 0:
self.gfut.AddMultiDimComponent(self.gfu.vec)
t += tau
def visualize(self, animate=True):
"""
Visualize the simulation results.
Parameters:
animate (bool): If True, show animation of results
"""
Draw(self.velocity, self.mesh, "velocity", vectors=True)
Draw(self.density, self.mesh, "density")