DFN+fluid compressibility not using the correct reference volume

[Bruno Chareyre]

In the basic PFV scheme the incremental change of density of a pore fluid after a change of pore volume is dependent on dV/Vo where Vo is the reference pore space within a tetrahedral cell [1], i.e. V(tetrahedron)-V(spheres).

In DFNFlow "V" should reflect the fact that porosity is mainly due to cracks, i.e. it is much smaller than the numerical porosity calculated between spheres.
Typically V=opening*area+matrixPorosity*cellVolume
Currently it is using the same formula, hence overestimating the compressibility effect (because the DEM porosity is larger than a typical rock porosity). What should be the reference "opening" in above formula is to be clarified though, it has physical as well as numerical implications. Maybe slotInitialAperture is a candidate? Let you DFN people tell what it should be.

Note that the "dV" is less a problem because it is a difference (independent on the choice of the reference volume), so the incompressible scheme is not affected.

Bruno

[1] https://github.com/yade/trunk/blob/master/pkg/pfv/FlowEngine.ipp.in#L439

[Luc Scholtès (luc)]

Hi Bruno,

If we want to compute the volume associated to a crack, we should use directly crackAperture*crackArea as Vo in dV/Vo, no? Of course, this is not as straightforward because a cell is not associated to a unique crack in our scheme but still, why did you specifically propose V(tetrahedron)-crackAperture*crackArea-matrixPorosity? "aperture" defined in DFNFlow is the right candidate for crackAperture as it is defined as a function of the mechanical crack opening (particle-particle normal displacements) and hydraulic residual aperture. We should then probably change the way it is declared (locally in a function of DFNFlow right now).

Now, regarding the case of an pore contained in the intact rock matrix, are you suggesting that we need to "scale down" the porosity to adapt the scheme for rock materials? It makes sense of course but then we would need to define a scaling factor depending on the rock type and on the deviation of their microstructure with respect to the one of a granular assembly (e.g. a sandstone and a shale rock). The expression would then be something like Vo=scalingFactor*(Vtetrahedron-Vspheres) with scalingFactor=rockPorosity/granularAssemblyPorosity. Is that what you meant?

Luc

[Bruno Chareyre (bruno-chareyre)]

> we should use directly crackAperture*crackArea as Vo in dV/Vo, no?

That's more or less what I wanted to suggest. I can't explain why I wrote that strange formula.
However I also wanted to skip the calculation of crackAperture*crackArea, namely because there can be multiple cracks, as you point out.
So my idea was to store the initial value of the numerical porosity in one cell (call it V(t=0)) and to define the cracked volume as the difference V(t)-V(t=0). It makes multiple cracks less a problem I think.

However, an obvious problem of this is that it result in a reference equal to zero when the crack occurs, and we need to divide by this volume... that was the reason to introduce an initial amount of liquid via matrixPorosity (could also depend on slotInitialAperture). Just an idea anyway. The problem of fluid flow inside an advancing crack tip seems fundamentally tough... I can't imagine a perfect answer yet.

Bruno

[Robert Caulk (rcaulk)]

Ok, I finally got some time to digest this. Thank you for pointing this
out, Bruno, it may have pretty severe impacts on DFNFlow simulations.

In my opinion, we would just *add* the crack volume to the original pore
volume. Let me see if I can explain why I think this:

It seems we need to agree on how the crackArea and crackAperture are
computed for these volume calcs. Each cell is currently associated with a
unique crackArea [1], but I am unsure what the geometry of the current C++
lines represent [2]. I am particularly confused by [3]. I guess I need
someone to translate that for me.

If I were to imagine a plane that represents the crackArea, I would
probably change the existing code to look something like this:

if (calcCrackArea) {
    CVector edgeMidPoint = (
ed_it->first->vertex(ed_it->second)->point().point()
+ed_it->first->vertex(ed_it->third)->point.point() ) / 2;
    CVector v1 = edgeMidPoint - CellCentre1;
    CVector v2 = edgeMidPoint - CellCentre2;
    Real halfCrackArea = 0.25*sqrt(std::abs(cross_product(v1,
v2).squared_length()));
    cell1->info().crackVol += halfCrackArea*aperture
    cell2->info().crackVol += halfCrackArea*aperture;
}

Which would give us a crack plane that looks like this following poorly
drawn figure:

[image: Inline image 1]

In this scheme, the volume of the crack (crackVol in code above) within a
cell is simply the sum of all halfCrackArea_i*aperture_i where i indexes
the broken edge. Then we can add this volume to the original pore volume.
It's as if we are artificially creating volume space that is not reflected
by the packing geometry. Which is kind of what this all boils down to,
right? Or maybe I am totally lost. If you guys think this is worth testing,
I will implement, test, and report back.

Cheers,

Robert

[1] https://github.com/yade/trunk/blob/master/pkg/pfv/DFNFlow.cpp#L269
[2] https://github.com/yade/trunk/blob/master/pkg/pfv/DFNFlow.cpp#L263
[3] https://github.com/yade/trunk/blob/master/pkg/pfv/DFNFlow.cpp#L266

