Flux bc1
From Gerris
Contents |
Revision of Cut-Cells Flux Boundary Condition
Problem statement
If you use non-homogeneous Neumann boundary condition on solid surfaces (cut cells), or generally use filed dependent flux boundary condition) on solid surfaces, you easily figure out Gerris problems in this regard (BC for passive tracers or momentum components). Based on applied case un-reasonable results, heavily mesh-dependent results or divergence of diffusion solver (usually by run-time error) could be observed. It should be noted that such boundary conditions are not formally used in simple fluid flow problem and due to this this obvious bug is not fixed so far (gerris-snapshot-080130 is checked).
The problem has several sources: the main source of error is because Gerris does not correctly compute flux-bc across solid interfaces at each time step, what Gerris done is computing flux based on user defined function or fixed defined user and add its contribution in rhs of diffusion equation (and residual during relaxation), while at least this flux should be multiplied by time step and area of cutting solid.
Furthermore the “Beta” blending factor (switching between 2nd and 1nd order time accuracy) is not correctly taken into account. This is because, unlike Dirichlet BC, Gerris treat flux BC explicitly during relaxation (problematic for field dependent BC), so if we accept problems of this treatment, “Beta” factor should not be considered for flux-BC while Gerris consider them (see: “diffusion_rhs”). The better treatment is to consider flux BC implicitly (and considering blending factor as well) during relaxations.
As stated above applying field dependent BC is problematic in Gerris (instability/accuracy issues), particularly when time step sizes are large with respect to stability of explicit diffusion solution. Note that for Dirichlet BC or homogeneous flux BC none of these problems happen and Gerris works fine.
We perform solution for this problem in this note and present validation against exact solution.
Another important issue is time step of diffusion solver. Of course Gerris solve its implicitly and so with sufficient number of relaxations instability is not issue, but large time step (with respect to stability bound of explicit method) leads to time inaccuracy. Since Gerris does not consider any criterion for diffusion terms, time accuracy is a-priori un-known (not predictable). In convection dominated flow it is not usually matter since stability bound of convection (CFL) is usually an order of magnitude smaller than stability bound of explicit diffusion solver. But for diffusion dominated flow (creeping flow, sub-surface flow, conduction in solids), where velocity filed is small (or zero) it leads to problem. We suggest to add time step criteria for diffusion solver, e.g. using time step bound of diffusion solver as 20 times of stability bound of explicit diffusion solver (based on author experience it is reasonble).
Important Note: another problematic issue of Gerris diffusion solver which is not related to flux BC only is its default blending factor “Beta=0.5”, theoretically it could lead to 2nd order accuracy and is critical bound for stability, i.e., “Beta” smaller than “0.5” leads to instability. But to author knowledge factor “0.5” is derived for pure structured or unstructured grids and not for immersed boundary grids, so it could leads to instability particularly in vicinity of cut-cells (author personally observe such behavior), so we strongly recommend a factor strictly larger than “0.5”, e.g. “0.51”.
What is Field Dependent BC
We mean filed dependent BC to cover BCs that are function of filed variable, e.g. when we solve for temperature, the temperature boundary condition be itself depend on temperature. This is mainly happen in flux BCs and to authors knowledge there is not field dependent pure Dirichlet problems (solvability could be issue). In practice convective or third kind boundary condition is very common in heat transport phenomena (sometimes is called as Rubin BC). Due to imperfect match of materials interfaces, interface works as a thermal barrier and could be fairly modeled by a convective boundary condition (thermal resistance is proportional to inverse of interfacial heat transfer coefficient). This boundary condition is defined as follows:
Flux = h (T-T_a)
Where h is interfacial heat transfer coefficient and T_a is ambient temperature.
Another famous field dependent BC is radiative BC which is defined as follows:
Flux = sigma (T^4-T_a^4)
Fixing Problem
The first correction is in function “gfs_surface_bc_bc”, corrected code is as follows (note that since time step is required in this function and is not accessible formally, we define a global variable in “solid.h” header, gdouble dt_temp;, and load its in “advance_tracers” function, you should do its yourself):
static void gfs_surface_bc_bc (FttCell * cell, GfsSurfaceGenericBc * b)
{
GfsSurfaceBc * bc = GFS_SURFACE_BC (b);
gdouble val = gfs_function_value (bc->val, cell);
if (gfs_function_value (bc->type, cell) > 0.) {
cell->flags |= GFS_FLAG_DIRICHLET;
GFS_STATE (cell)->solid->fv = val;
}
else {
//rt@modif
gdouble h = ftt_cell_size (cell);
gdouble solid_area;
FttVector n;
/*
get averaged unit normal vector in mixed cell
multiplied by area of solid bc
Warning: this should be scaled by "h*h" in 2D
and "h*h*h" in 3D to be real value
*/
gfs_solid_normal (cell, &n);
/*
compute area of solid boundary across its
flux is transported
Note that we do not use correct scaling, correct
value need one more scaling with h, but we do not
consider its to keep code consistence with other palce
where this value is used. Using correct value just
increases number of floating point operations
*/
#if FTT_2D
solid_area = sqrt(n.x * n.x + n.y * n.y);
#else
solid_area = h*sqrt(n.x * n.x + n.y * n.y + n.z * n.z);
#endif
cell->flags &= ~GFS_FLAG_DIRICHLET;
GFS_STATE (cell)->solid->fv = val*dt_temp*solid_area;
}
}
The second correction is related to blending factor and implicit treatment of flux BC, the changed functions are as follows, note that implementation of this issue (implicit treatment of flux BC) could not be performed uniquely but now we select a method which be consistent with current instruction of Gerris (solvers and IO ports) so it be quit general and minor code modification be required. For this purpose we should apply flux-BC, "gfs_domain_surface_bc", (plus considering correct blending) prior to each call of “gfs_diffusion_residual” and also fix issue related to blending factor in "diffusion_rhs" and "diffusion_residual" function (hope i do not forget to state any other required changes).
So corrected version of "diffusion_rhs" is:
static void diffusion_rhs (FttCell * cell, RelaxParams * p)
{
gdouble a, f, h, val;
FttCellNeighbors neighbor;
FttCellFace face;
if (GFS_IS_MIXED (cell)) {
a = GFS_STATE (cell)->solid->a*GFS_VARIABLE (cell, p->dia);
if (((cell)->flags & GFS_FLAG_DIRICHLET) != 0)
f = gfs_cell_dirichlet_gradient_flux (cell, p->u, -1, GFS_STATE (cell)->solid->fv);
else
// for neumann or user defined neumann bc
// neumann value or function value is solid->fv
f = GFS_STATE (cell)->solid->fv/(1.+p->beta);
}
else {
a = GFS_VARIABLE (cell, p->dia);
f = 0.; /* Neumann condition by default */
}
h = ftt_cell_size (cell);
val = GFS_VARIABLE (cell, p->u);
face.cell = cell;
ftt_cell_neighbors (cell, &neighbor);
for (face.d = 0; face.d < FTT_NEIGHBORS; face.d++) {
GfsGradient g;
face.neighbor = neighbor.c[face.d];
/* @@@ rt@comment
Note that the face gradient is scaled by "h" of cell and so
you should devide its by "h" to get real gradient
*/
gfs_face_gradient_flux (&face, &g, p->u, -1);
f += g.b - g.a*val;
}
corrected version of "diffusion_residual"
static void diffusion_residual (FttCell * cell, RelaxParams * p)
{
gdouble a;
GfsGradient g = { 0., 0. };
gdouble h;
FttCellNeighbors neighbor;
FttCellFace face;
h = ftt_cell_size (cell);
if (GFS_IS_MIXED (cell)) {
a = GFS_STATE (cell)->solid->a*GFS_VARIABLE (cell, p->dia);
if (((cell)->flags & GFS_FLAG_DIRICHLET) != 0)
g.b = gfs_cell_dirichlet_gradient_flux (cell, p->u, -1, GFS_STATE (cell)->solid->fv);
else
g.b = GFS_STATE (cell)->solid->fv * p->beta;
}
else
a = GFS_VARIABLE (cell, p->dia);
face.cell = cell;
ftt_cell_neighbors (cell, &neighbor);
for (face.d = 0; face.d < FTT_NEIGHBORS; face.d++) {
GfsGradient ng;
face.neighbor = neighbor.c[face.d];
/* @@@ rt@comment
Note that the face gradient is scaled by "h" of cell and so
you should devide its by "h" to get real gradient
*/
gfs_face_gradient_flux (&face, &ng, p->u, -1);
g.a += ng.a;
g.b += ng.b;
}
/* @@@ rt@comment
note that one "h" is scale gradient to correct
value and another is due to area/volume relation
*/
a *= h*h;
g_assert (a > 0.);
g.a = 1. + g.a/a;
g.b = GFS_VARIABLE (cell, p->rhs) + g.b/a;
GFS_VARIABLE (cell, p->res) = g.b - g.a*GFS_VARIABLE (cell, p->u);
}
Corrected version of "gfs_diffusion_cycle" is as follows
void gfs_diffusion_cycle (GfsDomain * domain,
guint levelmin,
guint depth,
guint nrelax,
GfsVariable * u,
GfsVariable * rhs,
GfsVariable * dia,
GfsVariable * res)
{
guint n;
GfsVariable * dp;
RelaxParams p;
gpointer data[2];
g_return_if_fail (domain != NULL);
g_return_if_fail (u != NULL);
g_return_if_fail (rhs != NULL);
g_return_if_fail (dia != NULL);
g_return_if_fail (res != NULL);
dp = gfs_temporary_variable (domain);
dp->component = u->component;
/* compute residual on non-leafs cells */
gfs_domain_cell_traverse (domain,
FTT_POST_ORDER, FTT_TRAVERSE_NON_LEAFS, -1,
(FttCellTraverseFunc) gfs_get_from_below_intensive, res);
/* relax top level */
gfs_domain_cell_traverse (domain,
FTT_PRE_ORDER, FTT_TRAVERSE_LEVEL, levelmin,
(FttCellTraverseFunc) gfs_cell_reset, dp);
p.maxlevel = levelmin;
p.u = dp->i;
p.res = res->i;
p.dia = dia->i;
p.beta = my_beta;
for (n = 0; n < 10*nrelax; n++) {
gfs_domain_homogeneous_bc (domain,
FTT_TRAVERSE_LEVEL | FTT_TRAVERSE_LEAFS,
levelmin, dp, u);
gfs_domain_cell_traverse (domain, FTT_PRE_ORDER,
FTT_TRAVERSE_LEVEL | FTT_TRAVERSE_LEAFS,
levelmin,
(FttCellTraverseFunc) diffusion_relax, &p);
}
/* relax from top to bottom */
for (p.maxlevel = levelmin + 1; p.maxlevel <= depth; p.maxlevel++) {
/* get initial guess from coarser grid */
gfs_domain_cell_traverse (domain,
FTT_PRE_ORDER, FTT_TRAVERSE_LEVEL | FTT_TRAVERSE_NON_LEAFS,
p.maxlevel - 1,
(FttCellTraverseFunc) get_from_above, dp);
for (n = 0; n < nrelax; n++) {
gfs_domain_homogeneous_bc (domain,
FTT_TRAVERSE_LEVEL | FTT_TRAVERSE_LEAFS,
p.maxlevel, dp, u);
gfs_domain_cell_traverse (domain, FTT_PRE_ORDER,
FTT_TRAVERSE_LEVEL | FTT_TRAVERSE_LEAFS, p.maxlevel,
(FttCellTraverseFunc) diffusion_relax, &p);
}
}
/* correct on leaf cells */
data[0] = u;
data[1] = dp;
gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
(FttCellTraverseFunc) correct, data);
gfs_domain_bc (domain, FTT_TRAVERSE_LEAFS, -1, u);
// rt@modif
gfs_domain_surface_bc (domain, u);
/* compute new residual on leaf cells */
// this call account field dependent bc
gfs_diffusion_residual (domain, u, rhs, dia, res);
gts_object_destroy (GTS_OBJECT (dp));
}
Note that better treatment for linear filed dependent flux-BC is to include its contribution in relaxation (better convergence due to increase diagonal dominancy!), for non-linear case (e.g. radiative heat transfer, happen in high temperatures and liquid metal treatments) it is not as easy but it could be possible to linearize its and add some contributions to diagonal coefficient, e.g., For radiation BC we have: Flux = sigma * (T^4 – T_a^4) = sigma * (T-T_a) * (T+T_a) * (T^2+T_a^2) And so taking first term implicitly (like a convective BC)
Time Step Determination
In our examples we consider pure diffusion problem (zero velocity filed ), we take time step size as 10-20 times of stability bound of explicit diffusion which is computed by the following relation (modification of “gfs_simulation_set_timestep”):
void gfs_simulation_set_timestep (GfsSimulation * sim)
{
...
//rt@modif : stability of diffusion
domain = GFS_DOMAIN (sim);
h_min = G_MAXDOUBLE; /* h_min */
gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
(FttCellTraverseFunc) find_h_min, &h_min);
// assume diffusion coefficient is equal to 1.
dt_diffusin = 10. * h_min*h_min/(2. * FTT_DIMENSION * 1.);
sim->advection_params.dt = MIN(sim->advection_params.dt, dt_diffusin);
tnext = G_MAXINT;
i = sim->events->items;
while (i) {
GfsEvent * event = i->data;
if (t < event->t && event->t < tnext)
tnext = event->t + 1e-9;
i = i->next;
}
....
}
Validation against exact solution
Exact solution
There are a little number of exact solution for heat flow across complex geometry, for validation purpose we take pure conduction inside an sphere which is placed in environment with unit interfacial heat transfer coefficient and zero ambient temperature.
Source code to compute exact solution is as follows:
#ifndef __EXACT_SOLUTION_H__
#define __EXACT_SOLUTION_H__
//#include "exact_solution.h"
/* exact solution of transient heat conduction within sphere
constant initial temperature and convective bc with ambient
temperature zero
k : conductivity
h : interfacial heat transfer coefficient
R : radius
T0: initial temperature
r : desired location
time : desired time
*/
double sphere_exact_solution_third_kind_bc
(const double k, const double h, const double R, const double T0, const double r, const double time) {
/* first six roots of mu cot(mu) + hR = 1, for hR = 0.4 */
double mu[6] = { 1.0528, 4.5822, 7.7770, 10.9403, 14.0946, 17.2440};
double hR = h*R;
double T = 0.;
int n;
if(fabs(hR-0.4)>1.e-3) {
printf("\n exact solution is not available for your case \n");
exit(0);
}
for(n=0; n<6; n++)
T += 2.*h*R*R*T0/r*
(mu[n]*mu[n]+(hR-1.)*(hR-1.))/(mu[n]*mu[n]*(mu[n]*mu[n]+hR*(hR-1.)))*
sin(mu[n])*sin(mu[n]*r/R)*exp(-1.*k*mu[n]*mu[n]*time/R/R);
return T;
}
#endif
Numerical Results
For validation purpose we take an sphere with radius 0.4 and thermal diffusion coefficient unity and compute temperature history at two probe locations: at r=0.3 inside, and r=0.4 on the surface of sphere. Note that we have used tri-linear interpolation (like Hex finite element shape function) to compute numerical results at probe location. The blending factor "Beta" is taken equal to "0.51" in this example. The following input file shows details of problem (note that we did not use fine mesh for performance issue, this case was completed below 1 minute):
1 0 GfsSimulation GfsBox GfsGEdge {} {
Time { end = 0.5 }
Solid (-x*x - y*y - z*z + 0.4*0.4)
# RefineSolid 4
GfsRefine 4
VariableTracer {} T
SourceDiffusion {} T 1.
Init {} { T = 1.0 }
SurfaceBc T Neumann (-1.*T)
OutputTime { istep = 1 } stderr
# OutputTecplot { istep = 10 } stdout
}
GfsBox {}
Figure 2 and 3 show temperature history for probe location (r=0.3) and 0.4 respectively, it is obvious that in spite of coarse mesh numerical results fairly match with exact solution.
To grid resolution study, we use one additional reniement, i.e., solve for "GfsRefine 5", the following figures show results of this experiment. Plots show convergence to exact solution that confirm our revision.
It should be noted that we have almost 2nd order convergence in L2-norm, but local convergence near sloped boundaries are not as well, author belive that this is probably due to in-accuracy of computation of solid area in Gerris.
Convection-Diffusion Example
In this example a convection-diffusion problem is considered with a constant heat flux boundary condition on solid boundary. Input GFS file is as follows (only for grid resolution study number of refinement is changed (5, 6 and 7 level refinement of domain).
1 0 GfsSimulation GfsBox GfsGEdge {} {
Time { end = 1.2 }
# Refine 7
Refine 6
# Refine 5
Solid (x*x + y*y - 0.25*0.25)
VariableTracer {} T
SourceDiffusion {} T 0.001
SurfaceBc T Neumann 0.01
Init {} { U = 1 }
SourceDiffusion {} U 0.00078125
SourceDiffusion {} V 0.00078125
OutputTime { istep = 10 } stderr
OutputBalance { istep = 10 } stderr
OutputProjectionStats { istep = 10 } stderr
OutputTiming { start = end } stderr
}
GfsBox {
left = Boundary {
BcDirichlet U 1
}
right = BoundaryOutflow
}
Figure 6 and 7 shows temperature distribution at t=1.0 for old flux boundary condition and corrected version for various grid resolution. Note that scaling and range of colors in legend is automatically generated by “Tecplot” based on minimum and maximum temperature.
Plots show that old BC show instability near interface and we have not convergence by increasing grid resolution. In contrast revised flux BC shows good convergence without instability.