Skip to content
Permalink
Browse files

Merge pull request #13477 from WilkAndy/plsink_multiplasticity_12962

MultiParameterPlasticityStressUpdate now works with DiracKernels
  • Loading branch information...
permcody committed Jun 3, 2019
2 parents 51d07de + 1f4d3f2 commit 2f849e3fe80be6878c0b0c50ee5850a6af6fc507
Binary file not shown.
@@ -0,0 +1,238 @@
# Example: Injection into a uniform aquifer 10 x 10 x 5 km
# Drucker-Prager deformation
# Darcy flow
gravity = -9.81
solid_density = 2350
fluid_density = 1000
porosity0 = 0.1
[Mesh]
type = GeneratedMesh
dim = 3
xmin = 0
xmax = 1e4
ymin = 0
ymax = 1e4
zmax = 0
zmin = -5e3
nx = 2
ny = 2
nz = 2
[]

[GlobalParams]
PorousFlowDictator = dictator
gravity = '0 0 ${gravity}'
displacements = 'disp_x disp_y disp_z'
strain_at_nearest_qp = true
[]

[Modules]
[./FluidProperties]
[./simple_fluid]
type = SimpleFluidProperties
thermal_expansion = 0 # Not doing a thermal simulation
bulk_modulus = 2E9
density0 = ${fluid_density}
viscosity = 5E-4
[../]
[../]
[]

[PorousFlowFullySaturated]
coupling_type = HydroMechanical
porepressure = pp
dictator_name = dictator
fp = simple_fluid
[]

[Variables]
[./disp_x]
[../]
[./disp_y]
[../]
[./disp_z]
[../]
[./pp]
scaling = 1E6
[./InitialCondition]
type = FunctionIC
function = ini_pp
[../]
[../]
[]

[Functions]
[./ini_stress]
type = ParsedFunction
value = '-${gravity} * z * (${solid_density} - ${fluid_density}) * (1.0 - ${porosity0})' # initial effective stress that should result from weight force
[../]

[./ini_pp]
type = ParsedFunction
value = '${gravity} * z * ${fluid_density} + 1E5'
[../]
[]

[BCs]
[./p_top]
type = FunctionPresetBC
variable = pp
boundary = front
function = ini_pp
[../]

[./x_roller]
type = PresetBC
variable = disp_x
boundary = 'left right'
value = 0
[../]
[./y_roller]
type = PresetBC
variable = disp_y
boundary = 'top bottom'
value = 0
[../]
[./z_confined]
type = PresetBC
variable = disp_z
boundary = 'back front'
value = 0
[../]
[]

[UserObjects]
[./pls_total_outflow_mass]
type = PorousFlowSumQuantity
[../]

# Cohesion
[./mc_coh]
type = TensorMechanicsHardeningConstant
value = 6.0E6
[../]

# Friction angle
[./mc_phi]
type = TensorMechanicsHardeningConstant
value = 35.0
convert_to_radians = true
[../]

# Dilation angle
[./mc_psi]
type = TensorMechanicsHardeningConstant
value = 2
convert_to_radians = true
[../]

# Drucker-Prager objects
[./dp]
type = TensorMechanicsPlasticDruckerPragerHyperbolic
mc_cohesion = mc_coh
mc_friction_angle = mc_phi
mc_dilation_angle = mc_psi
yield_function_tolerance = 1E-3
internal_constraint_tolerance = 1E-6
[../]

# Tensile strength
[./tens]
type = TensorMechanicsHardeningConstant
value = 3.0E6
[../]

# Compressive strength (cap on yield envelope)
[./compr_all]
type = TensorMechanicsHardeningConstant
value = 1E10
[../]
[]

[Materials]
[./strain]
type = ComputeIncrementalSmallStrain
eigenstrain_names = eigenstrain_all
[../]

[./eigenstrain_all]
type = ComputeEigenstrainFromInitialStress
initial_stress = 'ini_stress 0 0 0 ini_stress 0 0 0 ini_stress'
eigenstrain_name = eigenstrain_all
[../]