On Wed, Dec 6, 2017 at 2:52 AM, Bruno Chareyre <email address hidden>
wrote:

> > we should use directly crackAperture*crackArea as Vo in dV/Vo, no?
>
> That's more or less what I wanted to suggest. I can't explain why I wrote
> that strange formula.
> However I also wanted to skip the calculation of crackAperture*crackArea,
> namely because there can be multiple cracks, as you point out.
> So my idea was to store the initial value of the numerical porosity in one
> cell (call it V(t=0)) and to define the cracked volume as the difference
> V(t)-V(t=0). It makes multiple cracks less a problem I think.
>
> However, an obvious problem of this is that it result in a reference
> equal to zero when the crack occurs, and we need to divide by this
> volume... that was the reason to introduce an initial amount of liquid
> via matrixPorosity (could also depend on slotInitialAperture). Just an
> idea anyway. The problem of fluid flow inside an advancing crack tip
> seems fundamentally tough... I can't imagine a perfect answer yet.
>
> Bruno
>
> --
> You received this bug notification because you are a member of Yade
> developers, which is subscribed to Yade.
> https://bugs.launchpa...

[Luc Scholtès (luc)]

1) Remark about the aperture of the cracks:

"However, an obvious problem of this is that it result in a reference equal to zero when the crack occurs, and we need to divide by this volume..."

Currently, if the crack occurs under tensile loading, its aperture is not null when it appears. It is actually set equal to the normal relative displacement needed to break the bond (interparticle distance at the moment of breakage). Your remark made me think about the physical relevance of such model and I tend to believe that this is OK if you imagine what is happening between two stretched out cemented grains.

2) Use of "aperture" in calculation of fluid volume:

The "aperture" used in the calculation of the local conductivity is equal to:
- mechanical aperture + residual aperture if mechanical aperture > residual aperture
- residual aperture if mechanical aperture <= residual aperture
Moreover, a facet is not tricked if mechanical aperture <= 0 and residual aperture <= 0 (e.g. for closed induced cracks).
"aperture" should thus never be equal to zero for a tricked facet and it can be used for computing the volume of a cell associated to a tricked facet.

3) Computation of the crack volume:

a) Currently, the local conductivity is computed considering that the flow occurs between parallel plates with same width and length equal to the distance between the connected cells centers. If this is still OK for everyone, is it compatible with Robert's proposal? I feel like there would be something wrong in the case of a facet with multiple broken edges since the conductivity is computed considering only one of the broken edge aperture while the crack area (and therefore the crack volume) is computed as the sum of the broken edges surface... Should we add (or average?) the crack dimensions for both the conductivity and the crack area? Why don't we use directly the size of the crack available in JCFPM as pi*min(R1,R2)^2 with R1 and R2 the particles radii?

b) We are interested in the volume of fluid contained in a cell which might have n facets associated to a crack and 4-n facets not associated to a crack. How do we define the fluid volume in such a case? V = (4-n)/4*Vcell + SumOveri_0TOn Vcracki? Or do we simply consider that Vcracki >> V and do what Robert is suggesting: V = Vcell + SumOveri_0TOn Vcracki?

lUC

[Bruno Chareyre (bruno-chareyre)]

>I am particularly confused by [3]

A local crack is a polygonal face of the Laguerre graph. I did not find online content showing a Laguerre graph in 3D but you can find inspiration by looking at a Voronoi graph [1] (some polygonal faces are visible on the outer contour).
The code at [3] is looping on the edges of the polygons, which means algorithmicaly a loop on all cells incident to a particular edge of the triangulation (each triangulation cell is the dual of a Laguerre vertex).
A wise definition of the "width" of the crack should account for this geometry more accurately (currently the width is just the distance between to points of the polygon).
It would also lead more naturally to an increase of the conductivity for multiple cracks in one cell.
Your code Robert is doing something along this line (you take from the polygon a triangular subdomain), however you apparently neglected the fact that there are not just two cells associated to a cracked edge, there are many of them. So you need a loop somewhere.

[1] http://math.lbl.gov/voro++/about.html

>It's as if we are artificially creating volume space that is not reflected
by the packing geometry. Which is kind of what this all boils down to,
right?

I don't think so. Instead the algorithm is relying strongly on the idea that the volume changes are exactly reflected by packing geometry. This being the case we know the cracked volume without a need to distinguish the number of cracks and their openings. Calculating volume changes is not an issue at all, as mentioned in the bug report. The only problem is for the reference volume used in the compressibility equation (where "reference" means "before any crack").It is enough to define it once since later it can always be updated by using the volume change of the cell (hence no need to consider sums over n or 4-n cracked edges).

You guys have been progressively diverging from the bug, and now I'm following you. ;)
The questions on tricking permeability are completely independent, calculating volume changes as well, let's leave the bug alone and continue in a different thread. To delineate the bug, the question here is:
"what the initial amount of fluid before any cracks should be?"

B

[Luc Scholtès (luc)]

You are right, we probably mixed up a few things... Sorry for that.

If the main question concerns the reference volume Vo, here is my understanding of your proposals:

