Plastic strain in reactive fluid transport rheology

fgueyd · March 11, 2025, 11:29am

Dear all,

I used the material model “reactive fluid transport “ based on the cookbook example developed by Daniel Douglas (Simplified Subduction Model with Parameterized Solid-Fluid Reactions — ASPECT 3.1.0-pre).

The base model is visco-plastic

I wish to track the plastic strain through time, as it is done in the classical visco-plastic rheology

I define a compositional field “plastic_strain” and define a Strain weakening mechanism so that the plastic strain is computed

The model run perfectly but the plastic strain is zero

It seems that when “reactive fluid transport” material model is used, even if the base model is visco-plastic, the plastic strain is not computed. I saw this message in the log
“Skipping plastic_strain composition solve because RHS is zero.”

Could you please let me know how to fix this problem

I attached hereafter two prm files
-one with a simple extension test with only visco-plastic rheology (plastic strain is well computed)
-one with the same parameters except the use of “reactive fluid transport” rheology (plastic strain is not computed)

Thank you for your feedbacks

Frédéric

-prm with a simple extension test with only visco-plastic rheology (plastic strain is well computed)

set Dimension = 2
set Use years in output instead of seconds = true
set Output directory = strain-VP-iterated
set End time = 0.5e6 #years

set Nonlinear solver scheme = iterated Advection and Stokes
set Nonlinear solver tolerance = 1e-1 # normaly 1e-5
set Max nonlinear iterations = 5 # 100 normaly

set Adiabatic surface temperature = 1600
#set Use operator splitting = true

subsection Geometry model
set Model name = box
subsection Box
set X repetitions = 4
set Y repetitions = 2
set X extent = 100e3
set Y extent = 50e3
end
end

subsection Mesh deformation
set Mesh deformation boundary indicators = top: free surface
subsection Free surface
set Surface velocity projection = vertical
end
set Additional tangential mesh velocity boundary indicators = left, right
end

set Pressure normalization = no

subsection Gravity model
set Model name = vertical
subsection Vertical
set Magnitude = 10.0
end
end

subsection Initial temperature model
set Model name = function
subsection Function
set Variable names = x,y,t
set Function constants = Ttop=273, Tbot=1473, ysize=50e3
set Function expression = Ttop+(Tbot-Ttop)*(ysize-y)/ysize
end
end

subsection Boundary temperature model
set Fixed temperature boundary indicators = top, bottom
set List of model names = function
subsection Function
set Function constants = Ttop=273, Tbot=1473, ysize=50e3
set Function expression = (y>ysize/2 ? Ttop : Tbot)
end
end

subsection Boundary velocity model
set Prescribed velocity boundary indicators = bottom y: function, left x: function, right x: function
subsection Function
set Variable names = x,y
set Function constants = v=0.06, xsize=100e3, ysize=50e3
set Function expression = (x<xsize/2 ? -v : v); (y<ysize/2 ? ysize/xsizev : -ysize/xsizev)
end
end

subsection Compositional fields
set Number of fields = 4
set Names of fields = peridotite, porosity, plastic_strain, seed
set Compositional field methods = field, field, field, field
end

subsection Initial composition model
set List of model names = function
subsection Function
set Variable names = x,y
set Function constants = xsize=100e3, ypos=30e3, rad=5e3, x1=40e3,x2=60e3,y1=5e3
set Function expression = 0; if(x>x1 && x<x2 && y<y1, 0.00, 0);0; if((x-xsize/2)^2+(y-ypos)^2<rad^2, 1, 0)
end
end

subsection Boundary composition model
set Fixed composition boundary indicators = top, bottom, left, right
set List of model names = initial composition
end

subsection Heating model
set List of model names = compositional heating

subsection Compositional heating
set Use compositional field for heat production averaging = 0
set Compositional heating values = 0.0
end
end

subsection Material model
set Model name = visco plastic
set Material averaging = harmonic average only viscosity

subsection Visco Plastic
set Densities = 3300 # kg/m3
set Thermal diffusivities = 1e-6
set Heat capacities = 1000
set Thermal expansivities = 0

set Viscous flow law                           = dislocation

set Prefactors for dislocation creep          = 2.7e-18, 1,1,1,  1e-15
set Activation energies for dislocation creep = 480e3,0,0,0,480e3
set Activation volumes for dislocation creep  = 1.1e-5,0,0,0,1.1e-5
set Stress exponents for dislocation creep    = 3.5,1,1,1,3.5

set Angles of internal friction                = 10.0, 10.,10, 1, 9.0 # degrees
set Cohesions                                  = 10e6 # Pa

set Minimum viscosity                          = 1e19 # Pa s
set Maximum viscosity                          = 1e25 # Pa s
set Strain weakening mechanism =    plastic weakening with plastic strain only
set Start plasticity strain weakening intervals  = 0.1
set End plasticity strain weakening intervals  = 1
set Cohesion strain weakening factors            =  0.25
set Friction strain weakening factors            = 0.25

end
end

subsection Mesh refinement
set Initial global refinement = 2 # normally 4
set Initial adaptive refinement = 2
set Time steps between mesh refinement = 1
set Strategy = composition, strain rate
set Refinement criteria merge operation = max
set Refinement criteria scaling factors = 1, 2
end

subsection Melt settings
#set Include melt transport = true
end

subsection Postprocess
set List of postprocessors = visualization
subsection Visualization
set Time between graphical output = 1
set Output mesh velocity = true
set List of output variables = material properties, strain rate, heating
subsection Material properties
set List of material properties = density, viscosity
end
end
end

-prm with the same parameters except the use of “reactive fluid transport” rheology (plastic strain is not computed)

set Dimension = 2
set Use years in output instead of seconds = true
set Output directory = strain-ReactiveFluidTransport-iterated
set End time = 0.5e6 #years