Recommended approach for temperature modeling in fracture propagation with 2D DFN elements embedded in 3D rock

[tanaka792]

I would like to ask a question. I have been working on fracture propagation analysis using 2D DFN elements embedded in a 3D rock matrix, following the procedure described below, but I have reached a dead end and am unsure how temperature boundaries and mesh handling should be treated.

My goal is to perform fracture propagation analysis associated with hydraulic fracturing around a wellbore. For this purpose, I am solving heat transfer (HT) using MOOSE and coupling it with an FDEM-based mechanical analysis. However, I am struggling with how to handle temperature analysis in the fracture region. My current workflow is to solve HT in MOOSE for one second, then perform the mechanical calculation with FDEM, identify regions where fractures occur, treat those regions as blocks consisting of lower-dimensional elements (TRI3 embedded in TET4), and then repeat the MOOSE calculation.

First, I tried using a single temperature variable and assigning material properties multiplied by the fracture aperture only in the fracture region. However, as fractures propagate, the solution no longer converges. In addition, heat transfer from fractures to the matrix can only be represented by adjusting the magnitude of thermal diffusion, rather than by an explicit heat exchange term.

I then switched to using two temperature variables. However, in the PorousFlow module, especially under MPI parallelization, it is not possible to assign a temperature variable Tf only to the fracture region, and the variable must instead be defined over the entire domain. Even when a NullKernel is applied on the matrix side, nonphysical temperature increases appear in both the fracture and matrix regions. A similar issue occurs when solving a variable only on DFN elements embedded in the matrix, where the nonphysical behavior appears in regions where the NullKernel is applied, although I have not been able to fully determine the cause.

What I want to represent is a situation in which water flows at low temperature inside fractures that are directly connected to a well, and the temperature rises quickly due to heat exchange. I tried limiting the temperature range using Bounds, but the boundary conditions could not be set properly. When using PorousFlowEnthalpySink, temperature changes occur in an unnatural way in regions without fractures. On the other hand, even when fixing the temperature using Dirichlet boundary conditions, I was not able to observe the temperature of the water being advected through the fracture when heat exchange and diffusion were turned off and only advection was considered.

On the other hand, with an approach in which only the fracture mesh is extracted and handled using a multi-app setup, the fractures appear as discrete meshes, which is likely to cause convergence problems. In some cases, the fracture consists of only a single TRI3 element. Fractures do not necessarily originate from the well location, so boundary conditions cannot be applied at the well. Even when using point injection, problems arise in that it is unclear which point should be selected and how much should be injected at that point, and such settings cannot be properly specified.

Given these issues, I would like to ask if there is any approach that could represent cold fluid flowing inside fractures.

[tanaka792]

With PorousFlowEnthalpySink, fluid flow inside fractures can be represented as shown below; however, before fractures develop sufficiently, a temperature distribution like that shown in the second figure appears.

This is part of the input file configuration.

[Variables]
 [pwater]
 []
 [temperature]
 []
 [temp_f]
 []
[]

[Kernels]
 [mass_water_dot]
   type = PorousFlowMassTimeDerivative
   fluid_component = 0
   variable = pwater
 []
 [flux_water]
   type = PorousFlowAdvectiveFlux
   fluid_component = 0
   use_displaced_mesh = false
   variable = pwater
 []
 [energy_dot2]
   type = PorousFlowEnergyTimeDerivative
   variable = temperature
   block = box
 []
  [heat_conduction2]
  type = PorousFlowHeatConduction
   variable = temperature
       block = box
  []
  [Q_frac2]
   type = PorousFlowHeatMassTransfer
   block = 2  
     variable = temperature     
     v = temp_f
     transfer_coefficient = 300         
     save_in = joules_per_f
  [] 
 [no]
   type = NullKernel
   variable = temp_f
     block = 'box'
 []
 [no]
   type = NullKernel
   variable = temperature
     block = 2
 []
 [energy_dot]
   type = PorousFlowEnergyTimeDerivative
   variable = temp_f
   block = 2
 []
 [convection]
   type = PorousFlowFullySaturatedUpwindHeatAdvection
   variable = temp_f
   block = 2
   fallback_scheme = harmonic
   full_upwind_threshold = 10000   
   []
  [heat_conduction]
  type = PorousFlowHeatConduction
  variable = temp_f
  block = 2      
   []
 [Q_frac]
  type = PorousFlowHeatMassTransfer
  block = 2  
  variable = temp_f          # u = T_f
  v = temperature
  transfer_coefficient = heattransfer_rate  #300* aparture        
  save_in = joules_per_s
 [] 
