# RootBeam
Analysis of coupled axial and lateral deformation of roots in soil

Notes 27/09/2018, GJMeijer  
Python code to analyse the behaviour of roots loaded in direct shear by solving the (large deformation) bending equation, taking i to account both axial and bending stiffness.

## Publication details
This code belong to journal publication:  

> Analysis of coupled axial and lateral deformation of roots in soil  
> GJ Meijer, D Muir Wood, JA Knappett, AG Bengough, T Liang  
> International Journal for Analytical and Numerical Methods in Geomechanics  
> 2018, volume x, issue x, page xx-xx  
> DOI: xxx

This code can be used to reproduce some of the shear tests with ABS described in this paper (Liang et al.)


## Quickstart
Ensure that all python files (`.py`) and input files (`.csv`) are placed in the same directory.
Analysis can be run by executing the `Main.py` script.
In the first lines of `Main.py`, some user-defined settings are required
* The name of the input parameter file, including extension
* Some booleans to decide whether results should be plotted, saved etc.

## Python requirements
Code is written in Python 3
Executing required installed modules (apart from standard modules):
* `SciPy`
* `NumPy`
* `MatPlotLib`


# Notes on data format

## Units
Throughout the code, units used:
* Length in millimeters (mm)
* Force in newtons (N)
* Angles in radians (rad)
Note that therefore, in constrast to common geotech practice:
* Stesses and stiffnesses are in (MPa)
* Unit weights are in (N/mm3)

## Input parameter file - required parameters
Required columns:

