modelling of general plate/shell

[duhd1993]

Hi, I'm trying to model a general Kirchhoff plate. The weak form looks like

By engineering convention, the moment and section force are commonly represented in matrix form

where the epsilon and kappa are vector form of the strain and curvature. All terms in ABD matrix could be populated so bending/in plane behavior is coupled.

Here is my seemingly working code.

mesh = MeshQuad1().refine(3)

e = ElementQuadBFS()
ev = ElementVector(e, 3)
ib = Basis(mesh, ev)

ABD = np.eye(6)

dofs = {
    'left' : ib.get_dofs(lambda x: x[0] == 0.0),
    'right': ib.get_dofs(lambda x: x[0] == 1.0),
}

@BilinearForm
def bilinf(u, v, w):
    def getStrain(u):
        kappa = -dd(u)
        eps = grad(u)

        strain = np.array([
            eps[0, 0],
            eps[1, 1],
            eps[0, 1] + eps[1, 0],
            kappa[2, 0, 0],
            kappa[2, 1, 1],
            kappa[2, 0, 1] + kappa[2, 1, 0]
        ])
        return strain
    strain_u = getStrain(u)
    strain_v = getStrain(v)

    return - dot(mul(ABD, strain_u),strain_v)

K = asm(bilinf, ib)

u = ib.zeros()
u[dofs['right'].nodal['u^1']] = 0.3

I = ib.complement_dofs(dofs)

u = solve(*condense(K, x=u, I=I))

I have two questions:

1. In order to solve this 3D problem instead of the original pure bending problem, I want to stack ElementQuadBFS with ElementVector. Am I doing it correctly? Besides, the in plane u1 and u2 requires less continuity. Maybe I should use ElementComposite and use normal element for u1 and u2?
2. I noticed that BilinearForm is best to be expressed in terms of tensors. When I slice trial/test functions, it becomes a numpy array, and grads no longer works. But it's more common to use matrix and vector form of stiffness and strain. I have to take grads of all components and then get the grads I need. This might cause unnecessary calculations. How should I approach this?

Any help is appreciated.

[kinnala]

Do you have a reference which has the governing equations for this formulation?

[duhd1993]

@kinnala It's more or less introduced in any composite materials textbook, but not easy to find this exact formulation. Check Mechanics of Laminated Composite Plates and Shells: Theory and Analysis by J. N. Reddy. Eq 3.3.19, neglecting all dynamic terms and complex boundary terms, becomes the form I showed. The constitutive equation is shown at 3.3.36 and 3.3.37.

The reason this form is not very commonly seen is in engineering practice people mostly use symmetric layups. Thus there is not tension-bending coupling and the B in the ABD matrix is zero. Therefore, the problem can be decomposed into a pure bending problem and 2D plane stress problem. But Abaqus has implementation for the general shell.

I was able to solve this problem using the code above. But as pointed out, it's not ideal to use C1 element for in plane behavior. I tried ElementComposite(ElementQuad0(), ElementQuad0(), ElementQuadBFS()), but it seems ElementQuadBFS is not compatible with ElementComposite.

[kinnala]

Thanks, I will check it out. By coincidence, Reddy was giving a series of lectures based on this book last week at our university. Unfortunately I wasn't able to attend due to my own lectures. I am not much of a plate theorist. I'm only familiar with the very basic Kirchhoff plate model which includes transverse deflection only.

I was not aware that ElementQuadBFS is incompatible with ElementComposite - perhaps it is a bug. If you have time, you could try to describe the issue that you are encountering and I can try to fix it.