1 - Vo=matrixPorosity*cellVolume for a cell in the intact matrix
2 - Vo=aperture*area+matrixPorosity*cellVolume for a cell cut by a pre-existing fracture (DFN).

In the latter case, jointInitialAperture and slotInitialAperture (or, more generally, residual aperture) are the right candidates for aperture (as suggested in your first message). There should not be any problem with aperture being equal to zero because of the +matrixPorosity*cellVolume.

Q1 - How do we deal with multiple prefractured edges? DFNs usually correspond to individual fracture planes which would respectively cross 1 or 2 edges of the facet. Also, several fracture planes can eventually cross the same facet with different orientations. If we forget the multiple planes configuration for the moment, aperture should be equal to 1*residualAperture even if 2 edges are concerned (no additivity), right? If you agree on this, then, we need to think about the area -> Q2.

Q2 - What about the area? Right now, conductivity is computed assuming squared parallel plates with reference length equal to the distance between the connected cells centers. Related to Q1, area should be "weighted" by the number of fractured edges, right? Using triangular surfaces joining the middle of the edge and the cell centers would probably be ideal with respect to the geometry of our problem but then the cubic law used for computing the conductivity would be kind of flawed, right? What about associating rectangular plates of dimensions 0.5*cell1CenterToCell2Center*middleOfBrokenEdge1ToMiddleOfBrokenEdge2 to each edges? The conductivity would then be adapted correspondingly to the number of fractured edges.

Q3 - Should we consider cellVolume as V(tetrahedron) or as V(tetrahedron)-V(spheres)? I guess it depends on how we compute dV=V-Vo. If V=V(tetrahedron) then cellVolume=V(tetrahedron), right?

Q4 - How do we define matrixPorosity? Any chance that we could relate matrixPorosity to the permeabilityFactor used to scale down the overall permeability or do we simply relate to the porosity of the material? This is probably stupid to relate permeability to porosity but I am asking because it is sometimes difficult to have the data.

Q5 - If statement 2 is confirmed, do we need to take into account the induced cracks contribution to V in the calculation of dV=V-Vo?

Luc

[Bruno Chareyre (bruno-chareyre)]

Q1 - I think it is possible to define a relevant geometry for multiple cracks; it is just the intersection of one tetrahedron with multiple polygonal faces. I think the problem of multiple fracture planes crossing one edge is very special, I would escape it on the basis that you can't break another time something that is already broken, so 2 planes = 1 crack locally if the two planes intersect.

> Q2. [...] area should be "weighted" by the number of fractured edges

As explained before it should be enough to accumulate the conductivities from different cracks. In pseudo code:
for e in crackedEdges:
    for facet incident to e:
        facet->conductivity += trickPermeability(edge) #currently the '+=' is a '='
The triangular surface you mention also appeared in Robert's proposal. It makes sense to use it in trickPermeability().
Problems of your equation are that 1/ it assumes the four points are in the same plane, not the case in practice, they define two triangles in two different planes and 2/ it cannot be used for 3 broken edges.
If 2 triangles are coplanar the above loops will produce your equation for 2 broken edges (it implies that the broken edges are parallel, which can only happen at a boundary of the problem).

Q3: I would say yes, but I don't really see the problem.
Except that V(tetrahedron)-V(spheres) doesn't make much sense when the spheres are just computational nodes.

Q4. How do we define matrixPorosity?
I would leave that to the user.

Q5. No. Why would we do that? The total volume of a cell is changing and it must be balanced by some fluid fluxes regardless of the cracks.

Bruno

[Luc Scholtès (luc)]

I believe we are somehow converging.

Just a couple of additional remarks.

1 - Regarding the area associated to the broken edges:

If you look at the attached drawing, you'll see that we have actually 2 cases:
1) the 2 broken edges are associated to 2 individual cracks
2) the 2 broken edges are associated to the same fracture plane
In 1), I agree with you and the triangular surfaces are not coplanar
In 2), the triangular surfaces should be coplanar to account for the fact the the broken edges belong to the same fracture plane.
If we are concerned with the calculation of Vo, case 2) should prevail since we are dealing with pre-existing fracture planes rather than individual cracks but it would probably be easier (and more consistent?) to use the same logic for fracture planes and cracks so we should probably go for case 1). What do you guys think about that?

2 - Regarding the cubic law for computing the conductivity:

The cubic law is based on the concept of rectangular parallel plates of length L and width W (cf. drawing) and I don't know how we could adapt it to the case of triangular plates. Considering rectangular plates instead of triangular ones seems like something reasonable in our case. If you agree on this assumption, then we need to decide how to compute W. As a first guess, I would go for W equal to the distance between the middle of the broken edge (location of the crack) and the barycentre of the facet. What do you guys think about that?

Luc

[Bruno Chareyre (bruno-chareyre)]

Hi Luc,
Please keep this discussion for the bug only.
Bruno

[Janek Kozicki (cosurgi)]

moved form https://bugs.launchpad.net/yade/+bug/1734653 to https://gitlab.com/yade-dev/trunk/issues/38