Commit ebe2ec0e authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

surface tension wip

parent b04303b0
Pipeline #4794 passed with stage
in 26 seconds
......@@ -748,6 +748,81 @@ static void fluid_problem_volume(FluidProblem *problem, const double *solution_o
free(h1);
}
static void fluid_problem_surface_tension(FluidProblem *problem, double *all_local_vector){
if (problem->n_fluids == 1) return;
const Mesh *mesh = problem->mesh;
for (int iel = 0; iel < mesh->n_elements; ++iel) {
const int *el = &mesh->elements[iel*N_N];
double dxdxi[D][D],dxidx[D][D];
for (int i=0; i<D; ++i)
for (int j=0; j<D; ++j)
dxdxi[i][j] = mesh->x[el[j+1]*3+i]-mesh->x[el[0]*3+i];
invDD(dxdxi, dxidx);
}
double sigma = 22.27*1e-9;
int n_fields = fluid_problem_n_fields(problem);
double *a_cg = malloc(sizeof(double)*mesh->n_nodes);
for (int i = 0; i< mesh->n_nodes; ++i) {
a_cg[i] = 0;
}
for (int iel = 0; iel < mesh->n_elements; ++iel) {
const int *el = &mesh->elements[iel*N_N];
double dxdxi[D][D],dxidx[D][D];
for (int i=0; i<D; ++i)
for (int j=0; j<D; ++j)
dxdxi[i][j] = mesh->x[el[j+1]*3+i]-mesh->x[el[0]*3+i];
double det = invDD(dxdxi, dxidx);
#if DIMENSION == 2
double vol = det/6;
#else
double vol = det/24;
#endif
for (int i = 0; i< N_N; ++i) {
a_cg[el[i]] += vol * problem->concentration[iel*N_N+i];
}
}
for (int i = 0; i< mesh->n_nodes; ++i) {
a_cg[i] /= problem->node_volume[i];
}
for (int iel = 0; iel < mesh->n_elements; ++iel) {
const int *el = &mesh->elements[iel*N_N];
double dxdxi[D][D],dxidx[D][D];
for (int i=0; i<D; ++i)
for (int j=0; j<D; ++j)
dxdxi[i][j] = mesh->x[el[j+1]*3+i]-mesh->x[el[0]*3+i];
invDD(dxdxi, dxidx);
double dphi[N_N][D];
grad_shape_functions(dxidx, dphi);
double a = 0;
double da[D] = {0};
double nda = 0;
for (int i = 0; i < N_SF; ++i) {
a += a_cg[el[i]];
}
a /= N_N;
//double fa = a>0.5 ? -1 : 1;
double fa = -1;
for (int k=0; k<D; ++k){
for (int i = 0; i < N_SF; ++i) {
da[k] += a_cg[el[i]]*dphi[i][k];
}
nda += da[k]*da[k];
}
nda = fmax(sqrt(nda),1e-8);
double *local_vector = &all_local_vector[iel*n_fields*N_N];
for (int iphi = 0; iphi < N_N; ++iphi) {
for (int id = 0; id < D; ++id) {
for (int jd = 0; jd < D; ++jd) {
local_vector[(U+id)*N_N+iphi] += da[id]*da[jd]*dphi[iphi][jd]/nda*sigma*fa;
}
local_vector[(U+id)*N_N+iphi] -= nda*dphi[iphi][id]*sigma*fa;
}
}
}
free(a_cg);
}
static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, const double *solution_old, double dt, int reduced_gravity, double stab_param)
{
HXTLinearSystem *lsys = problem->linear_system;
......@@ -785,6 +860,7 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
node_force_volume(problem, solution_old, dt, all_local_vector, all_local_matrix,reduced_gravity);
compute_weak_boundary_conditions(problem, dt, all_local_vector, all_local_matrix);
fluid_problem_volume(problem, solution_old, dt, all_local_vector, all_local_matrix,reduced_gravity,stab_param);
fluid_problem_surface_tension(problem, all_local_vector);
for (int iel=0; iel < mesh->n_elements; ++iel) {
double *local_vector = &all_local_vector[local_size*iel];
double *local_matrix = &all_local_matrix[local_size*local_size*iel];
......@@ -1273,6 +1349,7 @@ static void fluid_problem_set_mesh(FluidProblem *problem, Mesh *mesh) {
problem->boundaries = mesh_boundaries_create(problem->mesh);
}
void fluid_problem_adapt_mesh(FluidProblem *problem, double lcmax, double lcmin, double n_el, double cmax, double cmin)
{
struct timespec tic, toc;
......
# 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.
# 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:
......@@ -38,11 +37,11 @@ ii = 0
#physical parameters
g = -9.81
g = 0# -9.81
rhof = 1030 # fluid density
rhop = 2450 # grains density
r=154e-6 # grains radii
R = 3.3e-3 # drop radius
R = 0.01 # drop radius
compacity = .20 # solid volume fraction in the drop
phi = 1 - compacity
nuf = 1.17/rhof # kinematic viscosity
......@@ -57,8 +56,8 @@ print("rhom = %g, num = %g" % (rhom,num))
#numerical parameters
tEnd = 200 # final time
dt = 1e-3 # time step
outf = 10 # number of iterations between output files
dt = 1e-2 # time step
outf = 5 # number of iterations between output files
def outerBndV(x) :
print(.4*max(np.sin(t*np.pi*2./1),0))
......@@ -76,12 +75,16 @@ s = fluid.solution()
c = np.ndarray((fluid.n_nodes()))
x = fluid.coordinates()
for i in range(len(x[:,0])):
z = (x[i,0])**2+(x[i,1]-.52)**2
R1 = (0.9*R)**2
R2 = (1.1*R)**2
#s[i,1] = -0.01
c[i] = min(max(0,1/(R1-R2)*z-R2/(R1-R2)),1)
c[:] = 0
c[np.logical_and(np.abs(x[:,0])< R,np.abs(x[:,1])<R)] = 1
#for i in range(len(x[:,0])):
# z = (x[i,0])**2+(x[i,1])**2
# R1 = (0.7*R)**2
# R2 = (1.3*R)**2
# #s[i,1] = -0.01
# c[i] = min(max(0,1/(R1-R2)*z-R2/(R1-R2)),1)
# c[: ]
# c[i] = np.where(np.abs(x[i,0] < R)
fluid.set_concentration_cg(c)
......@@ -95,8 +98,8 @@ tic = time.time()
while t < tEnd :
#Fluid solver
fluid.implicit_euler(dt)
if (ii%11==0 and ii != 0):
fluid.adapt_mesh(1e-2,1e-4,5000)
#if (ii%11==0 and ii != 0):
# fluid.adapt_mesh(1e-2,1e-4,5000)
t += dt
#Output files writting
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))
......
L = 0.05;
H = 0.6;
H = 0.05;
y = 0;
lc = 0.001;
Point(1) = {-L, H, 0};
Point(2) = {-L, -H, 0};
Point(3) = {L, -H, 0};
Point(4) = {L, H, 0};
Point(1) = {-L, H, 0,lc};
Point(2) = {-L, -H, 0,lc};
Point(3) = {L, -H, 0,lc};
Point(4) = {L, H, 0,lc};
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
......@@ -20,20 +20,4 @@ Physical Line("Bottom") = {2};
Physical Line("Lateral") = {1,3};
Physical Line("Top") = {4};
Physical Surface("Domain") = {1};
Physical Point("PtFix") = {1};
Field[1] = Attractor;
Field[1].NodesList = {5};
Field[2] = Threshold;
Field[2].DistMax = 0.08;
Field[2].DistMin = 0.02;
Field[2].LcMax = 0.008;
Field[2].LcMin = 0.0005;
Field[2].IField = 1;
Background Field = 2;
Mesh.CharacteristicLengthFromPoints = 0;
Mesh.CharacteristicLengthExtendFromBoundary = 0;
Mesh.MshFileVersion = 2;
\ No newline at end of file
Mesh.MshFileVersion = 2;
Supports Markdown
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