Commit 695429b5 authored by Michel Henry's avatar Michel Henry
Browse files

wip

parent 9000a7fa
......@@ -419,8 +419,8 @@ static void particle_force_f(FluidProblem *problem, double *f, double *dfdu, dou
static void fluid_problem_f(const FluidProblem *problem, const double *sol, double *dsol, const double *solold, const double *dsolold, const double* mesh_velocity, const double c, const double *dc, const double *bf, const double cold, const double rho, const double mu, double dt, int iel, double *f0, double *f1, double *f00, double *f10, double *f01, double *f11) {
double p = sol[P];
double taup = problem->taup[iel];
double tauc = problem->tauc[iel];
double taup = problem->taup[iel]/16.0;
double tauc = problem->tauc[iel]/16.0;
const int n_fields = fluid_problem_n_fields(problem);
double u[D], uold[D], du[D][D], duold[D][D], umesh[D], dp[D];
if(problem->stab_param != 0)
......
......@@ -147,13 +147,40 @@ def compute_drag():
# COMPUTATION LOOP
#
while t < tEnd :
# fluid.move_particles(p.position(), p.velocity(), p.contact_forces())
# fluid.implicit_euler(dt)
# fext = g*p.mass() + compute_drag()#fluid.compute_node_force(dt)
# p.iterate(dt, fext)
# time_integration.iterate(fluid, p, dt, min_nsub=5, external_particles_forces=g*p.mass())
fluid.move_particles(p.position(), p.velocity(), p.contact_forces())
fluid.implicit_euler(dt)
fext = g*p.mass() + compute_drag()#fluid.compute_node_force(dt)
p.iterate(dt, fext)
# time_integration.iterate(fluid, p, dt, min_nsub=5, external_particles_forces=g*p.mass())
t += dt
time_integration._advance_particles(p, fext, dt, 5, 1e-8)
# predicteur correcteur
# alpha = 0.5
# sf0 = np.copy(fluid.solution())
# p.save_state()
# predicteur
# fluid.implicit_euler(dt)
# sf1 = np.copy(fluid.solution())
# f0 = g*p.mass() + compute_drag()
# time_integration._advance_particles(p,f0,dt,min_nsub=5,contact_tol=1e-8)
# fluid.particle_velocity()[:] = p.velocity()
# correcteur
# fluid.solution()[:,:] = sf0
# fluid.implicit_euler(dt)
# fluid.solution()[:,:] = (alpha*fluid.solution()+(1-alpha)*sf1)
# f1 = g*p.mass() + compute_drag()
# p.restore_state()
# time_integration._advance_particles(p,f1,dt,min_nsub=5,contact_tol=1e-8)
t += dt
# Output files writting
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
......
......@@ -33,8 +33,6 @@ from migflow import time_integration
import numpy as np
import os
import time
import shutil
import random
np.random.seed(0)
def genInitialPosition(filename, r, rout, rhop, compacity) :
......@@ -68,24 +66,24 @@ def genInitialPosition(filename, r, rout, rhop, compacity) :
p.write_vtk(filename,0,0)
# Define output directory
outputdir = "output"
outputdir = "output_adapt_tau_by_16"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
# Physical parameters
g = np.array([0,-9.81]) # gravity
rhop = 2450 # particles density
r = 154e-6/2.5 # particles radii
compacity = 0.20 # solid volume fraction in the drop
r = 25e-6#154e-6/2.5 # particles radii
compacity = 0.05#0.20 # solid volume fraction in the drop
rho = 1030 # fluid density
nu = 1.17/rho # kinematic viscosity
rout = 3.3e-3 # cloud radius
mu = nu*rho # dynamic viscosity
print("RHOP = %g, r = %g, RHO = %g" % (rhop,r,rho))
print("RHOP = %g, r = %g, RHO = %g, MU = %g" % (rhop,r,rho, mu))
# Numerical parameters
outf = 20 # number of iterations between output files
dt = 2e-3 # time step
outf = 10 # number of iterations between output files
dt = 2.5e-3 # time step
tEnd = 100 # final time
#
......@@ -115,9 +113,9 @@ fluid.set_wall_boundary("Lateral")
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.contact_forces())
t = 0
fluid.solution()[:,2] = (fluid.coordinates()[:,1]-0.6)*rho*g[1]
fluid.write_vtk(outputdir,0,0)
tic = time.time()
fluid.solution()[:,2] = fluid.coordinates()[:,1]*rho*g[1]
def write_points(p):
with open("adapt/p.geo","w") as f:
......@@ -144,8 +142,7 @@ exit(0)
"""
def compute_drag():
pv, pdp, pc = fluid.solution_at_particles()
dv = p.velocity()-pv
dv = (p.velocity()-pv)/pc
d = 2*(p.volume()/np.pi)**0.5
normu = np.hypot(dv[0],dv[1])
#Reynolds/|u_p-u_f|
......@@ -154,6 +151,7 @@ def compute_drag():
Cd_u = Cd_u*Cd_u
f = pc**(-(1.8-0.65*np.exp(-(1.5-np.log(Re_O_u*normu))*(1.5-np.log(Re_O_u*normu))/2.)))
gamma = f*Cd_u*rho/2*d
# v is dv
# m(v'-v)/dt = mg - V*dp - gamma * v'
# v' (m/dt + gamma) = mg -V*dp + vm/dt
......@@ -171,26 +169,21 @@ def compute_drag():
print(np.mean(np.linalg.norm(gammastar*dvstar, axis=1)), np.mean(np.linalg.norm(gamma*dv,axis=1)))
return -gammastar*dvstar-p.volume()*pdp
while t < tEnd :
if ((ii+1)%50==0):
if ((ii+1)%200==0):
number_p = fluid.n_particles
position_p = fluid.particle_position()
volume_p = fluid.particle_volume()
if (ii%50==0 and ii != 0):
if (ii%200==0 and ii != 0):
write_points(p)
fluid.adapt_mesh(2e-2,1e-3,10000,old_n_particles=number_p,old_particle_position=position_p,old_particle_volume=volume_p)
#time_integration.iterate(fluid, p, dt, external_particles_forces=g*p.mass())
fluid.move_particles(p.position(), p.velocity(), p.contact_forces())
fluid.implicit_euler(dt)
fext = g*p.mass() + compute_drag()#fluid.compute_node_force(dt)
time_integration._advance_particles(p, fext, dt, 5, 1e-8)
# p.iterate(dt, fext)
#p.velocity()[:] += fext/p.mass()*dt
#p.position()[:] += p.velocity()*dt
#time_integration._advance_particles(p, fext, dt, 5, 1e-8)
t += dt
# Output files writting
......
This diff is collapsed.
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