| Parameter | Type | Unit | Description |
|----------------------------------------------------------|-------------|------------------------------------------------------------------------------------------------------------------------------|---------------------------------------------------------------------------------------------------------|
| `RunID` | integer | - | Run identifier (ensure these are all unique values for every row) |
| `ModelName` | string | - | Name for this analysis. Name is used to create output directory with results |
| `NodeFile` | string | - | File name with node input data, including .csv extension |
| `SegmentFile` | string | - | File name with segment input data, including .csv extension |
| `phi` | double | rad | Soil peak angle of internal friction (used to calculate lateral resistance) |
| `phi_cs` | double | rad | Soil critical state angle (used to calculate effect of soil confinement on additional shear resistance) |
| `delta` | double | rad | Interface friction angle between soil and root |
| `gamma` | double | N/mm3 | Soil effective unit weight |
| `K` | double | - | Coefficient of lateral earth pressure |
| `sigmav0` | double | MPa | Vertical effective stress at the soil surface |
| `shearplane_param` | double | mm | Shear plane thickness parameter (b_sh) in paper |
| `shearplane_depth` | double | mm | Shear plane depth, at Y=0 |
| `shearplane_angle` | double | rad | Shear plane angle |
| `n_node` | integer | - | Number of nodes per segment |
| `n_step_max` | integer | - | Maximum number of displacement steps |
| `n_iter_max` | integer | - | maximum number of iterations (to find solution for a single displacement step) |
| `uext_inc0` | double | mm | Initial shear displacement step size |
| `uext_max` | double | mm | Maximum shear displacement (analysis stops once this value is reached) |
| `uext_factor_increase` | double | - | Multiplaction factor to automatically increase displacement step size if suitable |
| `uext_factor_decrease` | double | - | Multiplaction factor to decrease displacement step if solution if solution not found |
| `uext_incmin` | double | mm | Minimum displacement step size (if lower required, analysis stops) |
| `uext_incmax` | double | mm | Maximum displacement step size |
| `solve_tolerance` | double | - | Solver tolerance of <scipy.integrate.solve_bvp> solver in Python |
| `SaveStep` | integer | - | Save full root and soil data every <SaveStep> displacement step |`

## Root geometry
The root architecture is described by two seperate files:
* Nodal data   - This defined the position and displacement boundary conditions of individual nodes (points) in the geometry
* Segment data - This defines the connections between two nodes. For each connection, segment properties such as stiffness of diameter are given
Theoretically, there is no limit to the amount of segments, although it can be expected that the code will become less efficient with many segments (as the number of differential equations to solve increases)

### Node input data format
Required columns:

| Parameter | Type | Unit | Description |
|----------------------------------------------------------|-------------|------------------------------------------------------------------------------------------------------------------------------|--------------------------------------------------------------------------------------|
| `NodeID` | integer | - | Node identifier (ensure these are all unique values) |
| `X` | double | mm | Global X-position of node |
| `Y` | double | mm | Global Y-position of node |
| `bound_X` | integer | - | X-displacement constrained (value=1) or free (0) |
| `bound_Y` | integer | - | Y-displacement constrained (value=1) or free (0) |
| `bound_Theta` | integer | - | rotation constrained (value=1) or free (0) |

### Segment input data format
Required columns:

| Parameter | Type | Unit | Description |
|----------------------------------------------------------|-------------|------------------------------------------------------------------------------------------------------------------------------|--------------------------------------------------------------------------------------|
| `SegmentID` | integer | - | Segment identifier (ensure these are all unique values) |
| `NodeID1` | double | - | Node identifier of node at start of segment |
| `NodeID2` | double | - | Node identifier of node at end of segment |
| `d` | double | mm | Root diameter |
| `Et1` | double | MPa | Coefficient in root tensile stiffness (\( \xi_1 \) in paper) |
| `Et2` | double | - | Coefficient in root tensile stiffness (\xi_2 in paper) |
| `Et3` | double | MPa | Coefficient in root tensile stiffness (\xi_3 in paper) |
| `Eb1` | double | MPa | Coefficient in root bending stiffness (\xi_1 in paper) |
| `Eb2` | double | - | Coefficient in root bending stiffness (\xi_2 in paper) |
| `Eb3` | double | MPa | Coefficient in root bending stiffness (\xi_3 in paper) |

# General remarks on how the code works
* The program loads all input data in input parameter file specified
Every line in this file represents a single analysis (i.e. a single direct shear tests)
Required parameters are given above, as are node and segment parameters for the root architecture
Subsequently, the code runs through every run in sequence
* The set of differential equations and boundary conditions is solved for every displacement step. Initial step size is `uext_inc0`
The initial guess required is based on the last solution
When no solutions can be found, the displacement step is decreased (multiplying current step times `uext_factor_decrease`) to increase the chance that the system will be solved.
* If the displacement step becomes too small (smaller than `uext_incmin`), the analysis for this run will automatically stop.
The displacement step is multiplied by `uext_factor_increase` when steady solutions are found systematically, but will never be larger than `uext_incmax`.
The analysis stops when no solution can be found within acceptable step sizes, or when the final displacement `uext_max` is reached
* The set of differential equations is written in a set of first order differential equations, as is a requirement for the `solve_bvp` solver in SciPy. All differential equations are normalised by segment length.
* The order of variables is: theta, dtheta/ds. d^2theta/ds^2, eps, u, w. When there are multiple segments, these are added in sequence of occurance in the `dsegm` segment data dictionary (so first all variables for the the first segments, than all those for the second etc.)

# Outputs
Depending on the output settings in `Main.py`, output will be generated and saved. 
When `data_save` equals `True`, the code prints output `.csv` files to a new directory `Results`. 
Subfolders will be created based on the `ModelName` and `RunID` given in the input parameter file.
Data will be saved for every `SaveStep` displacement step for both soil and root, as well as summary data for the whole analyses.

## Output - Root
For every displacement step, a `.csv` is created containing the following fields:

| Parameter | Type | Unit | Description |
|----------------------------------------------------------|-------------|------------------------------------------------------------------------------------------------------------------------------|--------------------------------------------------------------------------------------|
| `s` | double | mm | Local coordinate along the (displaced) root axis |
| `t` | double | rad | Root rotation (root coordinate system) |
| `dt` | double | rad/mm | Root curvature (dt/ds) |
| `ddt` | double | rad/mm2 | Derivative of root curvature (d^2t/ds^2) |
| `eps` | double | - | Root axial strain |
| `u` | double | mm | Root axial displacement (root coordinate system) |
| `w` | double | mm | Root lateral displacement (root coordinate system) |
| `x` | double | mm | Root x-position (root coordinate system) |
| `y` | double | mm | Root y-position (root coordinate system) |
| `Xr` | double | mm | Root X-position (global coordinate system) |
| `Yr` | double | mm | Root Y-position (global coordinate system) |
| `umob` | double | mm | Relative axial soil-root displacement (delta u in paper) |
| `wmob` | double | mm | Relative lateral soil-root displacement (delta w in paper) |
| `SegmentID` | integer | - | Segment identifier |

## Output - Soil
For every displacement step, a `.csv` is created containing the following fields:

| Parameter | Type | Unit | Description |
|----------------------------------------------------------|-------------|------------------------------------------------------------------------------------------------------------------------------|--------------------------------------------------------------------------------------|
| `Xs` | double | mm | Soil X-position (global coordinate system) |
| `Ys` | double | mm | Soil Y-position (global coordinate system) |
| `SegmentID` | integer | - | Segment identifier |

## Output - Summary
Three summary data files are outputted
* `ModelName` + `_Results_ShearForceDisp.csv`:

| Parameter | Type | Unit | Description |
|----------------------------------------------------------|-------------|------------------------------------------------------------------------------------------------------------------------------|--------------------------------------------------------------------------------------|
| `StepID` | integer | mm | Displacement step identifier |
| `uext` | double | mm | Shear displacement |
| `F_total` | double | N | Total root-reinforcement in middle of shear plane |

* `ModelName` + `_Results_ShearReinforcement.csv`:

| Parameter | Type | Unit | Description |
|----------------------------------------------------------|-------------|------------------------------------------------------------------------------------------------------------------------------|---------------------------------------------------------------------------------------------------------------|
| `StepID` | integer | mm | Displacement step identifier |
| `SegmentID` | integer | - | Segment identifier |
| `F_direct` | double | N | Direct contribution of force in root to shear resistance (force component along shear displacement direction) |
| `F_comfine` | double | N | Additional normal force on shear plane due to roots (positive force compresses soil further) |
| `F_total` | double | N | Total root-reinforcement in middle of shear plane |

* `ModelName` + `_Results_Stresses.csv`:

| Parameter | Type | Unit | Description |
|----------------------------------------------------------|-------------|------------------------------------------------------------------------------------------------------------------------------|--------------------------------------------------------------------------------------|
| `StepID` | integer | mm | Displacement step identifier |
| `SegmentID` | integer | - | Segment identifier |
| `sig_ax_shearplane` | double | MPa | Axial stress at shear plane |
| `sig_be_shearplane` | double | MPa | Bending stress at shear plane |
| `sig_sh_shearplane` | double | MPa | Max. shear stress at shear plane |
| `sig_ax_max` | double | MPa | Max axial stress in segment |
| `sig_be_max` | double | MPa | Max bending stress in segment |
| `sig_sh_max` | double | MPa | Max shear stress in segment |
| `sig_fi_max` | double | MPa | Max stress in ultimate fibre in segment |
| `sig_fi_max_ax` | double | MPa | Axial stress component of `sig_fi_max` |
| `sig_fi_max_be` | double | MPa | Bending stress component of `sig_fi_max` |
| `X_sig_ax_max` | double | mm | Global (displaced) X-position where `sig_ax_max` occurs |
| `X_sig_be_max` | double | mm | Global (displaced) X-position where `sig_be_max` occurs |
| `X_sig_sh_max` | double | mm | Global (displaced) X-position where `sig_sh_max` occurs |
| `X_sig_fi_max` | double | mm | Global (displaced) X-position where `sig_fi_max` occurs |