functionGridGenerator

Information

1. Purpose of model

Pure function that generates optimized discretization grid parameters from physical aquifer properties and desired characteristic grid spacing. The function implements grid optimization strategies from Gillner, Jin, and Speerforck (2025) using logistic functions for radial discretization and quadratic polynomials for vertical discretization. Automatically computes the number of radial control volumes based on the desired minimum grid spacing.

2. Level of detail, physical effects considered, and physical insight

L1: Automatic grid generation based on physical parameters for a given number of radial control volumes.

• The number of radial control volumes N_r is provided as an input. It is computed beforehand by GridGenerator_SizeOnly from the desired first cell width dx_min; this function only builds the grid for that N_r.
• Radial discretization uses a logistic (S-shaped) function to concentrate fine control volumes near the well (high temperature gradients) and coarsen them away from the well (low gradients). The normalized logistic widths (sum = 1) are scaled by the domain radius R_domain. Output arrays are sized to N_r values only (no padding).
• Thermal radius R_th (equation 10, Gillner et al. 2025) determines the horizontal extent of the thermal plume and scales the radial domain size R_domain = domain_factor * R_th.
• Vertical discretization in the aquifer uses a quadratic polynomial to concentrate fine control volumes at top and bottom boundaries (heat conduction to adjacent layers) and coarsen them in the middle.
• Aquitard discretization is uniform, as vertical gradients in confining layers are secondary.

3. Limits of validity

- Grid generation parameters (coefficients b, A, k for logistic; b0, b1, b2 for polynomial) are optimized for the thermal conditions described in Gillner et al. (2025) and may require adjustment for significantly different operating conditions.

- Domain radius is computed as domain_factor * R_th; default domain_factor = 3.0 provides a balance between domain size and computational cost. Increase if boundary effects are observed.

- N_r is an input constrained to [5, 500]. It should be obtained from GridGenerator_SizeOnly so that the first cell width matches the desired dx_min; remaining cells follow the logistic distribution scaled to the domain radius.

- Function assumes positive volumetric heat capacities and aquifer dimensions. Thermal radius R_th must be positive.

- First cell must have positive width; if dx_raw[1] (logistic-derived) would result in zero or negative width after scaling, the result becomes ill-defined. This should not occur for typical input ranges.

4. Interfaces

Inputs:

c_f — volumetric heat capacity of fluid [J/(m³·K)]

c_a — effective volumetric heat capacity of aquifer [J/(m³·K)]

H_a — thickness of aquifer [m]

H_c — thickness of aquitard (upper and lower) [m]

V_inj — injected thermal energy storage volume [m³]

domain_factor — multiplicative factor for R_th to determine domain extent [-] (default 3.0)

N_r — number of radial control volumes [-], range [5, 500] (computed beforehand by GridGenerator_SizeOnly)

N_z — number of vertical control volumes in aquifer [-]

N_z_c — number of vertical control volumes in aquitard [-]

Outputs (for use in Discretization function):

R_GS_H[:] — number of control volumes per radial grid section (all = 1, size = N_r)

dx_GS[:] — radial thickness of control volumes [m] (size = N_r, first cell ≈ dx_min used in GridGenerator_SizeOnly)

H_GS_VA[:] — number of control volumes per vertical aquifer section (all = 1, size = N_z)

dz_GS_A[:] — height of control volumes in aquifer [m] (size = N_z)

H_GS_VC[:] — number of control volumes per vertical aquitard section (all = 1, size = N_z_c)

dz_GS_C[:] — height of control volumes in aquitard [m] (size = N_z_c)

R_th — thermal radius of the storage plume [m]

5. Nomenclature

(See equation definitions below for parameter meanings.)

6. Governing Equations

Thermal radius and domain extent:

R_th = sqrt( (c_f / c_a) * (V_inj / (pi * H_a)) )

R_domain = domain_factor * R_th

(The number of radial control volumes N_r is an input; it is determined beforehand by GridGenerator_SizeOnly so that dx[1] ≈ dx_min.)

Radial discretization via logistic function (equation 13):

f(r) = b + A / (1 + exp(-k * (r_norm - 1)))

where r_norm = (i - 1) / N_r for i = 1..N_r+1

Normalized cell widths (sum = 1, from helper function LogisticGrid_): dx_raw[i] = (cumsum[i+1] - cumsum[i]) / cumsum[N_r+1]

Scaled cell width: dx[i] = dx_raw[i] * R_domain

Result: cells follow the logistic distribution, sum(dx) = R_domain, first cell ≈ dx_min (by choice of N_r)

Optimized coefficients (Table 2, Gillner et al. 2025):

b = 0.5095, A = 31.82, k = 23.10

Vertical discretization in aquifer via polynomial (equation 12):

f(z) = b2 * z_norm^2 + b1 * z_norm + b0

where z_norm = (i - 1) / N_z for i = 1..N_z

Cell height: dz[i] = (f(z[i]) / sum(f)) * H_a

Optimized coefficients (Table 2, Gillner et al. 2025):

b0 = 1.1, b1 = 0.75, b2 = -1.0

Vertical discretization in aquitard:

Uniform cell height: dz_c[i] = H_c / N_z_c

7. Remarks for Usage

Call this function to generate grid parameters from physical aquifer properties for a given number of radial control volumes N_r. It eliminates the need for manual specification of grid sections (R_GS_H, dx_GS) by generating optimized cell distributions. The matching N_r for a desired first cell width dx_min must be obtained beforehand from GridGenerator_SizeOnly (e.g. via the script AutoConfigureN_r.mos).

Larger N_r (finer grids) result in longer computation times; smaller N_r reduces computational cost but may reduce accuracy near the thermal plume center.

Output arrays have variable sizes: R_GS_H and dx_GS have length N_r; H_GS_VA and dz_GS_A have length N_z; H_GS_VC and dz_GS_C have length N_z_c. Arrays can be passed directly to Discretization without slicing.

Typical call pattern in Geometry.mo:

initial equation

  (R_GS_H_opt, dx_GS_opt, H_GS_VA_opt, dz_GS_A_opt, H_GS_VC_opt, dz_GS_C_opt) =

    Base.Functions.GridGenerator(C_w, C, setting.H_a, setting.H_c, setting.V_inj, setting.r_0, setting.domain_factor, setting.N_r, setting.N_z, setting.N_z_c);

8. Validation

no validation

9. References

no remarks

10. Version History

Model created by Markus Gillner (markus.gillner@tuhh.de) on 19.05.2026

Model revised by Markus Gillner (markus.gillner@tuhh.de), June 2026

- Updated documentation to reflect that N_r is now an input (computed by GridGenerator_SizeOnly) instead of being calculated internally from dx_min.

Inputs

Outputs