Commit 82f5a7e4 authored by Matthieu Constant's avatar Matthieu Constant
Browse files

predictor-corrector for drag in stab

parent b7deb87c
Pipeline #5386 passed with stage
in 27 seconds
# MigFlow - Copyright (C) <2010-2018>
# <Universite catholique de Louvain (UCL), Belgium
# Universite de Montpellier, France>
#
# List of the contributors to the development of MigFlow: see AUTHORS file.
# Description and complete License: see LICENSE file.
#
# This program (MigFlow) is free software:
# you can redistribute it and/or modify it under the terms of the GNU Lesser General
# Public License as published by the Free Software Foundation, either version
# 3 of the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with this program (see COPYING and COPYING.LESSER files). If not,
# see <http://www.gnu.org/licenses/>.
#!/usr/bin/env python
# TESTCASE DESCRIPTION
# Bidimensional particles sedimentation in fluid
from migflow import fluid
from migflow import scontact
import numpy as np
import os
import time
import shutil
import random
def genInitialPosition(filename, r, H, ly, lx, rhop) :
"""Set all the particles centre positions and create the particles objects to add in the computing structure
Keyword arguments:
filename -- name of the output file
r -- max radius of the particles
H -- domain height
ly - particles area height
lx -- particles area width
rhop -- particles density
"""
# Particles structure builder
p = scontact.ParticleProblem(2)
# Load mesh.msh file specifying physical boundaries names
p.load_msh_boundaries("mesh.msh", ["Top", "Lateral","Bottom"])
#Definition of the points where the particles are located
x = np.arange(-lx/2+r, lx/2-r, 2*r)
y = np.arange(H/2-r, H/2-ly+r, -2*r)
x, y = np.meshgrid(x, y)
x = x.flat
y = y.flat
# Add a grain at each centre position
for i in range(len(x)) :
p.add_particle((x[i], y[i]), r, r**2 * np.pi * rhop)
p.write_vtk(filename,0,0)
# Define output directory
outputdir = "output"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
# Physical parameters
g = -9.81 # gravity
r = 1e-3 # particles radius
rhop = 1500 # particles density
rho = 1000 # fluid density
nu = 1e-6 # kinematic viscosity
# Numerical parameters
outf = 1 # number of iterations between output files
dt = 0.25e-2 # time step
tEnd = 100 # final time
# Geometrical parameters
ly = 5e-2 # particles area height
lx = 4e-1 # particles area widht
H = 1.2 # domain height
#
# PARTICLE PROBLEM
#
# Initialise particles
genInitialPosition(outputdir, r, H, ly, lx, rhop)
p = scontact.ParticleProblem(2)
p.read_vtk(outputdir,0)
print("r = %g, m = %g\n" % (p.r()[0], p.mass()[0]))
print("RHOP = %g" % rhop)
# Initial time and iteration
t = 0
ii = 0
#
# FLUID PROBLEM
#
fluid = fluid.FluidProblem(2,g,[nu*rho],[rho],drag_in_stab=1)
# Set the mesh geometry for the fluid computation
fluid.load_msh("mesh.msh")
fluid.set_wall_boundary("Bottom")
fluid.set_wall_boundary("Lateral")
fluid.set_wall_boundary("Top",pressure=0)
fluid.set_strong_boundary("Top",0,0)
fluid.set_strong_boundary("Top",1,0)
fluid.set_strong_boundary("Top",2,0)
fluid.set_strong_boundary("Lateral",0,0)
fluid.set_strong_boundary("Lateral",1,0)
fluid.set_strong_boundary("Bottom",0,0)
fluid.set_strong_boundary("Bottom",1,0)
fluid.solution()[:,2] = (fluid.coordinates()[:,1]-0.6)*rho*g
# Set location of the particles in the mesh and compute the porosity in each computation cell
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.contact_forces())
fluid.export_vtk(outputdir,0,0)
tic = time.time()
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.contact_forces())
#
# COMPUTATION LOOP
#
while t < tEnd :
s = np.copy(fluid.solution())
vp = np.copy(p.velocity())
### predictor
# Fluid solver
fluid.implicit_euler(dt)
# Computation of the new velocities
forces = fluid.compute_node_force(dt)
p.velocity()[:,:] = p.velocity() + forces*dt/p.mass()
### corrector
fluid.solution()[:,:] = s[:,:]
# Fluid solver
fluid.implicit_euler(dt)
# Computation of the new velocities
forces = fluid.compute_node_force(dt)
p.velocity()[:,:] = vp[:,:]
vn = p.velocity() + forces*dt/p.mass()
vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
# Number of contact sub time step
nsub = max(1, int(np.ceil((vmax*dt*4)/min(p.r()))))
print("NSUB", nsub,"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(p.r()))
# NLGS iterations
for i in range(nsub) :
tol = 1e-8 #numerical tolerance for particles intersection
p.iterate(dt/nsub, forces, tol)
# Localisation of the particles in the fluid
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.contact_forces())
t += dt
# Output files writting
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
p.write_vtk(outputdir, ioutput, t)
fluid.export_vtk(outputdir, t, ioutput)
ii += 1
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment