glwagner
commented
3 years ago One extra consideration is that we also allow users to specify diffusive fluxes across boundaries. This is especially important for geophysical problems at LES scales and larger, where its appropriate to employ a "wall model" to predict momentum and tracer fluxes at boundaries where there's an unresolved or partially-resolved turbulent boundary layer (rather than prescribing a particular value or gradient). That said, I think if we are able to specify the gradients of a field across a boundary it will likely be straightforward to extend that implementation to specifying fluxes.
The paper
"Moving from momentum transfer to heat transfer – A comparative study of an advanced Graetz-Nusselt problem using immersed boundary methods" by Lu et al. (2019)
may also be relevant. Their conclusion is a bit confusing. They state
In all simulations, excellent agreement are reached between CFM and DFM results, with the deviation being below 10%.
which suggests that accuracy may not be an important factor in deciding which method to use. But the next sentence is
Considering the nature of capturing the discontinuity at the fluid-solid interface, DFM might offer a more accurate result, which however requires more follow-up simulations to give a solid investigation.
which is difficult to interpret. I suppose all they can say is that their results are similar to one another, but they cannot say which one is more accurate (and perhaps it doesn't matter which method is more accurate in their case, if both methods return similar results).
That said, I think time-step considerations are really important, and seem like a good reason to choose DFM over CFM.
Balaras (2004)'s pressure equation is
where, crucially, Ω is the computational domain (irrespective of the immersed boundary). So following Balaras (2004) would mean not modifying the pressure solver? That's good news!
Does it make sense to first implement an algorithm that assumes the boundary coincides with the grid?
I am not super worried about the algorithmic changes require to implement a DFM immersed boundary. I think we can isolate the immersed boundary implementation from the rest of the code and interface with the time-stepping routines with a single function call that calculates the discrete forcing and applies a correction to the velocity field. Using multiple dispatch, this change to the algorithm will have no cost to simulations that don't use immersed boundaries. The main challenge I think is figuring out abstractions that make it easy to create immersed boundaries and assign boundary conditions (@ali-ramadhan and my job to figure this out) --- and we'd have this challenge for any immersed boundary implementation, whether CFM or DFM.
johncmarshall54
weymouth
whitleyv
christophernhill
navidcy
wenegrat
I am also getting a small velocity within the cylinder. I can plot the velocity within the cylinder at a couple of points compared to a rough estimate of the RK3 last stage pressure correction, and they match pretty well. I'm fairly certain that the velocity errors within the IB scale with the NH pressure correction, so as long as the time stepping is kept reasonably small, then the velocity errors within the cylinder should be negligible. Is this the correct interpretation of the velocity error due to pressure correction? @weymouth @glwagner
francispoulin
hdrake
@christophernhill @glwagner @weymouth I have been looking into what is the best way to implement an immersed boundary method here, and I hope to be able to model complex topography and possibly moving boundaries with it, for at least Dirichlet and Neumann boundary conditions, within the next year.
The immersed boundary method (IBM) is a way to deal with complex topography without a complicated mesh or coordinate transformations. Instead of conforming the mesh to the fluid domain, a cartesian grid is generated over the whole area. Then, an added forcing term incorporates the boundary conditions into the equations. Mittal and Iaccarino (2003), outline several ways we could do this, which all fall into two large categories:
Continuous Forcing: We add a forcing term to the continuous equations, discretize the new system, and solve as normal. If the IB is set to coincide with the mesh nodes, then you get a “stairstep” like boundary, whereas if the IB intersects the mesh arbitrarily, you must distribute the forcing to several nearby nodes. Either way, it may not give you a sharp enough boundary, and you get stability constraints on rigid boundaries.
Discrete Forcing: We discretize, then determine what the forcing needs to be at each time step to satisfy the boundary conditions. This method usually takes a predictor-corrector type approach as described by Balaras (2004). You would use a predictor step to interpolate the correct values at the boundary with nearby nodes. Then, you can solve for the necessary forcing in the discretized time stepper. Finally, you recompute the true solution for the next time step, with the correct forcing term. This method does not have the stability constraints and can create a sharper boundary, but implementation will involve more changes to create this routine.
I'm leaning towards the discrete forcing, since it seems like continuous just won't work well for rigid boundaries. I haven't fully worked out how the predictor-corrector idea will work with the pressure solver, etc.
WRF uses Discrete forcing in their LES IBM code, while PALM uses a masking method that involves stair step representation for complex topography.