Commit 92e0df31 authored by Matthieu Constant's avatar Matthieu Constant
Browse files

old equations

parent a1c6bf02
Pipeline #4531 failed with stage
in 20 seconds
......@@ -193,7 +193,8 @@ static void node_force_f(FluidProblem *problem, double *f0, double *sol, double
}
double mass = problem->particle_mass[ip];
double vol = problem->particle_volume[ip];
double g = problem->g;
double rhop = mass/vol;
double g = problem->g*(rhop-problem->rho)/rhop;
double gamma;
if(problem->n_fluids == 2){
......@@ -419,7 +420,7 @@ static void f_wall_pressure(FluidProblem *problem,const double *n, double *f,con
f[P] = 0;
for (int id = 0; id < D; ++id) {
f[U+id] = p*n[id];
f[U+id] = 0*p*n[id];
}
if(problem->n_fluids == 2){
......@@ -545,7 +546,7 @@ static void fluid_problem_f(FluidProblem *problem, double *f0, double *f1, doubl
}
}
f0[P] = (c-cold)/dt;
f0[P] = (c-cold)/dt+divu;
if(problem->n_fluids == 2){
f0[A] = (a*c-aold*cold)/dt;
......@@ -554,12 +555,12 @@ static void fluid_problem_f(FluidProblem *problem, double *f0, double *f1, doubl
for (int i = 0; i < D; ++i) {
double dudt = (u[i]-uold[i])/dt;
double R = dudt + u[i]*drhodt/rho + dp[i]/rho + u[i]*divu/c + utau[i]/c + uugradrho[i]/rho + (i==1 ? -c*problem->g : 0.);
f0[U+i] = rho*dudt + u[i]*drhodt + (i==1 ? -c*rho*problem->g : 0.);
f0[U+i] = rho*dudt + u[i]*drhodt + (i==1 ? -0*c*rho*problem->g : 0.);
for (int j = 0; j < D; ++j) {
f1[(U+i)*D+j] = -rho*u[i]*u[j]/c + mu*(tau[i][j]+tau[j][i]) + (i==j ? -p : 0) + R*u[j]*rho*taup;
f1[(U+i)*D+j] = -rho*u[i]*u[j]/c + mu*(tau[i][j]+tau[j][i]) + (i==j ? -p : 0) + R*u[j]*rho*taup*0;
}
f1[(U+i)*D+i] += divu*tauc*rho;
f1[P*D+i] = -u[i] + taup*R;
f1[(U+i)*D+i] += 0*divu*tauc*rho;
f1[P*D+i] = -u[i]*0 + 1e-6*dp[i];
if(problem->n_fluids == 2){
f1[P*D+i] += taua*(divu+(c-cold)/dt)*u[i];
f1[A*D+i] = -u[i]*a + taua*Ra*u[i] + taup*R*a + da[i]*(1e-7+extraa);
......@@ -634,6 +635,7 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, double dt, d
dx[0][0]*dx[1][1] - dx[0][1]*dx[1][0]};
const double detbnd = sqrt(nn[0]*nn[0] + nn[1]*nn[1] + nn[2]*nn[2]);
const double n[3] = {-nn[0]/detbnd,-nn[1]/detbnd,-nn[2]/detbnd};
//printf("%g,%g,%g,%s\n",n[0],n[1],n[2],mesh->phys_name[iphys]);
#endif
for (int iqp = 0; iqp < N_LQP; ++iqp) {
double xi[D];
......@@ -877,7 +879,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);
compute_weak_boundary_conditions(problem, dt, all_local_vector, all_local_matrix);
//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);
for (int iel=0; iel < mesh->n_elements; ++iel) {
double *local_vector = &all_local_vector[local_size*iel];
......
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