Purefrac

Purefrac generates 3D discrete fracture networks (DFN) without matrix porosity.

Overview

This MATLAB program generates 3D discrete fracture networks (DFN) without matrix porosity. It's designed for simulating fractured rock systems where fractures are the dominant flow paths.

Key Features

1. Multiple Fracture Sets

• Up to 3 independent fracture sets
• Each set with unique orientation and properties
• Realistic geological patterns

2. Fracture Properties

• Length and width (rectangular planes)
• Aperture (thickness)
• Dip angle (0-90°)
• Strike angle (0-360°)
• 3D rotation matrices for accurate orientation

3. Advanced Analysis

• Topology: Vertices, edges, faces, cells
• Connectivity: Connected component analysis
• Percolation: Tests in X, Y, Z directions
• P32 Intensity: Surface area per unit volume
• Fracture spacing: Statistical analysis
• Skeletonization: Network backbone extraction

4. Geological Visualizations

• 3D fracture network rendering
• Stereographic projection (lower hemisphere)
• Rose diagram of strike directions
• Cross-sectional views
• Aperture distribution histograms

5. Export Capabilities

• STL files (full network and individual sets)
• CSV data (for external analysis)
• Detailed text reports
• MATLAB workspace

Predefined Scenarios

1. Standard DFN (standard_dfn)

Application: General fractured rock

• Sub-horizontal bedding (8 fractures)
• NE-SW steep joints (10 fractures)
• NW-SE steep joints (10 fractures)

2. Dense Network (dense_network)

Application: Highly deformed or weathered rock

• 45 total fractures
• Higher density, smaller apertures
• Complex connectivity

3. Sparse Network (sparse_network)

Application: Competent crystalline rock

• 14 total fractures
• Larger sizes, wider apertures
• Simple connectivity

4. Bedding Planes (bedding_planes)

Application: Sedimentary sequences

• 12 horizontal bedding planes
• 8 vertical joints (minor)
• Anisotropic structure

5. Conjugate Faults (conjugate_faults)

Application: Tectonic fault systems

• Two intersecting steep fault sets
• Large fractures with wide apertures
• Characteristic 60° conjugate angle

6. Single Set (single_set)

Application: Columnar jointing, testing

• 20 parallel vertical fractures
• Simple geometry for validation

Installation & Usage

Quick Start (Standalone)

% Run the complete program with default settings
pure_fracture_network_generation 

Recommended Workflow (Config-based)

% Step 1: Configure parametersconfigure_pure_fracture_network 
% Step 2: Modify scenario in the config script (optional)% Open configure_pure_fracture_network.m% Change: scenario = 'dense_network';
% Step 3: Generate networkpure_fracture_network_generation_from_config 

Parameter Reference

Domain Parameters

config.nx = 100;  % Grid size in X
config.ny = 100;  % Grid size in Y
config.nz = 100;  % Grid size in Z
config.randomSeed = 42;  % For reproducibility

Fracture Set Parameters

config.set1.number = 10;          % Number of fractures
config.set1.minLength = 40;       % Min fracture length (voxels)
config.set1.maxLength = 80;       % Max fracture length (voxels)
config.set1.minWidth = 40;        % Min fracture width (voxels)
config.set1.maxWidth = 80;        % Max fracture width (voxels)
config.set1.aperture = 2;         % Fracture thickness (voxels)
config.set1.dipAngle = [0, 15];   % Dip range [min, max] (degrees)
config.set1.strikeAngle = [0, 360]; % Strike range [min, max] (degrees)
config.set1.color = [1, 0, 0];    % RGB color for visualization
config.set1.name = 'Set Name';    % Descriptive name

Orientation Conventions

• Dip: Angle from horizontal (0° = horizontal, 90° = vertical)
• Strike: Azimuth direction (0° = North, 90° = East, 180° = South, 270° = West)
• Right-hand rule: Strike is perpendicular to dip direction

Output Files

1. STL Files

• pure_fracture_network.stl - Complete network
• pure_fracture_network_solid.stl - Solid matrix
• fracture_set_1.stl, fracture_set_2.stl, fracture_set_3.stl - Individual sets

2. Data Files

• fracture_network_data.txt - Comprehensive analysis report
• fracture_data.csv - Tabular data for Excel/Python
• pure_fracture_network_workspace.mat - Complete MATLAB workspace

3. Analysis Report Contents

• Network statistics (density, surface area, P32)
• Connectivity analysis
• Percolation status (X, Y, Z directions)
• Individual fracture properties (center, size, orientation)

Visualization Guide

Figure 1: Complete Fracture Network

3D isosurface of the entire fracture network with realistic lighting.

Figure 2: Network with Skeleton

Shows the fracture network centerline (skeleton) overlaid on transparent fractures. Useful for understanding connectivity topology.

Figure 3: Fracture Centers by Set

Color-coded markers showing fracture center locations. Each set has a different color for easy identification.