[./elasticity_tensor]
type = ComputeIsotropicElasticityTensor
bulk_modulus = 3.3333E9
shear_modulus = 2.5E9
[../]

[./dp_mat]
type = CappedDruckerPragerStressUpdate
DP_model = dp
tensile_strength = tens
compressive_strength = compr_all
smoothing_tol = 1E5
yield_function_tol = 1E-3
tip_smoother = 0
[../]

[./stress]
type = ComputeMultipleInelasticStress
inelastic_models = dp_mat
[../]

# Permeability
[./permeability]
type = PorousFlowPermeabilityConst
permeability = '1E-13 0 0 0 1E-13 0 0 0 1E-13'
[../]

# Porosity
[./porosity]
type = PorousFlowPorosity
porosity_zero = ${porosity0}
biot_coefficient = 1.0
solid_bulk = 1.0 # Required but irrelevant when biot_coefficient is unity
mechanical = true
fluid = true
[../]

# Density of saturated rock
[./density]
type = PorousFlowTotalGravitationalDensityFullySaturatedFromPorosity
rho_s = ${solid_density}
[../]
[]

[DiracKernels]
[./pls]
type = PorousFlowPolyLineSink
variable = pp
SumQuantityUO = pls_total_outflow_mass
point_file = two_nodes.bh
function_of = pressure
fluid_phase = 0
p_or_t_vals = '0 1E7'
fluxes = '-1.59 -1.59'
[../]
[]

[Preconditioning]
[./usual]
type = SMP
full = true
[../]
[]

[Executioner]
solve_type = Newton
type = Transient
dt = 1E6
end_time = 1E6
nl_rel_tol = 1E-7
[]

[Outputs]
exodus = true
[]
@@ -362,4 +362,13 @@
issues = "#12541"
design = 'porous_flow/tests/dirackernels/dirackernels_tests.md porous_flow/sinks.md'
[../]
[./injection_with_plasticity]
type = 'Exodiff'
input = 'injection_with_plasticity.i'
exodiff = 'injection_with_plasticity_out.e'
threading = '!pthreads'
requirement = 'It shall be possible to use PorousFlow DiracKernels in simulations involving plastic deformation'
design = 'porous_flow/sinks.md PorousFlowPeacemanBorehole.md'
issues = '#12962'
[../]
[]
@@ -0,0 +1,2 @@
1 5000.0 5000.0 -1400.0
1 5000.0 5000.0 -1500.0
@@ -8,7 +8,7 @@
//* https://www.gnu.org/licenses/lgpl-2.1.html

#include "MultiParameterPlasticityStressUpdate.h"
#include "Conversion.h" // for stringify
#include "Conversion.h" // for stringify

// libMesh includes
#include "libmesh/utility.h" // for Utility::pow
@@ -136,8 +136,7 @@ void
MultiParameterPlasticityStressUpdate::propagateQpStatefulProperties()
{
_plastic_strain[_qp] = _plastic_strain_old[_qp];
for (unsigned i = 0; i < _num_intnl; ++i)
_intnl[_qp][i] = _intnl_old[_qp][i];
std::copy(_intnl_old[_qp].begin(), _intnl_old[_qp].end(), _intnl[_qp].begin());
}

void
@@ -151,6 +150,14 @@ MultiParameterPlasticityStressUpdate::updateState(RankTwoTensor & strain_increme
bool compute_full_tangent_operator,
RankFourTensor & tangent_operator)
{
// Size _yf[_qp] appropriately
_yf[_qp].assign(_num_yf, 0);
// _plastic_strain and _intnl are usually sized appropriately because they are stateful, but this
// Material may be used from a DiracKernel where stateful materials are not allowed. The best we
// can do is:
if (_intnl[_qp].size() != _num_intnl)
initQpStatefulProperties();

initializeReturnProcess();

if (_t_step >= 2)

0 comments on commit 2f849e3

Please sign in to comment.
You can’t perform that action at this time.