[]

[AuxVariables]
 [heattransfer_rate]
    order = CONSTANT
   family = MONOMIAL
   block = 2
 []
[]

[Functions]
 [p_dirich_fun]
   type = ParsedFunction
   expression = '24.76e6+4*1e6*t'
 []
 [inj_flux]
   type = PiecewiseLinear
   x = '0   10   30   60'
   y = '0   -0.008  -0.008  -0.008'
 []
[]

[BCs]
#  [inlet_p]
#    type = FunctionDirichletBC
#    boundary = 'inlet inlet_upper inlet_lower'
#    variable = pwater
#    function = p_dirich_fun
# []
 [left_p]
   type = PorousFlowSink
   variable = pwater
   boundary = inlet
  flux_function = inj_flux
   fluid_phase = 0
[]
 [inlet_T]
   type = ADConvectionHeatTransferBC
   boundary = 'inlet' 
   variable = temperature
   T_ambient = 293
   htc_ambient = 60  #3/0.05
 []
 [inlet_T2]
   type = DirichletBC
   boundary = noflux
   variable = temperature
   value = 573
 []
 [left_T]
   type = PorousFlowEnthalpySink
   variable = temp_f
   boundary = inlet
   T_in = 293
   fp = water
   flux_function = inj_flux
   porepressure_var = pwater
 []
 [inlet_T3]
   type = DirichletBC
   boundary = noflux
   variable = temp_f
   value = 573
[]
#  [inlet_f]
#   type = DirichletBC
#   boundary = 'inlet inlet_upper inlet_lower'
#   variable = temp_f
#   value = 293
# []
[]

[GiudGiud]

Figure 2 results are unconverged right?

@lynnmunday @rpodgorney

[tanaka792]

Thanks for the reply.

While the pore pressure and matrix temperature converged well (1e-6 and 1e-12), the residual of the fracture temperature remained　constant at about 265.428.
This may be due to the ConstantBounds constraint, or because the fracture-adjacent regions use NullKernel for this variable, meaning　that no active equation was solved and the solver proceeded to the next time step.

Do you have any recommendations for resolving this issue?

[rpodgorney]

This is doable. We do it on a regular basis for geothermal programs. I don’t have the time to provide detail I’ve about approach right now. The crux of the thing is you’re sharing nodes between the fracture and the matrix so essentially have one temperature variable, and the fracture and matrix have to be in thermal equilibrium. This can cause some retardation and heat transport in the fracture because too much heat is lost to the matrix.

If you want a separate temperature variable, you need to mesh the inside of each fracture as a 3-D volume not a 2-D plane and then come up with a transfer coefficient between the two. This is a non-trivial process. Another approach is a solve something like a lubrication equation inside the fracture we’ve done that too with limited success.

[rpodgorney]

Look through the examples inFALCON and porous flow. I’m pretty certain there is something in there

[tanaka792]

Thank you for your response.
As you mentioned, representing fractures as 3-D volumes is difficult, and from my perspective it　would also require remeshing every time fractures propagate, so I would prefer to avoid that approach.
As you pointed out, it was helpful to confirm that using a single temperature variable indeed has many limitation. I was not familiar with the lubrication equation, but I will look into it.
In parallel, I would like to continue thinking about whether it is possible to separate the temperature into two variables

[tanaka792]

I have been looking at the FALCON examples for some time now (e.g., example01), but so far I have not found an example that I can directly apply to my current work. Most of the examples focus on connected fracture meshes and point injection.

[lynnmunday]

Are you doing this coupling in MOOSE using the multiapp system? Or do you manually run each simulation separately with postprocessing steps to create the discrete fractures from the DEM simulation. If it is the latter, maybe you can put in some graph algorithms to make sure your fractures are fully connected to one another and the well. But you'll still get the retarded heat transfer Rob mentions.

