You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
Copy file name to clipboardExpand all lines: cookbooks/CPO_induced_anisotropic_viscosity/doc/CPO_induced_anisotropic_viscosity.md
+20-17Lines changed: 20 additions & 17 deletions
Display the source diff
Display the rich diff
Original file line number
Diff line number
Diff line change
@@ -7,33 +7,33 @@ This cookbook explains how to use the CPO-induced anisotropic viscosity material
7
7
8
8
## Introduction
9
9
10
-
Individual crystals of the mineral olivine reorganized their orientations into the crystal-preferred orientations (CPO) under deformation. The directional-dependence (anisotropy) in the viscous properties of olivine crystals suggests that the effective viscosity for olivine rocks/aggregates is different when deformations occur in different directions relative to the CPO. This CPO-induced anisotropic viscosity material model computes anisotropic viscosity based on the CPO evolution predicted by D-Rex ({cite}`fraters_billen_2021_cpo`; {cite}`kaminski2004`) and modifies the subsequent deformation.
10
+
Individual crystals of the mineral olivine reorganize their orientations into the crystal-preferred orientations (CPO) under deformation. The viscous properties of olivine crystals are direction-dependent (anisotropic), which suggests that the effective viscosity for olivine rocks/aggregates is different when deformations occur in different directions relative to the CPO. This cookbook model computes an anisotropic viscosity based on the CPO evolution predicted by D-Rex ({cite}`fraters_billen_2021_cpo`; {cite}`kaminski2004`) and includes this information in the subsequent modeling process.
11
11
12
12
Our constitutive equation for the relationship between the strain rate and stress using the anisotropic viscosity tensor is adapted from {cite:t}`signorelli:etal:2021`:
$\gamma_0=1.1\times 10^{5}$ is the isotropic fluidity, $Q=530$ $kJ/mol$ is the activation energy, $R=8.314 m^3 \cdot Pa \cdot K^{−1} \cdot$ $mol^{−1}$ is the gas constant, $d=0.001$ $m$ is the grain size, and $m=0.73$ is the grain size exponent. These values for olivine are taken from rock experiments performed by {cite:t}`hansen:etal:2016` and {cite:t}`HK04`. $J(\sigma_{ij})$ is the equivalent yield stress, where $\sigma_{ij}$ is the deviatoric (anisotropic) stress computed using the tensorial and scalar component of the anisotropic viscosity:
and $A_{ij}$ is the anisotropic tensor of fluidity in Kelvin notation:
34
34
35
35
```{math}
36
-
:label: eqn:four
36
+
:label: eqn:anisotropic_fluidity
37
37
A_{ij}=\frac{2}{3} \left[
38
38
\begin{matrix}
39
39
F+H & -F & -H & 0 & 0 & 0 \\
@@ -47,17 +47,17 @@ F+H & -F & -H & 0 & 0 & 0 \\
47
47
48
48
$J(\sigma_{ij})$ and $A_{ij}$ are computed using Hill coefficients $H, J, K, L, M,$ and $N$ {cite}`hill:1948`, which are obtained from regression analysis of a texture database constructed with olivine textures from laboratory experiments, shear box models, and subduction models (Kiraly et al., in rev.).
49
49
50
-
In this material model plugin, strain rate, density, temperature, and other parameters are taken as input to compute the viscosity, which is passed into the Stokes system to compute the stress. As a result, we adapt {math:numref}`eqn:one` to be:
50
+
In this material model plugin, strain rate, density, temperature, and other parameters are taken as input to compute the anisotropic viscosity, which is passed into the Stokes system to compute the stress. As a result, we adapt {math:numref}`eqn:anisotropic_general_stress` to be:
Since the Hill coefficients are defined in the microscopic CPO reference frame, and parameters computed in ASPECT are in the macroscopic model reference frame, several reference frame conversions are needed. First, we need to rotate $\sigma_{ij}$ in {math:numref}`eqn:three` from the model reference frame to the CPO reference frame so that $J(\sigma_{ij})$ is in the CPO reference frame. This is achieved by constructing a matrix from the eigenvectors corresponding with the largest eigenvalues of the covariance matrix for the a-, b-, and c-axis of olivine textures and then we assign the nearest orthogonal matrix to be the rotation matrix R:
57
+
Since the Hill coefficients are defined in the microscopic CPO reference frame, and parameters computed in ASPECT are in the macroscopic model reference frame, several reference frame conversions are needed. First, we need to rotate $\sigma_{ij}$ in {math:numref}`eqn:equivalent_yield_stress` from the model reference frame to the CPO reference frame so that $J(\sigma_{ij})$ is in the CPO reference frame. This is achieved by constructing a matrix from the eigenvectors corresponding with the largest eigenvalues of the covariance matrix for the a-, b-, and c-axis of olivine textures and then we assign the nearest orthogonal matrix to be the rotation matrix R:
We compute the rotation matrix R on the particles and further convert it to Euler angles for computation and memory efficiency. These properties need to be interpolated from particles to fields to be used in the material model. As a result, the anisotropic viscosity material model requires at least one particle in each cell so that all cells can have the texture parameters (Euler angles and eigenvalues) for constructing the rotation matrix R and compute the Hill coefficients. In the material model, the interpolated Euler angles are converted to the rotation matrix again. We use the same notation R to describe the rotation matrix used in the material model in the following paragraphs.
70
70
71
-
The inverse of the A tensor then needs to be rotated to the model reference frame. Since $inv(A)$ is the Kelvin notation of the rank-4 tensor, we apply the Kelvin notation representation of the R rotation matrix, $R_K$, on $inv(A_{ij})$:
71
+
The inverse of the A tensor then needs to be rotated to the model reference frame. Since $A_{ij}^{-1}$ is the Kelvin notation of the rank-4 tensor, we apply the Kelvin notation representation of the R rotation matrix, $R_K$, on $A_{ij}^{-1}$:
We save $\frac{1}{\gamma J(R\prime*\sigma_{ij}* R)^{(n-1)}}$ as the material model viscosity output, which we call the scalar viscosity. The tensorial part of anisotropic viscosity, $R_K * inv(A_{ij}) * R_K\prime$ is called the stress-strain director and is stored in the additional outputs. Due to the dependence of the scalar viscosity on stress and the stress is determined by the scalar viscosity from the previous time step, the computation of the scalar viscosity is not stable and its value naturally fluctuates up and down. This fluctuation ultimately causes a numerical instability associated with the prediction of the anisotropic viscosity. Therefore, we damp the scalar viscosity using a non-linear Newtonian iteration, and in each iteration, only half of the change in the scalar viscosity is applied until the result is stable.
93
+
$R'$ and $R_K'$ is the transpose of matrix $R$ and $R_K$ respectively. We save $\frac{1}{\gamma J(R'*\sigma_{ij}* R)^{(n-1)}}$ as the material model viscosity output, which we call the scalar viscosity. The tensorial part of anisotropic viscosity, $R_K * A_{ij}^{-1} * R_K'$ is called the stress-strain director and is stored in the additional outputs. Due to the dependence of the scalar viscosity on stress and the stress is determined by the scalar viscosity from the previous time step, the computation of the scalar viscosity is not stable and its value naturally oscillates. This fluctuation ultimately causes a numerical instability associated with the prediction of the anisotropic viscosity. Therefore, we damp the scalar viscosity using a non-linear Newton iteration, and in each iteration, only half of the change in the scalar viscosity is applied until the result is stable.
94
94
95
95
96
96
97
97
## Compiling requirement
98
98
99
-
Since the determinant of $A_{ij}$ is 0, $A_{ij}$ is a singular, non-invertible matrix. We find the MoorePenrose pseudoinverse of the matrix $A_{ij}$ as an approximate of the inverse of $A_{ij}$, using the SCALAPACK package in deal.ii. Thus it is required to link ASPECT with a deal.ii version with the scalapack package included and you can turn on the scalapack package when compiling deal.ii, for example, with `./candi.sh -j 8 --packages="p4est trilinos sundials scalapack dealii"`.
99
+
Since the determinant of $A_{ij}$ is 0, $A_{ij}$ is a singular, non-invertible matrix. We find the MoorePenrose pseudoinverse of the matrix $A_{ij}$ as an approximate of the inverse of $A_{ij}$, using the SCALAPACK package provided in deal.II. Thus it is required to link ASPECT with a deal.II version with the scalapack package included in order to run this cookbook. When using candi you can enable the scalapack package by including `scalapack` in the list of installed packages or uncommenting the line in `candi.cfg` that corresponds to the scalapack installation. Alternatively, you can install ScaLAPACK yourself and enable the configuration variable `DEAL_II_WITH_SCALAPACK` during the cmake configuration of deal.II.
100
100
101
-
## Input
101
+
## Model setup
102
102
103
103
The usage of the AV material model is demonstrated with a 3d simple shear box model, where its dimension is $1 \times 1 \times 1 $ (non-dimensionalized). The shear strain rate is set to
104
104
$0.5$. The origin is the center of the box, and one Olivine particle with 1000 grains sits at the origin to track CPO developments for computation of anisotropic viscosity parameters.
@@ -123,3 +123,6 @@ Note: These settings are similar to those used for simulations involving CPO alo
123
123
```
124
124
125
125
In the `CPO induced Anisotropic Viscosity` material model subsection, all parameters have reasonable default values and do not need to be manually specified unless customization is needed.
126
+
127
+
The output includes a particle property, anisotropic stress, which is a 3-by-3 matrix that can be visualized as a tensor, showing the components of the stress tensor. As a result of using the anisotropic viscosity material model, the components will be different in different directions.
0 commit comments