Confusion of Units and Ability to Apply Fixed Displacement

[jmanthony3]

I have been playing around with this package for a couple of weeks now and have come to really enjoy using it. I am working on a project that is similar to the von Mises Plasticity example, and that this package is open-source is really attractive.

So far, I want to modify the von Mises Plasticity example to perform a representative volume element (RVE) test by applying a fixed displacement in each timestep (like a strain-controlled loading condition). But I am supremely confused whether each element of the tangent stiffness matrix, K has units of stress or energy (Ke[i, j] += δϵ ⊡ D ⊡ Δϵ * dΩ at line 150 in the example).

Below is a MWE which includes some modifications to the von Mises Plasticity example: e. g. adding ϵ as an attribute to MaterialState and printing a comparison of the applied and actual displacement in each timestep. #788 was instrumental to strafe through the global displacement vector, u. Maybe this could be officially implemented?

Questions

1. Is the traction vector units of force or stress?
2. Is each element of the tangent stiffness matrix, K units of stress or energy?
3. How would I modify the MWE such that the maximum displacement is 0.1 m (integral of velocity through time) instead of the measured 0.6 m because of the traction vector?

MWE

Code

# von Mises plasticity example (plain program)
# https://ferrite-fem.github.io/Ferrite.jl/stable/tutorials/plasticity/

using Ferrite, FerriteGmsh, FerriteMeshParser, Tensors, SparseArrays, LinearAlgebra, Printf
using NumericalIntegration, WriteVTK

struct J2Plasticity{T, S <: SymmetricTensor{4, 3, T}}
    G::T  # Shear modulus
    K::T  # Bulk modulus
    σ₀::T # Initial yield limit
    H::T  # Hardening modulus
    Dᵉ::S # Elastic stiffness tensor
end;

function J2Plasticity(E, ν, σ₀, H)
    δ(i,j) = i == j ? 1.0 : 0.0 # helper function
    G = E / 2(1 + ν)
    K = E / 3(1 - 2ν)

    Isymdev(i,j,k,l) = 0.5*(δ(i,k)*δ(j,l) + δ(i,l)*δ(j,k)) - 1.0/3.0*δ(i,j)*δ(k,l)
    temp(i,j,k,l) = 2.0G *( 0.5*(δ(i,k)*δ(j,l) + δ(i,l)*δ(j,k)) + ν/(1.0-2.0ν)*δ(i,j)*δ(k,l))
    Dᵉ = SymmetricTensor{4, 3}(temp)
    return J2Plasticity(G, K, σ₀, H, Dᵉ)
end;

struct MaterialState{T, S <: SecondOrderTensor{3, T}}
    # Store "converged" values
    ϵ::S # total strain
    ϵᵖ::S # plastic strain
    σ::S # stress
    k::T # hardening variable
end

function MaterialState()
    return MaterialState(
                zero(SymmetricTensor{2, 3}),
                zero(SymmetricTensor{2, 3}),
                zero(SymmetricTensor{2, 3}),
                0.0)
end

function vonMises(σ)
    s = dev(σ)
    return sqrt(3.0/2.0 * s ⊡ s)
end;

