Commit d2bcead3 authored by Michel Henry's avatar Michel Henry
Browse files

wip

parent 5809f12e
Pipeline #9738 passed with stages
in 3 minutes and 54 seconds
......@@ -15,212 +15,7 @@ import numpy as np
import os
import time
# Add polymers if available space
def genMorePolymers(fluid, p, joints, polymers, r_max, rhop, inter_dist, max_dist, jab, n, H):
I = 15
border = 0.9
x = p.position()
index = len(x)
for i in range(I):
if len(x[(x[:,0] < border) & (x[:,1] >= (i+1)*H/(I+2)) & (x[:,1] < (i+2)*H/(I+2))]) == 0:
p1 = np.array([border-2*r_max, (i+1)*H/(I+2)])
p.add_particle(p1, r_max, r_max**2*np.pi*rhop, "Glass")
index += 1
for y in range(1,n):
sp = p1-np.array([inter_dist*y,0])
p.add_particle(sp, r_max, r_max**2*np.pi*rhop, "Glass")
joints.append([index-1, index, jab, max_dist])
index += 1
polymers += [list(range(index-n, index))]
if index != len(x):
print("Add polymers")
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(), p.contact_forces())
# Remove polymers at the end of the domain
def removeUselessParticls(fluid, p, joints, polymers, L, n):
x = p.position()
v = p.velocity()
element_id = fluid.particle_element_id()[:]
to_remove = 0
flags = np.full(len(x), True, dtype=bool)
for polymer in polymers:
polymer = np.array(polymer)
if len(polymer[element_id[polymer] == -1]) > n/5 or any(np.abs(v[polymer,0]) > 100) or any(np.abs(v[polymer,1]) > 100):
to_remove += 1
flags[polymer] = False
if to_remove > 0:
p.remove_particles_flag(flags)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(), p.contact_forces())
polymers = polymers[:-to_remove]
joints = joints[:-to_remove*(n-1)]
return joints, polymers
# Function used to apply surface tension at the free surface
def apply_tension_weak(x):
fs_nodes = fs.fs_nodes
lapl_h = fs.compute_lapl_fs(fs_nodes)
n_e = len(lapl_h)-1
kappa = np.zeros_like(x[:,0])
X = fluid.coordinates()[fs_nodes,0]
h = fluid.coordinates()[fs_nodes,1]
dX = X[1:] - X[:-1]
dhdx = (h[1:] - h[:-1])/dX[:]
dl = np.sqrt(dhdx[:]*dhdx[:] + 1)**(-3)
ind = np.argsort(x[:,0])
for i in range(0,n_e):
ind_i = ind[2*i:2*i+2]
xi = (2*x[ind_i,0] - (X[i+1] + X[i]))/(X[i+1]- X[i])
phi = np.column_stack([(1-xi)/2, (1+xi)/2])
kappa[ind_i] = (phi[0,:]*lapl_h[i] + phi[1,:]*lapl_h[i+1])*dl[i]
tension = -gamma*kappa[:]
# print(f"\tSurface tension : {max(abs(tension))}")
return tension
# Define output directory
outputdir = "output_6"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
# Physical parameters
g = 0 # gravity
rho = 1 # fluid density
nu = 1 # kinematic viscosity
mu = 0.8 # dynamic viscosity
gamma = 0 # superficial tension coefficient
v = 8 # shearing velocity
# Particles
r_max = 6e-4 # max particles size
n_chain = 100 # number of particles/polymer
start_free_surface = 1 # die lip
noise_force = 4e-2 # Brownian motion
noise_variance = 25*np.pi/180
mbound = n_chain*(4*r_max)
Hi = 2e-3 # FENE
Li = 40*r_max
# Numerical parameters
outf = 5 # number of iterations between output files
dt = 1e-4 # time step
tEnd = 8 # final time
# Geometrical parameters
H = 0.2 # domain height
L = 3 # domain width
Ox, Oy = 0, 0 # leftmost and bottommost point
# Initial time and iteration
t = 0
ii = 0
#
# PARTICLES PROBLEM
#
p = scontact.ParticleProblem(2,True)
p.load_msh_boundaries("mesh.msh", ["TopLeft", "Left", "BottomLeft", "BottomRight"],material="PVC") ## On peut ne pas importer
joints = []
polymers = []
p.set_friction_coefficient(.3,"Glass","Glass") # between particles
p.set_friction_coefficient(0,"Glass","PVC") # between particles and boundaries
G = np.zeros((p.n_particles(),2))
G[:,1] = p.mass()[:,0]*-(20)
#
# FLUID PROBLEM
#
fluid = mbfluid.FluidProblem(2,[0,g],[mu],[rho])
# Set the mesh geometry for the fluid computation
fluid.load_msh("mesh.msh")
fluid.set_open_boundary("Left", velocity = [lambda x : v*(1-(x[:,1]/H)**2),0])
fluid.set_strong_boundary("Left",0, lambda x : v*(1-(x[:,1]/H)**2))
fluid.set_strong_boundary("Left",1,0)
fluid.set_wall_boundary("TopLeft", velocity = [0,0])
fluid.set_open_boundary("TopRight", pressure = apply_tension_weak, compute_viscous_term = 0)
fluid.set_symmetry_boundary("BottomLeft")
fluid.set_symmetry_boundary("BottomRight")
fluid.set_open_boundary("Right", pressure = 0)
#
# FREE SURFACE STRUCTURES
#
fs = time_integration.FreeSurface(fluid, "TopRight", "BottomRight", isCentered=False)
fs_nodes = np.array(fs.fs_nodes)
ending_nodes = fs_nodes[np.argsort(fluid.coordinates()[fs_nodes, 0])][-16:]
fluid.write_vtk(outputdir,0,0)
p.write_vtk(outputdir, 0, 0)
random_angles = (np.random.rand(n_chain, 2)*2-.5)*2*r_max # Perturbation
tic = time.time()
#
# COMPUTATION LOOP
#
while t < tEnd :
if ii >= 10:
genMorePolymers(fluid, p, joints, polymers, r_max, 2.0, 3*r_max, 40*r_max, 2e-3, n_chain, 0.2)
joints, polymers = removeUselessParticls(fluid, p, joints, polymers, L, n_chain)
Gp = np.zeros((p.n_particles(),2))
xs = p.position()
vs = p.velocity()
t4 = time.time()
for polymer in polymers:
gravity_center = np.mean(xs[polymer], axis=0)
if ii %30 == 0:
random_angles = (np.random.rand(n_chain, 2)*2-1)*5*r_max # Perturbation
vect_to_gc = gravity_center - xs[polymer]
vect_to_gc += random_angles
vect_to_gc[np.linalg.norm(vect_to_gc, axis=1) < r_max/10,:] = 0
if gravity_center[0] < start_free_surface: # To avoid too many polymers
vect_to_gc *= 0
Gp[polymer,:] += vect_to_gc*noise_force/mbound
t5 = time.time()
p.velocity()[:] = vs
joints_force1 = np.zeros((p.n_particles(),2))
joints_force2 = np.zeros((p.n_particles(),2))
translation = xs[1:,:] - xs[:-1,:]
joints_force1[1:,0] -= translation[:,0]*Hi/(1-np.power(np.linalg.norm(translation, axis=1)/Li, 2))
joints_force1[1:,1] -= translation[:,1]*Hi/(1-np.power(np.linalg.norm(translation, axis=1)/Li, 2))
joints_force2[:-1,0] += translation[:,0]*Hi/(1-np.power(np.linalg.norm(translation, axis=1)/Li, 2))
joints_force2[:-1,1] += translation[:,1]*Hi/(1-np.power(np.linalg.norm(translation, axis=1)/Li, 2))
joints_force1[::n_chain,:] = 0
joints_force2[n_chain-1::n_chain,:] = 0
Gp += joints_force1 + joints_force2
t6 = time.time()
p.velocity()[xs[:,0] < start_free_surface,1] = 0
Gp[xs[:,0] < start_free_surface] *= 0
t8 = time.time()
time_integration.iterate(fluid, p, dt, min_nsub=5, external_particles_forces=Gp)
t9 = time.time()
fs.iterate(dt)
fluid.move_particles(p.position(), p.velocity(), p.contact_forces())
t10 = time.time()
print("Noise : " + str(t5-t4))
print("Joints : " + str(t6-t5))
print("Iteration : " + str(t9-t8))
print("Surface libre : " + str(t10-t9))
else:
time_integration.iterate(fluid, None, dt*5)
fs.iterate(dt)
fluid.move_particles(p.position(), p.velocity(), p.contact_forces())
t += dt
# Output files writting
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
fluid.export_vtk(outputdir, t, ioutput)
p.write_vtk(outputdir, ioutput, t)
ii += 1
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))
"""
# Add polymers if available space
def genMorePolymers(fluid, p, joints, polymers, r_max, rhop, inter_dist, max_dist, jab, n, H):
I = 15
......@@ -441,5 +236,4 @@ while t < tEnd :
fluid.write_vtk(outputdir, ioutput, t)
p.write_vtk(outputdir, ioutput, t)
ii += 1
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))
"""
\ No newline at end of file
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))
\ No newline at end of file
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