Ferrite 2D problem unphysical results

General Usage
, ,
1

Hello,
I am using the latest version of Ferrite (Ferrite v0.3.14). I’ve simplified my problem to a small 2D rectangular structure (height = 1; width = 2). The left side is clamped, and on the right-hand side, there is stress in the normal direction. Therefore, I would expect a symmetrical elongation, and not the result shown in the picture below.
image
Here is the running Ferrite code. I am thankful for any ideas.

using Ferrite, StaticArrays, SparseArrays, LinearAlgebra
using Plots
using Infiltrator
function compute_u(cellvalues, facevalues,dh,K,f,C,dbc)

    Ke = zeros(Float64,8,8)
    fe = zeros(Float64,8)

    assembler = start_assemble(K,f)

    for cell in CellIterator(dh)
        
        Ke .= 0.0
        fe .= 0.0
        reinit!(cellvalues, cell)

        for q_point in getnquadpoints(cellvalues)
            dΩ = getdetJdV(cellvalues, q_point)

            for i in 1:getnbasefunctions(cellvalues)
                ∇ε = tovoigt(SVector,
                    shape_symmetric_gradient(cellvalues, q_point, i),
                        offdiagscale=2.0)

                for j in 1:getnbasefunctions(cellvalues)
                    ∇δε = tovoigt(SVector,
                        shape_symmetric_gradient(cellvalues, q_point, j),
                        offdiagscale=2.0)

                    Ke[i,j] += ∇ε'*C*∇δε*dΩ
                end
            end
        end



        for face_id in 1:nfaces(cell)
            if onboundary(cell,face_id) && (cellid(cell),face_id) ∈ getfaceset(dh.grid,"pressure")
                reinit!(facevalues, cell,face_id)

                @show (cellid(cell),face_id)   
                for q_point in 1:getnquadpoints(facevalues)
                    dΓ = getdetJdV(facevalues, q_point)

                    for i in 1:getnbasefunctions(facevalues)
                        fe[i] += (facevalues.normals[1] ⋅ shape_value(facevalues, q_point, i)) * dΓ
                    end
                end
            end
        end
        assemble!(assembler,celldofs(cell),Ke,fe)   
    end
    apply!(K,f,dbc)
    u = Symmetric(K)\f
    return u
end


function FEM_MMWE()

    ## Create Domain ############
    corners = [Vec{2}((0.0, 0.0)),
               Vec{2}((2.0, 0.0)),
               Vec{2}((2.0, 1.0)),
               Vec{2}((0.0, 1.0))]

    n = 2
    grid = generate_grid(Quadrilateral, (n, n), corners);

    # node-/facesets for boundary conditions
    addnodeset!(grid, "clamped", x -> x[1] ≈ 0.0)
    addfaceset!(grid, "pressure", x -> x[1] ≈ 2.0);


    ## Define Interplotion and QuadratureRule for cells and domains
    qr      = QuadratureRule{2,RefCube}(3)
    cellvalues = CellVectorValues(qr, Lagrange{2,RefCube,1}(), Lagrange{2,RefCube,1}())
    qr_face = QuadratureRule{1,RefCube}(3)
    facevalues = FaceVectorValues(qr_face, Lagrange{2,RefCube,1}(), Lagrange{2,RefCube,1}())


    # DOF and ConstraintHandler
    dh = DofHandler(grid)
    add!(dh, :u, 2) # displacement
    close!(dh)

    dbc = ConstraintHandler(dh)
    add!(dbc, Dirichlet(:u, getnodeset(dh.grid, "clamped"), (x,t) -> zero(Vec{2}), [1,2]))
    close!(dbc)
    
    ## Define Material Properties
    E = 210e03; ν = 0.3 

    λ = E/(1-ν^2)
    C = Tensor{2,3,Float64}(λ*[1.0 ν 0.0; ν 1.0 0.0; 0.0 0.0 0.5*(1.0-ν)])

    K = create_sparsity_pattern(dh)
    f = zeros(Float64,size(K,1))

    ## Compute u
    u = compute_u(cellvalues,facevalues,dh,K,f,C,dbc) #compute_u(cellvalues, facevalues,dh,K,f,C,dbc)


    ### POSTPROCESSING ###
    u_nodes = reshape_to_nodes(dh,u,:u)
    u_1 = u_nodes[1,:]
    u_2 = u_nodes[2,:]

    u_coords = getcoordinates.(grid.nodes)
    x = getindex.(u_coords,1)
    y = getindex.(u_coords,2)


    # remove duplicated x coor
    x = unique(round.(x,digits = 5))
    y = unique(round.(y,digits = 5))

    display(x)
    display(y)

    X = [x_i for x_i in x, _ in y]
    Y = [y_i for _ in x, y_i in y]
    u_1 = reshape(u_1,length(x),length(y))
    u_2 = reshape(u_2,length(x),length(y))

    plotly()
    surface(X,Y,u_1,title = "x displacement")
    display(surface!())
    surface(X,Y,u_2,title = "y displacement")
    display(surface!())
end

FEM_MMWE()
2

Hey there,

the problem is that you are unintentionally integrating with one Gauss point the stiffness matrix, due to the missing 1: in 1:getnquadpoints(cellvalues)

17c16
<         for q_point in getnquadpoints(cellvalues)
---
>         for q_point in 1:getnquadpoints(cellvalues)

Further, I’d like to point out that the cast into tovoigt can be avoided and you can directly use the entities returned by shape_symmetric_gradient if you change the line Ke[i,j] += ∇ε'*C*∇δε*dΩ to Ke[i,j] += ∇ε ⊡ C ⊡ ∇δε*dΩ. If you want to use it this way, you should change C to a SymmetricTensor{4,2} and construct it as advertised in the docs of Tensors or Ferrite, see e.g., Demos · Tensors.jl

4 Likes
3

I would also highly recommend to either use FerriteViz.jl or ParaView for visualization and not trying to do it manually with Plots.jl.

Also note that the Ferrite devs are more active on Slack in #ferrite-fem as well as on Zulip in #Ferrite.jl .

2 Likes