function compute_stress_tangent(ϵ::SymmetricTensor{2, 3}, material::J2Plasticity, state::MaterialState)
    # unpack some material parameters
    G = material.G
    H = material.H

    # We use (•)ᵗ to denote *trial*-values
    σᵗ = material.Dᵉ ⊡ (ϵ - state.ϵᵖ) # trial-stress
    sᵗ = dev(σᵗ)         # deviatoric part of trial-stress
    J₂ = 0.5 * sᵗ ⊡ sᵗ   # second invariant of sᵗ
    σᵗₑ = sqrt(3.0*J₂)   # effective trial-stress (von Mises stress)
    σʸ = material.σ₀ + H * state.k # Previous yield limit

    φᵗ  = σᵗₑ - σʸ # Trial-value of the yield surface

    if φᵗ < 0.0 # elastic loading
        return σᵗ, material.Dᵉ, MaterialState(ϵ, state.ϵᵖ, σᵗ, state.k)
    else # plastic loading
        h = H + 3G
        μ =  φᵗ / h   # plastic multiplier

        c1 = 1 - 3G * μ / σᵗₑ
        s = c1 * sᵗ           # updated deviatoric stress
        σ = s + vol(σᵗ)       # updated stress

        # Compute algorithmic tangent stiffness ``D = \frac{\Delta \sigma }{\Delta \epsilon}``
        κ = H * (state.k + μ) # drag stress
        σₑ = material.σ₀ + κ  # updated yield surface

        δ(i,j) = i == j ? 1.0 : 0.0
        Isymdev(i,j,k,l)  = 0.5*(δ(i,k)*δ(j,l) + δ(i,l)*δ(j,k)) - 1.0/3.0*δ(i,j)*δ(k,l)
        Q(i,j,k,l) = Isymdev(i,j,k,l) - 3.0 / (2.0*σₑ^2) * s[i,j]*s[k,l]
        b = (3G*μ/σₑ) / (1.0 + 3G*μ/σₑ)

        Dtemp(i,j,k,l) = -2G*b * Q(i,j,k,l) - 9G^2 / (h*σₑ^2) * s[i,j]*s[k,l]
        D = material.Dᵉ + SymmetricTensor{4, 3}(Dtemp)

        # Return new state
        Δϵᵖ = 3/2 * μ / σₑ * s # plastic strain
        ϵᵖ = state.ϵᵖ + Δϵᵖ    # plastic strain
        k = state.k + μ        # hardening variable
        return σ, D, MaterialState(ϵ, ϵᵖ, σ, k)
    end
end

function symmetrize_lower!(K)
    for i in 1:size(K,1)
        for j in i+1:size(K,1)
            K[i,j] = K[j,i]
        end
    end
end;

# https://github.com/Ferrite-FEM/Ferrite.jl/discussions/788
function vertexdofs(dh::DofHandler, vertexid::VertexIndex)

    cellid, lvidx = vertexid
    nfields = length(dh.field_names)
    sdh = dh.subdofhandlers[dh.cell_to_subdofhandler[cellid]]
    local_vertex_dofs = Int[]
    for ifield in 1:length(sdh.field_names)
        offset = Ferrite.field_offset(sdh, ifield)
        # field_dim = Ferrite.getfielddim(sdh, ifield) # deprecated
        field_dim = Ferrite.n_components(sdh, ifield)
        _field_ip = sdh.field_interpolations[ifield]
        if _field_ip isa Ferrite.VectorizedInterpolation
            field_ip = _field_ip.ip
        else
            field_ip = _field_ip
        end
        vert = Ferrite.vertexdof_indices(field_ip)[lvidx]
        
        for vdof in vert, d in 1:field_dim 
            push!(local_vertex_dofs, (vdof-1)*field_dim + d + offset)
        end
    end

    dofs = zeros(Int, ndofs_per_cell(dh, cellid))
    celldofs!(dofs, dh, cellid)

    return dofs[local_vertex_dofs]
end

function nodeid_to_vertexindex(grid::Grid, nodeid::Int)
    for (cellid, cell) in enumerate(grid.cells)
        for (i, nodeid2) in enumerate(cell.nodes)
            if nodeid == nodeid2
                return VertexIndex(cellid, i)
            end
        end
    end
    error("Node $(nodeid) does not belong to any cell")
end

vectorstrafe(grid, dh, i) = vertexdofs(dh, nodeid_to_vertexindex(grid, i))

