Commit e75c932c authored by Matthieu Constant's avatar Matthieu Constant
Browse files

add viscous boundary terms

parent b9044c0d
......@@ -163,7 +163,7 @@ $\langlC\bm{f}_{q_p}^0,\psi_i\ranglC$&&&$=0$\\[1em]
\begin{align*}
\bm{f}_{\bm{u}}^0&=\rho_k\left( \dpartial{\ub}{t}+ \ub\dfrac{\nb\cdot\ub}{q}+\dfrac{\bm{A}_k\cdot\ub}{q}\right)-p\nb q_k-q_k\rho_k\bm{g}\\[1em]
\bm{f}_{\bm{u}}^1&= -q_kp\ib+\mu \bm{B}_k\\[1em]
\bm{f}_{\bm{ub}}&=\bm{B}_k\cdot\bm{n}_k+q_k p\bm{n}\\[1em]
\bm{f}_{\bm{ub}}&=-\mu\bm{B}_k\cdot\bm{n}_k+q_k p\bm{n}\\[1em]
\bm{f}_{q_k}^0&=\dpartial{q_k}{t}+\varepsilon\nb\cdot( q_k\rho_k\bm{g})\\[1em]
\bm{f}_{q_k}^1&=\varepsilon\left( q_k\nb p\right)-\ub\\[1em]
\bm{f}_{q_k b}&=-\varepsilon q_k\nb p\cdot\bm{n}_k+\ub\\[1em]
......
......@@ -249,14 +249,34 @@ static void f_inflow(FluidProblem *problem,const double *n, double *f,const doub
dp[id] = ds[P*D+id];
}
for (int ifluid = 0; ifluid < n_fluids; ++ifluid) {
double mu = problem->mu[ifluid];
int U = ifluid*(D+1);
int Q = ifluid*(D+1)+D;
double q = s[Q];
f[Q] = 0;
double dq[D] = {ds[Q*D+0],ds[Q*D+1]};
double u[D] = {s[U],s[U+1]};
double du[D][D] = {{ds[U*D+0],ds[U*D+1]},{ds[U*D+D+0],ds[U*D+D+1]}};
double oq2 = 1./fmax(q,1e-5);
double oq = 1./fmax(q,1e-2);
double divu = 0;
double tau[D][D];
double utau[D]={0.};
for (int i = 0; i < D; ++i) {
divu += du[i][i];
for (int j = 0; j < D; ++j) {
tau[i][j] = du[i][j]-u[i]*dq[j]*oq;
utau[i] += u[j]*tau[i][j];
}
}
for (int id = 0; id < D; ++id) {
f[Q] += s[U+id]*n[id];
f[Q] += -dp[id]*n[id]*epsilon*q;
f[U+id] = 0;//p*q*n[id];
for (int jd = 0; jd < D; ++jd){
f[U+id] -= mu*(tau[id][jd]+tau[jd][id])*n[jd];
}
}
}
f[P] = 0;
......@@ -275,14 +295,34 @@ static void f_outflow(FluidProblem *problem,const double *n, double *f,const dou
dp[id] = ds[P*D+id];
}
for (int ifluid = 0; ifluid < n_fluids; ++ifluid) {
double mu = problem->mu[ifluid];
int U = ifluid*(D+1);
int Q = ifluid*(D+1)+D;
double q = s[Q];
f[Q] = 0;
double dq[D] = {ds[Q*D+0],ds[Q*D+1]};
double u[D] = {s[U],s[U+1]};
double du[D][D] = {{ds[U*D+0],ds[U*D+1]},{ds[U*D+D+0],ds[U*D+D+1]}};
double oq2 = 1./fmax(q,1e-5);
double oq = 1./fmax(q,1e-2);
double divu = 0;
double tau[D][D];
double utau[D]={0.};
for (int i = 0; i < D; ++i) {
divu += du[i][i];
for (int j = 0; j < D; ++j) {
tau[i][j] = du[i][j]-u[i]*dq[j]*oq;
utau[i] += u[j]*tau[i][j];
}
}
for (int id = 0; id < D; ++id) {
f[Q] += s[U+id]*n[id];
f[Q] += -dp[id]*n[id]*epsilon*q;
f[U+id] = 0;//p*q*n[id];
for (int jd = 0; jd < D; ++jd){
f[U+id] -= mu*(tau[id][jd]+tau[jd][id])*n[jd];
}
}
}
f[P] = 0;
......@@ -302,13 +342,33 @@ static void f_wall_pressure(FluidProblem *problem,const double *n, double *f,con
}
double p = s[P];
for (int ifluid = 0; ifluid < n_fluids; ++ifluid) {
double mu = problem->mu[ifluid];
int U = ifluid*(D+1);
int Q = ifluid*(D+1)+D;
double q = s[Q];
f[Q] = 0;
double dq[D] = {ds[Q*D+0],ds[Q*D+1]};
double u[D] = {s[U],s[U+1]};
double du[D][D] = {{ds[U*D+0],ds[U*D+1]},{ds[U*D+D+0],ds[U*D+D+1]}};
double oq2 = 1./fmax(q,1e-5);
double oq = 1./fmax(q,1e-2);
double divu = 0;
double tau[D][D];
double utau[D]={0.};
for (int i = 0; i < D; ++i) {
divu += du[i][i];
for (int j = 0; j < D; ++j) {
tau[i][j] = du[i][j]-u[i]*dq[j]*oq;
utau[i] += u[j]*tau[i][j];
}
}
for (int id = 0; id < D; ++id) {
f[Q] += -dp[id]*n[id]*epsilon*q;
f[U+id] = q*p*n[id];
for (int jd = 0; jd < D; ++jd){
f[U+id] -= mu*(tau[id][jd]+tau[jd][id])*n[jd];
}
}
}
f[P] = 0;
......
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