From c1cf3c09be7976f97cced6ab4913fd07fc3b7c65 Mon Sep 17 00:00:00 2001 From: albert-de-montserrat Date: Mon, 6 Jan 2025 14:44:20 +0100 Subject: [PATCH] up --- docs/make.jl | 7 +- docs/src/man/constitutive_equations.md | 134 ++++++++++++++++++----- docs/src/man/equations_APT.md | 14 +-- docs/src/man/equations_basic.md | 11 +- docs/src/man/equations_discretization.md | 14 ++- 5 files changed, 131 insertions(+), 49 deletions(-) diff --git a/docs/make.jl b/docs/make.jl index f5e9ba26f..3a1531576 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -79,7 +79,12 @@ makedocs(; "User guide"=> Any[ "Installation" => "man/installation.md", "Backend" => "man/backend.md", - "Equations" => "man/equations.md", + "Equations" => Any[ + "Governing equations" => "man/equations_basic.md", + "Constitutive equations" => "man/constitutive_equations.md", + "APT equations" => "man/equations_APT.md", + "Discretization" => "man/equations_discretization.md", + ], "Boundary conditions" => "man/boundary_conditions.md", "Advection" => "man/advection.md", ], diff --git a/docs/src/man/constitutive_equations.md b/docs/src/man/constitutive_equations.md index 4e8e4b03c..b79ae9419 100644 --- a/docs/src/man/constitutive_equations.md +++ b/docs/src/man/constitutive_equations.md @@ -17,12 +17,37 @@ $$ \begin{align} \boldsymbol{\dot\varepsilon} = \frac{1}{2\eta_{\text{eff}}}\boldsymbol{\tau} + -\frac{1}{2G} \frac{D\boldsymbol{\tau}}{Dt} + \dot\lambda\frac{\partial Q}{\partial \boldsymbol{\tau_{II}}} +\frac{1}{2G} \frac{D\boldsymbol{\tau}}{Dt} + \dot\lambda\frac{\partial Q}{\partial \boldsymbol{\tau}_{II}} \end{align} $$ -where $\eta_{\text{eff}}$ is the effective viscosity, $G$ is the elastic shear modulus, $\dot\lambda$ is the plastic multiplier, and $Q$ is the plastic flow potential. +where $\eta_{\text{eff}}$ is the effective viscosity, $G$ is the elastic shear modulus, $\dot\lambda$ is the plastic multiplier, and $Q$ is the plastic flow potential is: +$$ +\begin{align} +Q = \boldsymbol{\tau}_{II} - p\sin{\psi} +\end{align} +$$ +where $\psi$ is the dilation angle. + +## Effective viscosity + +The effective material of a material defined as: +$$ +\begin{align} +\eta_{\text{eff}} = \frac{1}{\frac{1}{\eta^{\text{diff}}} + \frac{1}{\eta^{\text{disl}}}} +\end{align} +$$ +where $\eta^{\text{diff}}$ and $\eta^{\text{disl}}$ are the diffusion and dislocation creep viscosities, which are computed from their respective strain rate equations: + +$$ +\begin{align} +\dot{ε}_{II}^{\text{diff}} = A^{\text{diff}} τ_{II}^{n^{\text{diff}}} d^{p} f_{H_2O}^{r^{\text{diff}}} \exp \left(- {{E^{\text{diff}} + PV^{\text{diff}}} \over RT} \right) \\ +\dot{ε}_{II}^{\text{disl}} = A^{\text{disl}} τ_{II}^{n^{\text{disl}}} f_{H_2O}^{r^{\text{disl}}} \exp \left(- {{E^{\text{disl}} + PV^{\text{disl}}} \over RT} \right) +\end{align} +$$ + +where $A$ material specific parameter, $n$ is the stress powerlaw exponent, $p$ is the negative defined grain size exponent, $f$ is the water fugacity, $r$ is the water fugacity exponent, $E$ is the activation energy, $PV$ is the activation volume, and $R$ is the universal gas constant. ## Elastic stress ### Method (1): Jaumann derivative @@ -46,24 +71,6 @@ $$ ### Method (2): Euler-Rodrigues rotation -```julia -# from Anton's talk -@inline Base.@propagate_inbounds function rotate_elastic_stress3D(ωi, τ, dt) - # vorticity - ω = √(sum(x^2 for x in ωi)) - # unit rotation axis - n = SVector{3, Float64}(inv(ω) * ωi[i] for i in 1:3) - # integrate rotation angle - θ = dt * 0.5 * ω - # Euler Rodrigues rotation matrix - R = rodrigues_euler(θ, n) - # rotate tensor - τij = voigt2tensor(τ) - τij_rot = R * (τij * R') - tensor2voigt(τij_rot) -end -``` - 1. Compute unit rotation axis $$ \begin{align} @@ -87,11 +94,11 @@ $$ \begin{align} \boldsymbol{c} = \boldsymbol{n}\sin{\theta} \\ \boldsymbol{R_1} = - \left[ - c0 -c_3 c_2 \\ - c_3 c0 -c_1 \\ - -c_2 c_1 c0 \\ - \right] \\ + \begin{bmatrix} + c0 & -c_3 & c_2 \\ + c_3 & c0 & -c_1 \\ + -c_2 & c_1 & c0 \\ + \end{bmatrix} \\ \boldsymbol{R_2} = (1-\cos{\theta}) \boldsymbol{n} \boldsymbol{n}^T \\ \boldsymbol{R}=\boldsymbol{R_1}+\boldsymbol{R_2} \end{align} @@ -107,11 +114,11 @@ $$ ## Plastic formulation -JustRelax.jl implements the regularised plasticity model from (refs). In this formulation, the yield function is given by +JustRelax.jl implements the regularised plasticity model from [Thibault et al 2021](https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2021GC009675). In this formulation, the yield function is given by $$ \begin{align} -F = \tau_y - \left( P \sin{\phi} + C \cos{\phi} + \eta_{\text{reg}} + \dot\lambda \right) \leq 0 +F = \tau_y - \left( P \sin{\phi} + C \cos{\phi} + \dot\lambda \eta_{\text{reg}}\right) \leq 0 \end{align} $$ @@ -119,8 +126,77 @@ where $\eta_{\text{reg}}$ is a regularization term and $$ \begin{align} -\dot\lambda = \frac{F}{x} +\dot\lambda = \frac{F^{\text{trial}}}{\eta + \eta_{\text{reg}} + K \Delta t \sin{\psi}\sin{\phi}} \end{align} $$ -Note that since we are using an iterative method to solve the APT Stokes equation, the non-linearities are dealt by the iterative scheme. Othwersie, one would need to solve a non-linear problem to compute $\dot\lambda$, which typically requires to compute $\frac{\partial F}{\partial \dot\lambda}$ \ No newline at end of file +Note that since we are using an iterative method to solve the APT Stokes equation, the non-linearities are dealt by the iterative scheme. Othwersie, one would need to solve a non-linear problem to compute $\dot\lambda$, which requires to compute $\frac{\partial F}{\partial \dot\lambda}$ + +# Selecting the constitutive model in JustRelax.jl + +All the local calculations corresponding to the effective rheology are implemented in GeoParams.jl. The composite rheology is implemented using the `CompositeRheology` object. An example of how to set up the a visco-elasto-viscoplastic rheology is shown below: + +```julia +# elasticity +el = ConstantElasticity(; + G = 40e9, # shear modulus [Pa] + ν = 0.45, # Poisson coeficient +) +# Olivine dislocation law from Hirth and Kohlstedt 2003 +disl_wet_olivine = SetDislocationCreep(Dislocation.wet_olivine1_Hirth_2003) +# Olivine diffusion law from Hirth and Kohlstedt 2003 +diff_wet_olivine = SetDiffusionCreep(Diffusion.wet_olivine_Hirth_2003) +# plasticity +ϕ = 30 # friction angle +Ψ = 5 # dilation angle +C = 1e6 # cohesion [Pa] +η_reg = 1e16 # viscosity regularization term [Pa s] +pl = DruckerPrager_regularised(; C = C, ϕ = ϕ, η_vp=η_reg, Ψ=Ψ) +# composite rheology +rheology = CompositeRheology( + (el, disl_wet_olivine, diff_wet_olivine, pl) +) +``` + +`rheology` then needs to be passed into a `MaterialParams` object with the help of `SetMaterialParams`: + +```julia +rheology = ( + SetMaterialParams(; + Phase = 1, + CompositeRheology = rheology, + # other material properties here + ), +) +``` + +### Multiple rheology phases + +It is common in geodynamic models to have multiple rheology phases. In this case, we just need to build a tuple containing every single material phase properties: +```julia +rheology = ( + SetMaterialParams(; + Phase = 1, + CompositeRheology = rheology_one, + # other material properties here + ), + SetMaterialParams(; + Phase = 2, + CompositeRheology = rheology_two, + # other material properties here + ), +) +``` + +### Computing the effective viscosity + +The effective viscosity is computed internally during the Stokes solver. However, it can also be computed externally with the `compute_viscosity!` function as follows: + +```julia +# Rheology +args = (T=T, P=P, dt = dt) # or (T=thermal.Tc, P=stokes.P, dt=dt) +viscosity_cutoff = (1e18, 1e23) +compute_viscosity!(stokes, phase_ratios, args, rheology, viscosity_cutoff) +``` + +where `T` and `P` are the temperature and pressure fields defined at the **cell centers**, `dt` is the time step, and `phase_ratios` is an object containing the phase ratios corresponding to each material phase, of each cell. \ No newline at end of file diff --git a/docs/src/man/equations_APT.md b/docs/src/man/equations_APT.md index 5282ef525..823b8b24d 100644 --- a/docs/src/man/equations_APT.md +++ b/docs/src/man/equations_APT.md @@ -1,7 +1,7 @@ -The pseudo-transient method consists in augmenting the right-hand-side of the target PDE with a pseudo-time derivative (where $\psi$ is the pseudo-time) of the primary variables. We then solve the resulting system of equations with an iterative method. The pseudo-time derivative is then gradually reduced, until the original PDE is solved and the changes in the primary variables are below a preset tolerance. +The Accelerated Pseudo-Transient (APT) method consists in augmenting the right-hand-side of the target PDE with a pseudo-time derivative (where $\psi$ is the pseudo-time) of the primary variables. We then solve the resulting system of equations with an iterative method. The pseudo-time derivative is then gradually reduced, until the original PDE is solved and the changes in the primary variables are below a preset tolerance. ## Heat diffusion -The pseudo-transient heat-diffusion equation is: +The APT heat-diffusion equation is: $$ \begin{align} @@ -9,7 +9,7 @@ $$ \end{align} $$ -We use a second order pseudo-transient scheme were continuation is also done on the flux, so that: +We use a second order APT scheme were continuation is also done on the flux, so that: $$ \begin{align} @@ -19,7 +19,7 @@ $$ ## Stokes equations -For example, the pseudo-transient formulation of the Stokes equations yields: +For example, the APT formulation of the Stokes equations yields: $$ \begin{align} @@ -34,7 +34,7 @@ $$ $$ ## Constitutive equations -A pseudo-transient continuation is also done on the constitutive law: +A APT continuation is also done on the constitutive law: $$ \begin{align} @@ -60,7 +60,7 @@ $$ where the P-wave $\widetilde{V}=V_p$ is the characteristic velocity scale for Stokes, and $Re$ is the Reynolds number. -### Physical parameters + ### Pseudo-transient parameters diff --git a/docs/src/man/equations_basic.md b/docs/src/man/equations_basic.md index 4ccc62917..89893671d 100644 --- a/docs/src/man/equations_basic.md +++ b/docs/src/man/equations_basic.md @@ -1,3 +1,5 @@ +# Governing equations + ## Stokes equations Deformation of compressible viscous flow is desribed by the equations of conservation of momentum and mass: @@ -16,6 +18,8 @@ $$ where $\boldsymbol{\tau}$ is the deviatoric stress tensor, $p$ is pressure, $f$ is the external forces vector, $\boldsymbol{v}$ is the velocity vector, $\beta$ is the compressibility coefficient, $\alpha$ is the thermal expansivity coefficient and $T$ is temperature. +## Constitutive equation + To close the previous systems of equations, we further need the constitutive relationship between stress and deformationl. In its simplest linear form this is: $$ @@ -25,7 +29,6 @@ $$ $$ where $\eta$ is the shear viscosity and $\boldsymbol{\dot\varepsilon}$ is the deviatoric strain tensor. - ## Heat diffusion The pseudo-transient heat-diffusion equation is: @@ -36,9 +39,3 @@ $$ $$ where $\rho$ is density, $C_p$ is specific heat capacity, $\kappa$ is thermal conductivity, $T$ is temperature $\boldsymbol{\tau}:\boldsymbol{\dot\varepsilon}$ is the energy dissipated by viscous deformation (shear heating), $\alpha T(\boldsymbol{v} \cdot \nabla P)$ is adiabatic heating, and $H$ is the sum any other source term, such as radiogenic heat production. - -# Discretization of the governing equations - -We discretize both the Stokes and heat diffussion equations using a Finite Differences approach on a staggered grid (ref Taras book here), as sketched below: - -![Staggered Velocity Grid](../assets/StaggeredVelocity.png) \ No newline at end of file diff --git a/docs/src/man/equations_discretization.md b/docs/src/man/equations_discretization.md index a17a229a4..2edf420a4 100644 --- a/docs/src/man/equations_discretization.md +++ b/docs/src/man/equations_discretization.md @@ -2,8 +2,6 @@ We discretize both the Stokes and heat diffussion equations using a Finite Differences approach on a staggered grid (ref Taras book here). - - ## Heat diffusion The heat diffusion equation is discretized on a staggered grid as sketched below: @@ -39,8 +37,7 @@ The equations of conservation of mass and momentum are discretized on a staggere ![](../assets/stokes_stag2D.png) -where dotted lines represent the velocity ghost nodes. The APT Stokes equations are then discretized as follows: - +where dotted lines represent the velocity ghost nodes. +### Conservation of momentum + $$ \begin{align} \widetilde{\rho}\frac{u^{n+1}_x - u^n_x}{\Delta\psi} + \nabla\cdot\boldsymbol{\tau} - @@ -72,4 +71,9 @@ $$ \frac{1}{\Delta t} \left( \beta (p^{t+\Delta t} - p^{t}) + \alpha (T^{t+\Delta t} - T^{t}) \right) \\ \end{align} -$$ \ No newline at end of file +$$ + +### Conservation of mass + + +### Constitutive equation