Postprocess extraction of gradients at boundaries

[SA8416]

Hi,

I have a 2.5D (cylindrical axisymmetric) domain, over which I solve Laplace's equation for heat transfer. I have obtained the temperature gradient by doing:
∇T = nodeflux(sys, T) , which returns a matrix with the same dimensions as T.

I would like to obtain the temperature gradient on the West boundary, but am unsure what indexing yields this. I have surmised that grid[BFaceNodes], grid[BFaceEdges] or grid[BEdgeNodes] contain the indices of all boundary nodes.

Further to this, once I have obtained ∇T_west, would I be able to simply input this as a boundary condition to another PDE on the same domain? For example:
boundary_neumann!(f, p, bnode; species = 1, region = 4, value = ∇T_west)

[j-fu]

Hi, thanks for your interest in the package!

Just to avoid misunderstandings:

• If your domain has cylindrical symmetry, you are using the corresponding feature from VoronoiFVM.jl ?
• the other PDE also assumes that the domain has cylindrical symmetry ?
• Is it true that the "western boundary" in your case is not at r=0 ? ( otherwise, the cylindrical treatment automatically assumes symmetry boundary conditions)
• How does the coupling of the temperature with the other process write in the strong formulation of your PDEs ?

[SA8416]

Thanks for the reply!

1. Yes I believe I am using the correct feature:

R = collect(r1:dr:r2)
Z = collect(0:dz:z2)
grid = VoronoiFVM.Grid(R, Z)
circular_symmetric!(grid)

2. Yes all PDEs I am solving assume cylindrical symmetry and use the same grid.
3. Yes the western boundary is not at r=0. It is set at a non-zero inner-radius of r1, though the length scales are on the order of microns so there may be issues with the 1/r term if VoronoiFVM doesn't include automatic numerical conditioning.
4. The strong formulation of the other PDE uses Darcy's law for pressure:
   ∇⋅((-ε*ρ*κ/μ)∇p) = 0 to find the velocity distribution, but can be simplified to Laplace's equation for the time being: ∇²p = 0
   I believe you have a similar example here https://j-fu.github.io/VoronoiFVM.jl/stable/nbhtml/problemcase/ but I have been unable to find the actual code that defines the problem. The boundary condition for Laplace's equation for pressure is -(κ/μ)*(∂p/∂r) = -(k/(ρ*hfg))*(∂T/∂r) on the western boundary, which is why I need the temperature gradient along this boundary.

As I am attempting to integrate VoronoiFVM.jl into my current code, each PDE is defined and solved seperately, i.e. a seperate VoronoiFVM.System is created and solved for each.

[j-fu]

Sorry for coming back late:

If you have $\nabla T$ on the whole grid, you can extract the x component ∇T_x as a function on the whole grid.
which gives you the flux in normal direction of your west boundary (assumed it is parallel to the z axis).
If you work on the same grid, you can make this visible to the other PDE.

For that PDE you can define a bcondition where you use

boundary_neumann!(f, p, bnode; species = 1, region = 4,
value = ∇T_x[bnode.index])

[SA8416]

Thank you!
I've implemented it this and it seems to work.
Yes the flux ∇T · n is in the direction of the outward-facing normal direction of the west boundary (4), but it is parallel to the r axis and perpendicular to the z axis.
Just to be sure, I have attached a diagram of the domain I am modelling (in gray), where the flux ∇T · n is represented by the red arrows.

I believe this is the same domain as the one given here?https://github.com/j-fu/VoronoiFVM.jl/blob/faba3b6edbaf9d090feb3df60a06f2a4a14d5aaf/examples/Example203_CoordinateSystems.jl#L151C20-L151C20

[j-fu]

Yes, this seems to fit.

[SA8416]

Further to this, is there a way to refer to the index of the boundary values of T?
I have several water property functions that require the values of T along the same Western boundary.

[j-fu]

Sorry, I lost a bit the plot...

If it hasn't been solved yet: I am not sure if I properly understand your question. If it is about the index of a boundary node in the grid, you can use bnode.index in oder to address a value in a global array. Otherwise: can you be a bit more spevfic ?

[SA8416]

No problem at all.

Yes I use bnode.index to refer to boundary nodes in any VoronoiFVM function that loop over bnode, but once I solve the PDE for T, I have several custom functions that depend on the value of T at the western boundary. I have found a workaround that seems to suffice:

function create_boundary_index(grid,ibc)
    sub = subgrid(grid,[ibc],boundary=true)
    v = [i for i = 1:num_nodes(grid)]
    subv = view(v,sub)
    ibc_vec = sort(subv)
    return ibc_vec
end

index_west = create_boundary_index(grid,4)

which gives index_west which contains all the boundary node indices at the west boundary.
I have not rigorously checked this workaround, so I was wondering if you had a suggestion for a better way to approach this?

On another note, do you have any examples of how VoronoiFVM treats spatially varying transport coefficients?

[j-fu]

No really good examples for this (note to self...). What is the characteristics of the spatial variation ? Different constant values in different subdomains or relatively smooth variation ?

[SA8416]

Thanks for your prompt reply.
No I've seen the examples you have for different constant values in separate subdomains. I'm referring to smooth variations throughout the entire domain, i.e. each node has a different value for the transport coefficient.

I have tried averaging the diffusion coefficient as such:

function c_flux(F, c, edge)
   for ispec = 1:Nspec
      De_ave = (De[ispec,1] + De[ispec,2])/2
      F_dif = De_ave*(c[ispec, 1] - c[ispec, 2])
      F[ispec] = F_dif
   end
end

for spatially-varying diffusion coefficients, and have similar formulations for spatially varying reaction rate constants for my reaction terms. I was wondering if you had any examples of these.

[j-fu]

No, and in fact it would be good to have an example.
While your approach certainly works, one can argue that it is more accurate to use an harmonic average of the diffusion coefficient. Basically, this would be the outcome of the approach in the Eymard/Fuhrmann/Gärtner paper cited in the docs. So yeah it makes sense to have such an example to make the point....