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

symmetric part viscous tensor

parent 1a689a76
......@@ -241,57 +241,66 @@ static void f_boundary(WeakBoundary *wbnd, FluidProblem *problem,const double *n
const int aid = wbnd->aid;
if (wbnd->type == BND_WALL) {
double p = s[P];
double u[D], du[D][D], dp[D], dun[D];
double u[D], du[D][D], dp[D], dun[D], dunt[D];
for (int iD = 0; iD < D; ++iD) {
u[iD] = s[U+iD];
dp[iD] = ds[P*D+iD];
dun[iD] = 0;
dunt[iD] = 0;
for (int jD = 0; jD < D; ++jD) du[iD][jD] = ds[(U+iD)*D+jD];
}
double un = 0;
double dcn = 0;
double s_c = 0;
for (int i = 0; i < D; ++i) {
un += u[i]*n[i];
dcn += dc[i]*n[i];
for (int j = 0; j < D; ++j) dun[i] += du[i][j]*n[j];
s_c += vid<0?0:(u[i]-data[vid+i])*n[i];
for (int j = 0; j < D; ++j) {
dun[i] += du[i][j]*n[j];
dunt[i] += du[j][i]*n[j];
}
}
double h = problem->element_size[eid];
f[P] = (pid<0 ? 0 : -(s[P]-data[pid])/h*problem->taup[eid]);
for (int id = 0; id < D; ++id) {
f[U+id] = (pid<0?p:data[pid])*n[id] + (un> 0 ? rho*u[id]*un/c : 0) + (vid<0?0:(1+D)*(s[U+id]-data[vid+id])/(D*h)*(mu+problem->tauc[eid])*N_N) - mu/c*(dun[id]-u[id]/c*dcn);
f[U+id] = (pid<0?p:data[pid])*n[id] + (un> 0 ? rho*u[id]*un/c : 0) + (vid<0?0:(1+D)/(D*h)*mu*(1+problem->tauc[eid])*N_N*((u[id]-data[vid+id])+s_c*n[id])) - mu/c*(dun[id]-u[id]/c*dcn+dunt[id]-un/c*dc[id]);
}
//printf("%s:velocity=[%g,%g]\n",wbnd->tag,data[vid],data[vid+1]);
}
else {
double p = pid<0 ? s[P] : data[pid];
double u[D], du[D][D], dp[D];
double u[D], du[D][D], dp[D], dun[D], dunt[D];
for (int iD = 0; iD < D; ++iD) {
u[iD] = s[U+iD];
dp[iD] = ds[P*D+iD];
dun[iD] = 0;
dunt[iD] = 0;
for (int jD = 0; jD < D; ++jD) du[iD][jD] = ds[(U+iD)*D+jD];
}
double divu = 0;
double tau[D][D];
double utau[D]={0.};
double un = 0;
double unext = 0;
double dcn = 0;
const double *uext = vid < 0 ? &s[U] : &data[vid];
double cext = cid<0?c:data[cid];
double s_c = 0;
for (int i = 0; i < D; ++i) {
un += u[i]*n[i];
unext += uext[i]*n[i];
divu += du[i][i];
dcn += dc[i]*n[i];
s_c += vid<0?0:(u[i]-data[vid+i])*n[i];
for (int j = 0; j < D; ++j) {
tau[i][j] = du[i][j]-u[i]*dc[j]/c;
utau[i] += u[j]*tau[i][j];
dun[i] += du[i][j]*n[j];
dunt[i] += du[j][i]*n[j];
}
}
double rho_ext =(problem->n_fluids == 1 ? rho : problem->rho[0]*data[aid]+problem->rho[1]*(1-data[aid]));
double h = problem->element_size[eid];
f[P] = unext;
for (int id = 0; id < D; ++id) {
f[U+id] = p*n[id] + (un> 0 ? rho*u[id]*un/c : rho_ext*unext*uext[id]/data[cid]);
f[U+id] = p*n[id] + (un> 0 ? rho*u[id]*un/c : rho_ext*unext*uext[id]/cext) + (vid<0?0:(1+D)/(D*h)*mu*(1+problem->tauc[eid])*N_N*((u[id]-uext[id])+s_c*n[id])) - mu/c*(dun[id]-u[id]/c*dcn+dunt[id]-un/c*dc[id]);
}
}
}
......
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