Tutorial: Non-Linear dynamic analysis using Kirchhoff-Love shells

This example shows how to perform a non-linear dynamic analysis using the Dynamic solvers submodule of Structural Analysis Module. The dynamic solvers provided by this submodule work on user-defined operators from gsStructuralAnalysisOps, hence can be used with any element and in any application. In this example, we use the gsKLShell. The goals of the tutorial are the following:

• Solve a non-linear quasi-static problem using gsDynamicNewmark

For detailed information about the construction of gsThinShellAssembler and gsMaterialMatrixBase, please consult the tutorials Tutorial: Linear Kirchhoff-Love shell analysis and Tutorial: Non-Linear Kirchhoff-Love shell analysis.

Problem set-up

The set-up of the problem to be solved for this example is identical to the one in Tutorial: Non-Linear Kirchhoff-Love shell analysis. The reader is referred to that example for a full overview of the code.

Mass matrix assembly

The assembly of the linear stiffness matrix and the external force vector are similar as in Tutorial: Non-Linear Kirchhoff-Love shell analysis. To assemble the mass matrix, an additional assembly is performed using assembleMass(). The full assembly is given by:

    assembler.assemble();
    gsSparseMatrix<> K = assembler.matrix();
    gsVector<> F = assembler.rhs();
 
    assembler.assembleMass();
    gsSparseMatrix<> M = assembler.matrix();

Jacobian and Residual

To use the Newmark time integrator or any other time integrator from Dynamic solvers, we need to define a Jacobian matrix and a residual, optionally time-independent. Both are defined from the gsStructuralAnalysisOps, similar to Tutorial: Non-Linear Kirchhoff-Love shell analysis using the gsStructuralAnalysis module.

In addition, an operator for mass and damping needs to be defined. These operators are simple lambda functions returning a constant matrix:

    gsSparseMatrix<> C = gsSparseMatrix<>(assembler.numDofs(),assembler.numDofs());
    gsStructuralAnalysisOps<real_t>::Damping_t Damping = [&C](const gsVector<real_t> &, gsSparseMatrix<real_t> & m) { m = C; return true; };
    gsStructuralAnalysisOps<real_t>::Mass_t    Mass    = [&M](                          gsSparseMatrix<real_t> & m) { m = M; return true; };

Define the dynamic solver

Using the mass, damping and Jacobian matrices \(M\), \(C\) and \(K_{NL}(\mathbf{u})\) and a residual vector \(R(\mathbf{u})\), an non-linear dynamic solver can be defined. In the following, a gsDynamicNewmark is defined:

    gsInfo<<"Solving system with "<<assembler.numDofs()<<" DoFs\n";
    gsDynamicNewmark<real_t,true> solver(Mass,Damping,Jacobian,Residual);
    solver.options().setSwitch("Verbose",true);

The solution is initialized with zero-vectors, as follows:

    size_t N = assembler.numDofs();
    gsVector<> U(N), V(N), A(N);
    U.setZero();
    V.setZero();
    A.setZero();
    solver.setU(U);
    solver.setV(V);
    solver.setA(A);

Solve the non-linear dynamic problem

The non-linear dynamic problem is solved by stepping over time. After some initializations, a simple time stepping loop is used:

    index_t step = 50;
    gsParaviewCollection collection("Deformation");
    gsMultiPatch<> displ, def;
 
    real_t time = 0;
    real_t dt = 1e-2;
    solver.setTimeStep(dt);
    for (index_t k=0; k<step; k++)
    {
        gsInfo<<"Load step "<< k<<"\n";
        gsStatus status = solver.step();
        GISMO_ENSURE(status==gsStatus::Success,"Time integrator did not succeed");
 
        U = solver.solutionU();
        displ = assembler.constructDisplacement(U);
        def = assembler.constructSolution(U);
 
        // ! [Export visualization in ParaView]
        if (plot)
        {
            // Plot the displacements on the deformed geometry
            gsField<> solField(def, displ);
            std::string outputName = "Deformation" + std::to_string(k) + "_";
            gsWriteParaview<>( solField, outputName, 1000, true);
            collection.addPart(outputName + "0.vts",time);
        }
        time = solver.time();
    }