I think the method outlined in this paper (https://inldigitallibrary.inl.gov/sites/sti/sti/Sort_88359.pdf) is how you would need to model this to the resolution you are wanting. I think this is the lubrication method Rob is referring to. Its a hard simulation to run.

It looks like a high density of fractures, would it make more sense to homogenize them when you transfer them to the porous flow simulation. The porousflow simulation would no longer use discrete 2D fractures, instead the fractures would be represented by changing the permeability of a volume element. This should be a lower bound on what you are really trying to do.

[tanaka792]

Thank you very much for your valuable feedback.

While we employ the MultiApp system for hydrothermal analysis simulations using MOOSE, the generation of discrete fracture structures involves transferring relevant data to FDEM for separate processing.

The graph algorithm seems like an effective and promising approach. I believe this could help resolve issues related to defining boundary conditions, which had been challenging due to the fracture state. Ultimately, since we need to consider fracture closure, this is a technology that will likely be necessary at some point, so I want to investigate it.

At this stage, our goal is to represent heat exchange. We are exploring whether an approach exists to treat fractured zone temperature and matrix temperature as separate variables. On the other hand, applying a homogenization approach after fracturing has progressed to a certain extent would likely result in a single temperature variable. However, since this could potentially be determined based on the number of fracture elements, it is attractive from an implementation perspective. If a feasible two-variable approach cannot be found, we would like to try this option.But still, implementing lubrication theory seems difficult.

[tanaka792]

Thank you all for the valuable comments and suggestions.
I am currently considering a model with two temperature variables defined on the same mesh: one for the porous solid phase and one for the pore water, based on a two-variable formulation.
The benchmark study is not yet complete, so I cannot judge the correctness of this model at this stage.

With this formulation, heat transfer into the rock matrix can be expressed through the temperature difference between the rock and the injected cold water, while the inflow of cold water into the fracture can be represented via boundary conditions. In this setup, assigning fracture-like material properties to a subset of blocks representing the water domain did not cause convergence issues.

I understand the modeling approaches suggested above. At this stage, I cannot yet assess their suitability for my case, but I will keep them in mind as references when further developing

[tanaka792]

I would like to decide on the final answer after confirming the accuracy through numerical benchmarks, but I encountered an issue while performing a heat-exchange benchmark.

1. When running the FALCON heat-exchange example with MPI parallelization (40 ranks), why does the fracture temperature change while the matrix temperature does not change at all, even though the serial (MPI = 1) result matches the documentation?
2. Why do almost all official examples use DiracKernels for heat exchange? Is this considered the most appropriate or robust approach, especially under MPI parallelization?

As a reference for heat exchange, I downloaded the FALCON example and ran it with MOOSE version 2.1. When running with MPI = 1, or without specifying MPI, the results match those shown on the documentation page. However,No runtime errors are reported, when running in MPI parallel with 40 ranks, I observed that the temperature in the fracture domain changes as expected, while the temperature in the matrix domain does not change at all. I would appreciate any insight into the possible cause of this behavior.I also tried replacing older transfer implementations with newer ones, but this did not improve the situation.

Regarding modeling approaches for heat exchange, I understand that there are several possible options, such as transferring the fracture temperature to the matrix and using PorousFlowHeatMassTransfer, or storing quantities such as joule_per_s and applying them via a body force. Based on this idea, I tried representing the fracture location as a block and applying kernels in that region. However, compared to the official example, the heat did not spread into the matrix as much.I would therefore like to ask why almost all official examples use DiracKernels for heat exchange. Is this because DiracKernels are considered the most appropriate or robust approach for this type of problem?

[lynnmunday]

What benchmark example are you running and finding these MPI inconsistencies. I need to look at it.

[lynnmunday]

I use diracKernels because I model well injection and production. I'm modeling a much coarser scale than you. Andy Wilkins pointed out several issues with BCs in porousFlow simulations in another discussion a few weeks ago that you should look at and recommended this documentation about porousFlow bcs: https://mooseframework.inl.gov/modules/porous_flow/boundaries.html

[tanaka792]

This benchmark is based on the FALCON example (https://mooseframework.inl.gov/falcon/examples/example01.html).
https://github.com/tanaka792/FALCON
The detailed modifications are described in README.md; except for the Transfer settings and output configuration, the setup is mostly unchanged from the original example.
With these modifications, the case has been confirmed to run correctly in MOOSE using mpiexec -np 1. However, when running with MPI ≥ 2, the Matrix_T results in mrect1 do not change at all.