Implementing new constitutive model

[amirhosseinrezaei9776]

Hi!
I'm working on implementing a Mohr-Coulomb Strain Softening constitutive model in MaterialPointSolver.jl, which involves introducing a few hyperparameters, modifying the cohesion (c) formula, and adding a softening formulation for the friction angle (φ).
I tried to modify the mohrcoulomb.jl and property.jl files accordingly, but unfortunately the model does not run. I suspect I may have missed a key step or failed to properly integrate the changes into the solver framework.
Could you kindly guide me on the correct way to implement a custom constitutive model in your framework? Specifically, how should I introduce additional parameters and modify the strength parameters over time or strain history?
Any help or advice you can offer would be greatly appreciated.

[ZenanH]

Let's assume the constitutive model is like this:

$\Delta\sigma_{xx} = \Delta\epsilon_{xx} + c_1$
$\Delta\sigma_{yy} = \Delta\epsilon_{yy} + c_2$
$\Delta\sigma_{xy} = \Delta\epsilon_{xy} + c_3$

$\sigma_{xx}^{t+1} = \sigma_{xx}^{t+1} + \Delta\sigma_{xx}$
$\sigma_{yy}^{t+1} = \sigma_{yy}^{t+1} + \Delta\sigma_{yy}$
$\sigma_{xy}^{t+1} = \sigma_{xy}^{t+1} + \Delta\sigma_{xy}$

where $c_1$, $c_2$, $c_3$ are the new parameters you want to use.

---

Firstly, generate args, grid, mp, attr, and bc in the main.jl (main file).

If you don't know about them, please check the files under examples folder and documentation.

Now, you should execute materialpointsolver!(args, grid, mp, attr, bc).

However, to use your modified constitutive model, you should write them in another file, e.g., funcs.jl and add these (refer to here for the constitutive model and here for the MPM procedure):

using KernelAbstractions

# create extension for `attr` to save your new parameters
struct NewPara{T1<:AbstractArray} <: UserPropertyExtra
    c1::T1
    c2::T1
    c3::T1
end

@user_struct NewPara

@kernel function userdefined!(
    mp  ::DeviceParticle2D{T1, T2}, # normally you only need these two variables, but you can add more
    attr::  DeviceProperty{T1, T2}
) where {T1, T2}
    ix = @index(Global)
    if ix ≤ mp.np
        id = attr.nid[ix] # determine particle property
        c1 = attr.ext.c1[id] # note here is using `id`, not `ix` 
        c2 = attr.ext.c2[id] # if your case is easy, you can just define the values here...
        c3 = attr.ext.c3[id]    
        Δσxx = mp.Δϵijs[ix, 1] + c1
        Δσyy = mp.Δϵijs[ix, 2] + c2
        Δσxy = mp.Δϵijs[ix, 4] + c3 # mp.Δϵijs[ix, 3] is Δϵ_zz for 2d plane strain problem
        mp.σij[ix, 1] += Δσxx
        mp.σij[ix, 2] += Δσyy
        mp.σij[ix, 4] += Δσxy
        mp.σm[ix] = (mp.σij[ix, 1] + mp.σij[ix, 2]) * T2(0.5)
    end
end

function newprocedure!(
    args::     DeviceArgs{T1, T2}, 
    grid::     DeviceGrid{T1, T2}, 
    mp  :: DeviceParticle{T1, T2}, 
    attr:: DeviceProperty{T1, T2},
    bc  ::DeviceVBoundary{T1, T2},
    ΔT  ::T2,
    Ti  ::T2,
        ::Val{:OS},
        ::Val{:MUSL}
) where {T1, T2}
    G = Ti < args.Te ? args.gravity / args.Te * Ti : args.gravity
    dev = getBackend(Val(args.device))
    # MPM procedure
    resetgridstatus_OS!(dev)(ndrange=grid.ni, grid)
    resetmpstatus_OS!(dev)(ndrange=mp.np, grid, mp, Val(args.basis))
    P2G_OS!(dev)(ndrange=mp.np, grid, mp, G)
    solvegrid_OS!(dev)(ndrange=grid.ni, grid, bc, ΔT, args.ζs)
    doublemapping1_OS!(dev)(ndrange=mp.np, grid, mp, attr, ΔT, args.FLIP, args.PIC)
    doublemapping2_OS!(dev)(ndrange=mp.np, grid, mp)
    doublemapping3_OS!(dev)(ndrange=grid.ni, grid, bc, ΔT)
    # F-bar based volumetric locking elimination approach
    if args.MVL == false
        G2P_OS!(dev)(ndrange=mp.np, grid, mp, ΔT)
    else
        G2Pvl1_OS!(dev)(ndrange=mp.np, grid, mp)
        fastdiv!(dev)(ndrange=grid.ni, grid)
        G2Pvl2_OS!(dev)(ndrange=mp.np, grid, mp, ΔT)
    end
    # update stress status
    liE!(dev)(ndrange=mp.np, mp, attr)
    Ti ≥ args.Te && userdefined!(dev)(ndrange=mp.np, mp, attr) # only change here to use your model
    return nothing
end

Then, in the main.jl:

using MaterialPointGenerator
using MaterialPointSolver
using MaterialPointVisualizer

include(joinpath(@__DIR__, "funcs.jl"))

args = ...
grid = ...
mp   = ...

tmp_attr = NewPara{AbstractArray}([c1], [c2], [c3]) # adjust them if you have different particle partitions
ext_attr = user_adapt(Array, tmp_attr)
attr = UserProperty(
    ...
    ext = ext
) # then your parameters can be used in your model

bc = ...

materialpointsolver!(attr, grid, mp, attr, bc, workflow=newprocedure!) # it is the function's name

[ZenanH]

Note this will only be effective when you are using MUSL, otherwise it will cause error since the workflow, newprocedure!, has no definition for other stress update schemes. However, you can do it in the funcs.jl if you want.