It can be seen that the export of the solution field is similar as in Tutorial: Linear Kirchhoff-Love shell analysis Tutorial: Non-Linear Kirchhoff-Love shell analysis, but the Paraview files are added to a gsParaviewCollection. This is saved at the end of the simulation as follows:

    if (plot)
        collection.save();

Annotated source file

Here is the full file nonlinear_shell_dynamic.cpp. Clicking on a function or class name will lead you to its reference documentation.

 
#include <gismo.h>
 
#ifdef gsKLShell_ENABLED
#include <gsKLShell/gsKLShell.h>
#endif
 
#include <gsStructuralAnalysis/src/gsDynamicSolvers/gsDynamicNewmark.h>
#include <gsStructuralAnalysis/src/gsStructuralAnalysisTools/gsStructuralAnalysisUtils.h>
 
 
using namespace gismo;
 
// Choose among various shell examples, default = Thin Plate
int main(int argc, char *argv[])
{
#ifdef gsKLShell_ENABLED
    bool plot  = false;
    index_t numRefine  = 1;
    index_t numElevate = 1;
 
    gsCmdLine cmd("Linear shell tutorial.");
    cmd.addInt( "e", "degreeElevation","Number of degree elevation steps to perform before solving (0: equalize degree in all directions)", numElevate );
    cmd.addInt( "r", "uniformRefine", "Number of Uniform h-refinement steps to perform before solving",  numRefine );
    cmd.addSwitch("plot", "Create a ParaView visualization file with the solution", plot);
    try { cmd.getValues(argc,argv); } catch (int rv) { return rv; }
 
    // Initialize [ori]ginal and [def]ormed geometry
    gsMultiPatch<> ori;
    gsReadFile<>("surfaces/paraboloid.xml", ori);
    ori.computeTopology();
 
    // p-refine
    if (numElevate!=0)
        ori.degreeElevate(numElevate);
    // h-refine
    for (int r =0; r < numRefine; ++r)
        ori.uniformRefine();
 
    gsMultiBasis<> bases(ori);
 
    // Define the boundary conditions object
    gsBoundaryConditions<> bc;
    // Set the geometry map for computation of the Dirichlet BCs
    bc.setGeoMap(ori);
 
    // Set the boundary conditions
    bc.addCornerValue(boundary::southwest, 0.0, 0, -1); // (corner,value, patch, unknown)
    bc.addCornerValue(boundary::southeast, 0.0, 0, -1); // (corner,value, patch, unknown)
    bc.addCornerValue(boundary::northwest, 0.0, 0, -1); // (corner,value, patch, unknown)
    bc.addCornerValue(boundary::northeast, 0.0, 0, -1); // (corner,value, patch, unknown)
 
    gsPointLoads<real_t> pLoads = gsPointLoads<real_t>();
    gsVector<> point(2);
    gsVector<> load (3);
    point<< 0.5, 0.5 ; load << 0, 0.0, -1e4 ;
    pLoads.addLoad(point, load, 0 );
 
    // The surface force is defined in the physical space, i.e. 3D
    gsFunctionExpr<> force("0","0","0",3);
 
 
    // Define material parameters
    // The material parameters are defined in the physical domain as well (!)
    gsConstantFunction<> E  (1.E9,3);
    gsConstantFunction<> nu (0.45,3);
    gsConstantFunction<> t  (1E-2,3);
    gsConstantFunction<> rho(1.E5,3);
 
    // Define a linear material, see \ref gsMaterialMatrixLinear.h
    // The first parameter is the physical domain dimension
 
    // gsMaterialMatrixLinear<3,real_t> materialMatrix(ori,t,E,nu);
 
    std::vector<gsFunctionSet<>*> parameters{&E,&nu};
    gsOptionList options;
    options.addInt("Material","Material model: (0): SvK | (1): NH | (2): NH_ext | (3): MR | (4): Ogden",1);
    options.addInt("Implementation","Implementation: (0): Composites | (1): Analytical | (2): Generalized | (3): Spectral",1);
    gsMaterialMatrixBase<real_t>::uPtr materialMatrix = getMaterialMatrix<3,real_t>(ori,t,parameters,rho,options);
 
    // Define the assembler.
    // The first parameter is the physical domain dimension
    // The second parameter is the real type
    // The third parameter controls the bending stiffness
    gsThinShellAssembler<3, real_t, true > assembler(ori,bases,bc,force,materialMatrix);
    assembler.setPointLoads(pLoads);
 
    assembler.assemble();
    gsSparseMatrix<> K = assembler.matrix();
    gsVector<> F = assembler.rhs();
 
    assembler.assembleMass();
    gsSparseMatrix<> M = assembler.matrix();
 
    // Function for the Jacobian
    gsStructuralAnalysisOps<real_t>::Jacobian_t Jacobian = [&assembler](gsVector<real_t> const &x, gsSparseMatrix<real_t> & m)
    {
        gsMultiPatch<> def;
        assembler.constructSolution(x,def);
        ThinShellAssemblerStatus status = assembler.assembleMatrix(def);
        m = assembler.matrix();
        return status==ThinShellAssemblerStatus::Success;
    };
    // Function for the Residual
    gsStructuralAnalysisOps<real_t>::Residual_t Residual = [&assembler](gsVector<real_t> const &x, gsVector<real_t> & v)
    {
        gsMultiPatch<> def;
        assembler.constructSolution(x,def);
        ThinShellAssemblerStatus status = assembler.assembleVector(def);
        v = assembler.rhs();
        return status==ThinShellAssemblerStatus::Success;
    };
 
    gsSparseMatrix<> C = gsSparseMatrix<>(assembler.numDofs(),assembler.numDofs());
    gsStructuralAnalysisOps<real_t>::Damping_t Damping = [&C](const gsVector<real_t> &, gsSparseMatrix<real_t> & m) { m = C; return true; };
    gsStructuralAnalysisOps<real_t>::Mass_t    Mass    = [&M](                          gsSparseMatrix<real_t> & m) { m = M; return true; };
 
 
    gsInfo<<"Solving system with "<<assembler.numDofs()<<" DoFs\n";
    gsDynamicNewmark<real_t,true> solver(Mass,Damping,Jacobian,Residual);
    solver.options().setSwitch("Verbose",true);
 
    size_t N = assembler.numDofs();
    gsVector<> U(N), V(N), A(N);
    U.setZero();
    V.setZero();
    A.setZero();
    solver.setU(U);
    solver.setV(V);
    solver.setA(A);
 
    index_t step = 50;
    gsParaviewCollection collection("Deformation");
    gsMultiPatch<> displ, def;
 
    real_t time = 0;
    real_t dt = 1e-2;
    solver.setTimeStep(dt);
    for (index_t k=0; k<step; k++)
    {
        gsInfo<<"Load step "<< k<<"\n";
        gsStatus status = solver.step();
        GISMO_ENSURE(status==gsStatus::Success,"Time integrator did not succeed");
 
        U = solver.solutionU();
        displ = assembler.constructDisplacement(U);
        def = assembler.constructSolution(U);
 
        // ! [Export visualization in ParaView]
        if (plot)
        {
            // Plot the displacements on the deformed geometry
            gsField<> solField(def, displ);
            std::string outputName = "Deformation" + std::to_string(k) + "_";
            gsWriteParaview<>( solField, outputName, 1000, true);
            collection.addPart(outputName + "0.vts",time);
        }
        time = solver.time();
    }
 
    // ![Save the paraview collection]
    if (plot)
        collection.save();
    // ![Save the paraview collection]
 
    return EXIT_SUCCESS;
#else
    GISMO_ERROR("The tutorial needs to be compiled with gsKLShell enabled");
    return EXIT_FAILED;
#endif
}// end main