function solve()
    # Define material parameters
    E = 200.0e3     # [MPa]
    H = E/20        # [MPa]
    ν = 0.3         # [-]
    σ₀ = 200.0      # [MPa]
    material = J2Plasticity(E, ν, σ₀, H)

    n_timesteps = 11
    u_max = zeros(n_timesteps)
    traction_magnitude = 3000 * range(0., 1.0, length=n_timesteps)
    V = 1
    VV = zeros(n_timesteps)

    # Create geometry, dofs and boundary conditions
    grid = generate_grid(Hexahedron, (1, 1, 1),
        Vec{3}([0., 0., 0.]), Vec{3}([1., 1., 1.]))
    addfacetset!(grid, "back_yz", x-> x[1] ≈ 0.)
    addfacetset!(grid, "back_xz", x-> x[2] ≈ 0.)
    addfacetset!(grid, "back_xy", x-> x[3] ≈ 0.)
    addfacetset!(grid, "front_xz", x-> x[2] ≈ 1.)
    addcellset!(grid, "front_xz", x-> x[2] ≈ 1.)
    addnodeset!(grid, "rve", x->x[1] >= 0.)
    addnodeset!(grid, "front_xz", x->x[2] ≈ 1.)
    gridnodeset = getnodeset(grid, "rve")
    refshape = RefHexahedron
    interpolation = Lagrange{refshape, 1}()^3

    # dh = create_dofhandler(grid, interpolation) # JuaFEM helper function
    dh = DofHandler(grid)
    add!(dh, :u, interpolation) # add a displacement field with 3 components
    close!(dh)
    # dbcs = create_bc(dh, grid) # create Dirichlet boundary-conditions
    dbcs = ConstraintHandler(dh)
    # Clamped on the left side
    add!(dbcs, Dirichlet(:u, getfacetset(grid, "back_xz"), (x,t) -> [0.0], [2]))
    add!(dbcs, Dirichlet(:u, getfacetset(grid, "back_yz"), (x,t) -> [0.0], [1]))
    add!(dbcs, Dirichlet(:u, getfacetset(grid, "back_xy"), (x,t) -> [0.0], [3]))
    close!(dbcs)
    vectorstrafe_i(i) = vectorstrafe(grid, dh, i)

    # cellvalues, facetvalues = create_values(interpolation)
    # setup quadrature rules
    qr      = QuadratureRule{refshape}(2)
    facet_qr = FacetQuadratureRule{refshape}(3)
    # cell and facetvalues for u
    cellvalues = CellValues(qr, interpolation)
    facetvalues = FacetValues(facet_qr, interpolation)

    # Pre-allocate solution vectors, etc.
    n_dofs = ndofs(dh)  # total number of dofs
    u  = zeros(n_dofs)  # solution vector
    Δu = zeros(n_dofs)  # displacement correction
    u₀ = zeros(n_dofs)
    r = zeros(n_dofs)   # residual
    K = allocate_matrix(dh); # tangent stiffness matrix

    # Create material states. One array for each cell, where each element is an array of material-
    # states - one for each integration point
    nqp = getnquadpoints(cellvalues)
    states = [MaterialState() for _ in 1:nqp, _ in 1:getncells(grid)]
    states_old = [MaterialState() for _ in 1:nqp, _ in 1:getncells(grid)]

    # pvd = paraview_collection("rve")
    pvd = paraview_collection("rve-quad")

    # Newton-Raphson loop
    NEWTON_TOL = 1 # 1 N
    print("\n Starting Netwon iterations:\n")
    gridnodes = deepcopy(getnodes(grid))
    gridnodes_range = gridnodes[[gridnodeset...]]
    gridnodes_range_view = @view gridnodes[[gridnodeset...]]
    Δt = 0.1 / (n_timesteps - 1)

    for (timestep, t) in enumerate(collect(range(0, 0.1, n_timesteps)))
        # t = timestep # actual time (used for evaluating d-bndc)
        traction = Vec((0., traction_magnitude[timestep], 0.))
        # newton_itr = -1
        @printf("\n Time step %d (t = %.4f s) @ %.4f m/s:\n", timestep - 1, t, V) # -1 to match ParaView
        i = 8; println((i, get_node_coordinate(gridnodes[i])[2], u[vectorstrafe_i(i)][2], u[vectorstrafe_i(i)][2] - u₀[vectorstrafe_i(i)][2]))
        VV[timestep] = V
        ΔV = 0.
        u₀ = u
        update!(dbcs, t) # evaluates the D-bndc at time t
        apply!(u, dbcs)  # set the prescribed values in the solution vector
        if t > 0
            # # apply_analytical!(u, dh, :u, x -> x + [0, V * Δt, 0], getcellset(grid, "front_xz"))
            # # [20250127T1655] (JMA3) - integration is needed for total displacement
            # println((V * Δt, integrate(VV[begin:timestep], collect(range(0, 0.1, n_timesteps))[begin:timestep])))
            # apply_analytical!(u, dh, :u, x -> [0., integrate(VV[begin:timestep], collect(range(0, 0.1, n_timesteps))[begin:timestep]), 0.], getcellset(grid, "front_xz"))
            dbcs_t = ConstraintHandler(dh)
            # add!(dbcs_t, Dirichlet(:u, getfacetset(grid, "front_xz"), (x,t) -> x + [0, V * Δt, 0], [1, 2, 3]))
            add!(dbcs_t, Dirichlet(:u, getfacetset(grid, "front_xz"), (x,t) -> [0., integrate(collect(range(0, 0.1, n_timesteps))[begin:timestep], VV[begin:timestep]), 0.], [1, 2, 3]))
            close!(dbcs_t)
            apply!(u, dbcs_t)  # set the prescribed values in the solution vector
        end

        i = 8; println((i, get_node_coordinate(gridnodes[i])[2], u[vectorstrafe_i(i)][2], u[vectorstrafe_i(i)][2] - u₀[vectorstrafe_i(i)][2]))

        # while true; newton_itr += 1
        newton_itr = 1
        while newton_itr <= 10
            # Tangent and residual contribution from the cells (volume integral)
            # doassemble!(K, r, cellvalues, dh, material, u, states, states_old);
            assembler = start_assemble(K, r)
            nu = getnbasefunctions(cellvalues)
            re = zeros(nu)     # element residual vector
            ke = zeros(nu, nu) # element tangent matrix

            for cell in CellIterator(dh)
                fill!(ke, 0)
                fill!(re, 0)
                eldofs = celldofs(cell)
                ue = u[eldofs]
                ue₀ = u₀[eldofs]
                state = @view states[:, cell.cellid]
                state_old = @view states_old[:, cell.cellid]
                # assemble_cell!(ke, re, cell, cellvalues, material, ue, state, state_old)
                n_basefuncs = getnbasefunctions(cellvalues)
                reinit!(cellvalues, cell)
                for q_point in 1:getnquadpoints(cellvalues)
                    # For each integration point, compute stress and material stiffness
                    ϵ = function_symmetric_gradient(cellvalues, q_point, ue) # Total strain
                    # ? [20250103T1654] (JMA3) - is σ SecondPiolaKirchoff?
                    σ, D, state[q_point] = compute_stress_tangent(ϵ, material, state_old[q_point])
                    dΩ = getdetJdV(cellvalues, q_point)
                    for i in 1:n_basefuncs
                        δϵ = shape_symmetric_gradient(cellvalues, q_point, i) # ? [20250128T1356] (JMA3) - is this the gradient matrix [B]
                        re[i] += 0.5 * (δϵ ⊡ σ) * dΩ # add internal force to residual
                        for j in 1:i # loop only over lower half
                            Δϵ = shape_symmetric_gradient(cellvalues, q_point, j) # ? [20250128T1356] (JMA3) - is this the gradient matrix [B]
                            ke[i, j] += 0.5 * δϵ ⊡ D ⊡ Δϵ * dΩ # ? [20250128T1356] (JMA3) - is this the correct form to assemble tangent stiffness matrix
                        end
                    end
                end
                symmetrize_lower!(ke)
                assemble!(assembler, eldofs, ke, re)
            end

            # Residual contribution from the Neumann boundary (surface integral)
            # doassemble_neumann!(r, dh, getfacetset(grid, "right"), facetvalues, traction)
            n_basefuncs = getnbasefunctions(facetvalues)
            re = zeros(n_basefuncs)                      # element residual vector
            for fc in FacetIterator(dh, getfacetset(grid, "front_xz"))
                # Add traction as a negative contribution to the element residual `re`:
                reinit!(facetvalues, fc)
                fill!(re, 0)
                for q_point in 1:getnquadpoints(facetvalues)
                    dΓ = getdetJdV(facetvalues, q_point)
                    for i in 1:n_basefuncs
                        δu = shape_value(facetvalues, q_point, i)
                        re[i] -= (δu ⋅ traction) * dΓ
                    end
                end
                assemble!(r, celldofs(fc), re)
            end
            norm_r = norm(r[Ferrite.free_dofs(dbcs)])

            print("Iteration: $newton_itr \tresidual: $(@sprintf("%.8f", norm_r))\n")
            if norm_r < NEWTON_TOL
                break
            end

            apply_zero!(K, r, dbcs)
            Δu = Symmetric(K) \ r
            u -= Δu
            newton_itr += 1
        end
        if newton_itr > 10
            error("Reached maximum Newton iterations, aborting")
        end

        # * [20250121T1533] (JMA3) - keep the `for`-loop below for properly strides through ⃗u to update node coordinates
        for (i, j, n) in zip(gridnodeset, 1:length(gridnodeset), gridnodes_range)
            gridnodes_range_view[j] = Node(Ferrite.Vec{3}(get_node_coordinate(n) + (u[vectorstrafe_i(i)]))) #  - u₀[vectorstrafe_i(i)])))
        end
        i = 8; println((i, get_node_coordinate(gridnodes[i])[2], u[vectorstrafe_i(i)][2], u[vectorstrafe_i(i)][2] - u₀[vectorstrafe_i(i)][2]))

        # Update the old states with the converged values for next timestep
        states_old .= states
        V += ΔV
        # u₀ = u

        u_max[timestep] = maximum(abs, u) # maximum displacement in current timestep
        # ## Postprocessing
        # Only a vtu-file corresponding to the last time-step is exported.
        #
        # The following is a quick (and dirty) way of extracting average cell data for export.
        strain = zeros(getncells(grid))
        strain_11 = zeros(getncells(grid))
        strain_22 = zeros(getncells(grid))
        strain_33 = zeros(getncells(grid))
        strain_12 = zeros(getncells(grid))
        strain_13 = zeros(getncells(grid))
        strain_23 = zeros(getncells(grid))
        strainᵖ = zeros(getncells(grid))
        mises_values = zeros(getncells(grid))
        κ_values = zeros(getncells(grid))
        for (el, cell_states) in enumerate(eachcol(states))
            for state in cell_states
                strain[el] += vonMises(state.ϵ)
                strain_11[el] += state.ϵ[1, 1]
                strain_22[el] += state.ϵ[2, 2]
                strain_33[el] += state.ϵ[3, 3]
                strain_12[el] += state.ϵ[1, 2]
                strain_13[el] += state.ϵ[1, 3]
                strain_23[el] += state.ϵ[2, 3]
                strainᵖ[el] += vonMises(state.ϵᵖ)
                mises_values[el] += vonMises(state.σ)
                κ_values[el] += state.k*material.H
            end
            strain[el] /= length(cell_states) # average von Mises stress
            strain_11[el] /= length(cell_states) # average von Mises stress
            strain_22[el] /= length(cell_states) # average von Mises stress
            strain_33[el] /= length(cell_states) # average von Mises stress
            strain_12[el] /= length(cell_states) # average von Mises stress
            strain_13[el] /= length(cell_states) # average von Mises stress
            strain_23[el] /= length(cell_states) # average von Mises stress
            strainᵖ[el] /= length(cell_states) # average von Mises stress
            mises_values[el] /= length(cell_states) # average von Mises stress
            κ_values[el] /= length(cell_states)     # average drag stress
        end
        VTKGridFile("rve-quad" * "-t$timestep", dh) do vtk
            write_solution(vtk, dh, u) # displacement field
            write_cell_data(vtk, strain, "strain [total]")
            write_cell_data(vtk, strain_11, "strain_11")
            write_cell_data(vtk, strain_22, "strain_22")
            write_cell_data(vtk, strain_33, "strain_33")
            write_cell_data(vtk, strain_12, "strain_12")
            write_cell_data(vtk, strain_13, "strain_13")
            write_cell_data(vtk, strain_23, "strain_23")
            write_cell_data(vtk, strainᵖ, "strain [PEEQ]")
            write_cell_data(vtk, mises_values, "von Mises [Pa]")
            write_cell_data(vtk, κ_values, "Drag stress [Pa]")
            collection_add_timestep(pvd, vtk, t)
        end
    end

    vtk_save(pvd);


    return u_max, traction_magnitude