Figure 4: Aperture Distribution

Histogram showing the statistical distribution of fracture apertures throughout the network.

Figure 5: Stereographic Projection

Lower hemisphere projection of fracture poles (normals). Standard structural geology representation for analyzing orientation patterns.

Figure 6: Rose Diagram

Circular histogram showing the distribution of fracture strike directions. Helps identify preferred orientations.

Figure 7: Cross-sections

Three orthogonal slices (XY, XZ, YZ) through the domain showing fracture distribution patterns.

Advanced Analysis

P32 Fracture Intensity

P32 = Total Fracture Surface Area / Domain Volume

• Typical values: 0.001 - 0.1 [1/length]
• Higher P32 = more fractured rock
• Critical parameter for flow modeling

Percolation Analysis

Tests if connected fracture pathways exist across the domain:

• Yes: Flow can traverse the entire domain
• No: Isolated fracture clusters only

Connectivity Components

• Reports number of disconnected fracture clusters
• Largest component indicates main flow pathway
• Small isolated clusters may represent dead-end porosity

Applications

1. Geothermal Energy

• Model fractured geothermal reservoirs
• Analyze flow pathways and heat exchange
• Estimate permeability tensors

2. Oil & Gas

• Characterize naturally fractured reservoirs
• Hydraulic fracture network simulation
• Production optimization

3. Groundwater Hydrology

• Contaminant transport in fractured aquifers
• Fracture-dominated flow systems
• Well capture zone analysis

4. Nuclear Waste Storage

• Repository safety assessment
• Radionuclide transport modeling
• Long-term barrier integrity

5. CO2 Sequestration

• Storage in fractured formations
• Leakage pathway identification
• Caprock integrity assessment

6. Computational Fluid Dynamics

• Export STL for mesh generation
• Direct flow simulation
• Permeability calculation

Customization Tips

Creating Realistic Networks

Sedimentary Rocks:

- Horizontal bedding: dipAngle = [0, 10] 
- Sparse vertical joints: 4-6 fractures 
- Wide bedding planes: aperture = 2-3

Crystalline Rocks:

- Steep joints: dipAngle = [70, 90] 
- Multiple intersecting sets 
- Smaller apertures: aperture = 1-2

Fault Zones:

- Large fractures: length = 60-90
- Wide apertures: aperture = 3-5
- Fewer fractures but bigger: number = 4-8Adjusting Fracture Density

Increase Density:

• Increase number for each set
• Reduce minLength and maxLength
• Use smaller domain (faster generation)

Decrease Density:

• Reduce number for each set
• Increase minLength and maxLength
• Use larger domain

Optimization Tips

1. Use config-based workflow for multiple runs
2. Start small (50³) to test parameters
3. Reduce visualization for large domains (comment out figure commands)
4. Use parallel computing for parameter studies (future implementation)

Troubleshooting

Issue: No Percolation Detected

Causes:

• Too few fractures
• Fractures too small
• Isolated sets don't intersect

Solutions:

• Increase number in each set
• Increase maxLength and maxWidth
• Ensure sets have intersecting orientations

Issue: Memory Error

Solutions:

• Reduce domain size (nx, ny, nz)
• Close other MATLAB figures
• Reduce number of fractures
• Use 64-bit MATLAB

Issue: Fractures Not Visible in Visualization

Solutions:

• Check aperture is at least 1-2 voxels
• Increase fracture length and width
• Check orientation ranges aren't too narrow
• Verify number > 0 for each set

Issue: Generation Takes Too Long

Solutions:

• Reduce domain size
• Reduce total number of fractures
• Comment out complex visualizations (stereographic projection)

Validation & Quality Control

Check Network Quality

1. Fracture Density: Should be 0.001-0.1 (0.1%-10%)
2. Percolation: At least one direction should percolate
3. Connectivity: Largest component should be >50% of total
4. P32 Value: Compare with field data (typically 0.01-0.5)

Visual Inspection

1. Check cross-sections for uniform distribution
2. Verify fracture orientations in stereographic plot
3. Ensure no clustering artifacts
4. Confirm aperture distribution is reasonable

Citation

If you use this code in your research, please cite:

Jimmy X. Li (2025). Purefrac (https://www.mathworks.com/matlabcentral/fileexchange/<...>), MATLAB Central File Exchange. Retrieved October 22, 2025.

License

MIT License - See script header for complete license text

Support & Further Development

Future Enhancements

• Stochastic fracture length distributions (power law)
• Fracture clustering and spacing models
• Mechanical aperture vs hydraulic aperture
• Fracture intersection analysis
• Direct permeability tensor calculation
• Export to external flow simulators (TOUGH2, FEHM, etc.)

Related Tools

• Export to OpenFOAM: Use STL files for CFD mesh
• Python integration: Use CSV files with pandas/numpy
• Geostatistical analysis: Import into GeoStats.jl or SGeMS