Commit 508ed905 authored by Matthieu Constant's avatar Matthieu Constant
Browse files

new equations tested on two testcases

parent 82501b31
......@@ -256,11 +256,12 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
grad_shape_functions(dxidx, dphi);
double *local_vector = &all_local_vector[N_SF*n_fields*iel];
double *local_matrix = &all_local_matrix[N_SF*N_SF*n_fields*n_fields*iel];
double dc[DIMENSION]={0}, du[DIMENSION][DIMENSION]={{0}}, dp[DIMENSION]={0};
double dc[DIMENSION]={0}, du[DIMENSION][DIMENSION]={{0}}, dp[DIMENSION]={0}, dcdq[DIMENSION]={0};
for (int i = 0; i< N_SF; ++i) {
for (int j = 0; j < DIMENSION; ++j) {
dc[j] += dphi[i][j]*porosity[el[i]];
dp[j] += dphi[i][j]*solution[el[i]*n_fields+DIMENSION];
dc[j] += dphi[i][j]*solution[el[i]*n_fields+DIMENSION+1];
dp[j] += dphi[i][j]*solution[el[i]*n_fields+DIMENSION];
dcdq[j] += dphi[i][j];
for (int k = 0; k < DIMENSION; ++k)
du[k][j] += dphi[i][j]*solution[el[i]*n_fields+k];
}
......@@ -289,7 +290,7 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
double tau[DIMENSION][DIMENSION];
for (int i = 0; i < DIMENSION; ++i)
for (int j = 0; j < DIMENSION; ++j)
tau[i][j] = du[i][j]-u[i]*dc[j]/c;
tau[i][j] = du[i][j]-u[i]*dc[j]/q;
double dphidp = 0, divu = 0;
for (int k = 0; k < DIMENSION; ++k){
dphidp += dphii[k]*dp[k];
......@@ -326,7 +327,7 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
double udtau = 0;
double dphiidtau = 0;
for (int j = 0; j < DIMENSION; ++j) {
dtau[j] = dphij[j]-phij*dc[j]/c;
dtau[j] = dphij[j]-phij*dc[j]/q;
dphiidphij += dphii[j]*dphij[j];
udtau += u[j]*dtau[j];
dphiidtau += dphii[j]*dtau[j];
......@@ -336,9 +337,13 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
LOCAL_MATRIX(U,P) += -jw*dphii[j]*phij;
LOCAL_MATRIX(P,U) += jw*phii*dphij[j];
double utau =0;
double dphisigmadq = 0;
double dutaudq = 0;
for (int k = 0; k < DIMENSION; ++k) {
int V = k;
utau += u[k]*tau[j][k];
dutaudq += u[k]*u[j]*(dc[k]*phij/(q*q)-dcdq[k]/q);
dphisigmadq += 2*mu*dphii[k]*0.5*(u[k]*dc[j]*phij/(q*q)+u[j]*dc[k]*phij/(q*q)-u[k]*dcdq[j]/q-u[j]*dcdq[k]/q);
//LOCAL_MATRIX(U,V) += -jw*rho*u[j]*dphii[k]*phij/c+jw*2*mu*dphii[k]*dtau[j]*0.5;
LOCAL_MATRIX(U,V) +=
jw*rho*phii*phij*tau[j][k]/q
......@@ -347,7 +352,7 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
}
//LOCAL_MATRIX(U,U) += jw*mu*2*0.5*dphiidtau + jw*rho*(phii*phij/dt-phij*dphii[j]*u[j]/c);
LOCAL_MATRIX(U,U) += jw*mu*2*0.5*dphiidtau + jw*rho*phii*(phij/dt+phij/q*divu+udtau/q);
LOCAL_MATRIX(U,Q) += -jw*rho*phii*(u[j]*divu*phij/(q*q)+utau*phij/(q*q));
LOCAL_MATRIX(U,Q) += jw*rho*phii*(-u[j]*divu*phij/(q*q)-utau*phij/(q*q)+dutaudq/q)+jw*dphisigmadq;
}
LOCAL_MATRIX(P,P) += jw*epsilon*dphiidphij;
LOCAL_MATRIX(Q,Q) += -jw*phii*phij;
......
......@@ -58,7 +58,7 @@ def genInitialPosition(filename, r, rout, rhop, compacity) :
p.write_vtk(filename,0,0)
outputdir = "output"
outputdir = "outputTestNew"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
......
......@@ -45,7 +45,7 @@ def genInitialPosition(filename, r, ox, oy, lx, ly, rhop) :
p.write_vtk(filename,0,0)
outputdir = "output"
outputdir = "outputTest"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
......@@ -69,7 +69,7 @@ epsilon = 1e-4 # stabilization parameter
print("eps = %g" % ( epsilon))
#Enable the use of lmgc90
use_lmgc = True
use_lmgc = False
#Object particles creation
......@@ -90,12 +90,12 @@ outf = 20 # number of iterations between o
#Format: strong_boundaries = [(Boundary tag, Fluid field, Value)
strong_boundaries = [("Box",0,0),("Box",1,0),("Top",0,0.),("Top",1,0.),("Top",2,0.)]
fluid = fluid.fluid_problem("mesh.msh",g,nu*rho,rho,epsilon,strong_boundaries)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),True)
fluid.export(outputdir,0,0)
tic = time.clock()
#Computation loop
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),True)
while t < tEnd :
#Fluid solver
fluid.implicit_euler(dt)
......@@ -108,7 +108,7 @@ while t < tEnd :
print("NSUB", nsub,"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(p.r()))
#Contact solver
for i in range(nsub) :
tol = 1e-6
tol = 1e-8
p.iterate(dt/nsub, forces, tol)
#Localisation of the grains in the fluid
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
......
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