end

u_max, traction_magnitude = solve();

using Plots
plot(
    vcat(0.0, u_max),                # add the origin as a point
    vcat(0.0, traction_magnitude),
    linewidth=2,
    title="Traction-displacement",
    label=nothing,
    markershape=:auto
    )
ylabel!("Traction [Pa]")
xlabel!("Maximum deflection [m]")

Plot

[KnutAM]

1. In re[i] -= (δu ⋅ traction) * dΓ, we have the nodal force re, and the shape function, δu (unitless), the traction (force/area), and dΓ area.
2. Stress or energy per volume: N/m^2 = (Nm)/m^3 = J/m^3
3. Not exactly sure what you mean with "integral of velocity through time", but if you want to control the displacement instead of the traction, you can use a Dirichlet boundary condition and specify a fixed displacement on some boundary.

was instrumental to strafe through the global displacement vector, u.

Why do you want the displacement of node 8 in the example? If you want to get the displacement at a specific coordinate, you can use PointEvalHandler.

[jmanthony3]

Thank you for clarifying the units of (1) and (2). What I meant in (3) is the total displacement by integrating the applied velocity from time zero up to the current time.

To address your question, I wanted to measure the actual displacement of a particular node while I have been learning how to use this package and what is going on after each timestep. For example, using the method of #788, I learned that the coordinate positions of the grid do not change during loadings. This led to me finding #1027 to explicitly transform the coordinates. I wasn't sure if PointEvalHandler was measuring displacement at a coordinate in the reference configuration or after all applied loadings. If the latter after all loading, then I would have no idea what coordinate position to provide to measure displacement.

[KnutAM]

I see, yes, unless you transform the grid, the grid remains at the reference coordinates (Ferrite doesn't know that you are simulating the displacements).
If you want to have the displacements of all nodes, you can use
https://ferrite-fem.github.io/Ferrite.jl/stable/reference/export/#Ferrite.evaluate_at_grid_nodes also.

The only time you need the dof-number AFAIK is when you want to apply a constraint to the node. This usually occurs when wanting to constrain a single node to avoid rigid body motion, in which case applying Dirichlet on a node set is supported.

[jmanthony3]

The clarification of units was still helpful. But how do I modify the code if I wanted to do a representative volume element (RVE) test to verify the behavior of the model against experimental data? I figured the most straightforward way was to apply the same displacements onto the RVE cube as the experimental data to see if the evaluated stress and strain data matched.

So far, this is the best I have to try and perform a RVE test against a fully annealed copper, ASTM-E8 specimen 3 mm thick pulled in tension at $10^{-3}$ mm/mm/s until fracture. Even if the shape of the von Mises plasticity material model doesn't resemble that of the experimental data, I would still have anticipated that the evaluated strain would be within the same range but seems to come up short.