underworldcode / UWGeodynamics

Underworld Geodynamics
Other
80 stars 32 forks source link

Viscous strain softening #157

Closed raghugram closed 2 years ago

raghugram commented 4 years ago

I am working on viscous strain softening using UWGeo.

Similar to plastic strain softening except that the pre-exponential factor function as given in Ros et. al. 2017 supplementary material(page 5) is depenedent on the temperature field and accumulated strain

Rather than the total accumulated strain, we just consider the accumulated viscous strain in the model.

In _model.py

1) I first created a swarm particle to save the accumulated viscous strain. 2) Calculate viscous strain from the viscosity of viscousEta(from _viscFn which calls just the _ViscousEta() function.

    viscousStrainRate2I = fn.tensor.second_invariant(2. * self._viscFn * self.strainRate)
    viscousStrainRate2I =  (0.5 * viscousStrainRate2I /self._viscFn)
    viscousStrainIncrement = dt  * viscousStrainRate2I.evaluate(self.swarm)
    self.viscousStrain.data[:] += viscousStrainIncrement

In _rheology.py Added the initialize the viscousStrain wherever needed and 3) In _effectiveViscosity() function, I added

    I2v_acc = self.viscousStrain
    alpha1 = fn.Function.convert(1.)
    alpha2 = fn.Function.convert(30.)
    epsilon1 = fn.Function.convert(0.)
    epsilon2 = fn.Function.convert(1.)
    Tmax = nd(1473.* u.degK)
    Tmax= fn.Function.convert(Tmax)
    Tmin = nd(1073.* u.degK)
    Tmin = fn.Function.convert(Tmin)
    def _viscous_tempdependent_alpha(tField):
        Tmin = fn.Function.convert(Tmin)    
        Talphacondition = [(tField <= Tmin, alpha2),
           (tField >= Tmax, alpha1),
           (True, fn.math.exp(fn.math.log(alpha2)*(Tmax/(Tmin - Tmax)) * ((tField - Tmax)/(Tmax))))]             
        return fn.branching.conditional(Talphacondition)

    alphaValcondition = [(I2v_acc <= epsilon1, alpha1),
           (I2v_acc >= epsilon2, alpha2),
           (True, alpha1 +((_viscous_tempdependent_alpha(T)  - alpha1) /
                                  (epsilon2 - epsilon1)) *  (I2v_acc - epsilon1))]
    alphaVal = fn.branching.conditional(alphaValcondition)       

    mu_eff = f * 0.5 * A**(-1.0 / n)
    mu_eff *= alphaVal**(-1.0 / n)
    ...
    ....

I am not sure if this is working?I2v_acc is the swarmvariable and _viscous_tempdependent_alpha() function input is the mesh variable T.

Does the behind the scenes numpy handling of swarm and mesh variables make this compatible?

Also, can I do the same in the front end, when I develop the model? Like, define viscosityFn as a function of accumulated strain and temperature in UWGeo

Ref: 1)Ros, E., Pérez‐Gussinyé, M., Araújo, M., Thoaldo Romeiro, M., Andrés‐Martínez, M. and Morgan, J.P., 2017. Lower crustal strength controls on melting and serpentinization at magma‐poor margins: Potential implications for the South Atlantic. Geochemistry, Geophysics, Geosystems, 18(12), pp.4538-4557. 2)Andrés-Martínez, M., Pérez Gussinyé, M., Armitage, J. and Morgan, J.P., 2019. Thermomechanical implications of sediment transport for the architecture and evolution of continental rifts and margins. Tectonics, 38(2), pp.641-665.

Thanks, Raghu Ram

rbeucher commented 4 years ago

Hi @raghugram,

I have added the option to do viscous strain softening. Have a look at the viscous_strain_softening branch.

I have added a series of arguments to the ViscousCreep class. You can pass alpha1, alpha2, epsilon1, epsilon2 and alphaFn.

alphaFn would be your temperature function in that case. This is strain softening so the dependence on strain is handled internally. You should be able to pass your _viscous_tempdependent_alpha function as alphaFn.

It is untested...have a look. It is pretty much what you had suggested previously. Maybe we can work from that and see how it goes.

Romain

rbeucher commented 4 years ago

Hi @raghugram ,

Have you played with the changes I have made/ I am thinking of merging into dev and add it to the next release.

rbeucher commented 4 years ago

I had a look at your email, sorry about the delay. I am fixed the mistake with the viscousStrain initialisation. I looked at the example you sent me. I think it needs a bit of work. You can define the viscosity without initialising the temperature field.

Setting Model.temperature = True Initialise the temperature. Then the effective viscosity is calculated using the Model.temperature field, as long as you define your function using Model.temperature

Would be good to have a use case though. I think there is a paper from Huismans et al that looks at viscous strain weakening.

I am keeping things in the viscous_strain_weakening branch. Feel free to fork and make a pull request.

rbeucher commented 4 years ago

I have put a copy of the example you sent me by email. https://github.com/underworldcode/UWGeodynamics/blob/viscous_strain_softening/examples/1_33_Viscous_Strain_Softening.ipynb