Nonlinear thermoelastic plate bending

This example solves a cantilever rectangular plate with a uniformly distributed load on the top of the plate.

For direct solver, run with options: -ksp_type preonly -pc_type lu

Initialize libMesh library.

libMesh::LibMeshInit init(argc, argv);

Create Mesh object on default MPI communicator and generate a 2D mesh of QUAD4 elements. We discretize with 12 elements in the x-direction (0.0 to 36.0 inches) and 40 elements in the y-direction (0.0 to 120.0 inches). Note that in libMesh, all meshes are parallel by default in the sense that the equations on the mesh are solved in parallel by PETSc. A "ReplicatedMesh" is one where all MPI processes have the full mesh in memory, as compared to a "DistributedMesh" where the mesh is "chunked" up and distributed across processes, each having their own piece.

Real
length = 0.3,
width  = 0.3;
libMesh::ReplicatedMesh mesh(init.comm());
libMesh::MeshTools::Generation::build_square(mesh, 30, 30, 0.0, length, 0.0, width, libMesh::QUAD9);
mesh.print_info();


libMesh::Point
pt(length/2., width/2., 0.), // location of mid-point before shift
pt0,
dr1, dr2;

const libMesh::Node
*nd = nullptr;

if a finite radius is defined, change the mesh to a circular arc of specified radius

libMesh::MeshBase::node_iterator
n_it   = mesh.nodes_begin(),
n_end  = mesh.nodes_end();

initialize the pointer to a node

nd   = *n_it;
pt0  = *nd;

for ( ; n_it != n_end; n_it++) {
    
    dr1  = pt0;
    dr1 -= pt;
    
    dr2  = **n_it;
    dr2 -= pt;
    
    if (dr2.norm() < dr1.norm()) {
        
        nd  = *n_it;
        pt0 = *nd;
    }
}

std::cout << *nd << std::endl;

For later reference, the boundary and subdomain ID's generated by the libMesh mesh generation are sketched below.

*         y               (#) Boundary ID 
*         | (2)           [#] Subdomain ID 
*   120.0 O-----O 
*         |     |          x - right 
*         |     |          y - up 
*     (3) | [0] | (1)      z - out of screen 
*         |     | 
*         |     | 
*     0.0 O-----O--x 
*           (0) 
*        0.0   36.0 
*  

Create EquationSystems object, which is a container for multiple systems of equations that are defined on a given mesh.

libMesh::EquationSystems equation_systems(mesh);

Add system of type MAST::NonlinearSystem (which wraps libMesh::NonlinearImplicitSystem) to the EquationSystems container. We name the system "structural" and also get a reference to the system so we can easily reference it later.

MAST::NonlinearSystem & system = equation_systems.add_system<MAST::NonlinearSystem>("structural");

Create a finite element type for the system. Here we use first order Lagrangian-type finite elements.

libMesh::FEType fetype(libMesh::FIRST, libMesh::LAGRANGE);

Initialize the system to the correct set of variables for a structural analysis. In libMesh this is analogous to adding variables (each with specific finite element type/order to the system for a particular system of equations.

MAST::StructuralSystemInitialization structural_system(system,
                                                       system.name(),
                                                       fetype);

Initialize a new structural discipline using equation_systems.

MAST::PhysicsDisciplineBase discipline(equation_systems);

Create and add boundary conditions to the structural system. A Dirichlet BC adds fixed displacement BCs. Here we use the side boundary ID numbering created by the libMesh generator to clamp the edge of the mesh along x=0.0. We apply the BC to all variables on each node in the subdomain ID 0, which clamps this edge.

MAST::DirichletBoundaryCondition clamped_edge0, clamped_edge1, clamped_edge2, clamped_edge3;    // Create BC object.
std::vector<unsigned int> vars = {0, 1, 2};
clamped_edge0.init(0, vars);   // Assign boundary ID and variables to constrain
clamped_edge1.init(1, vars);   // Assign boundary ID and variables to constrain
clamped_edge2.init(2, vars);   // Assign boundary ID and variables to constrain
clamped_edge3.init(3, vars);   // Assign boundary ID and variables to constrain
discipline.add_dirichlet_bc(0, clamped_edge0);     // Attach boundary condition to discipline
discipline.add_dirichlet_bc(1, clamped_edge1);     // Attach boundary condition to discipline
discipline.add_dirichlet_bc(2, clamped_edge2);     // Attach boundary condition to discipline
discipline.add_dirichlet_bc(3, clamped_edge3);     // Attach boundary condition to discipline
discipline.init_system_dirichlet_bc(system);      // Initialize the BC in the system.

Initialize the equation system since we now know the size of our system matrices (based on mesh, element type, variables in the structural_system) as well as the setup of dirichlet boundary conditions. This initialization process is basically a pre-processing step to preallocate storage and spread it across processors.

equation_systems.init();
equation_systems.print_info();

Create parameters.

MAST::Parameter thickness("th",   0.001);
MAST::Parameter E("E",            72.e9);
MAST::Parameter nu("nu",           0.33);
MAST::Parameter rho("rho",        2700.);
MAST::Parameter kappa("kappa",    5./6.);
MAST::Parameter alpha("alpha",   1.5e-5);
MAST::Parameter zero("zero",        0.0);
MAST::Parameter pressure("p",       1.0);
MAST::Parameter temperature("T",    0.0);

Create ConstantFieldFunctions used to spread parameters throughout the model.

MAST::ConstantFieldFunction th_f("h", thickness);
MAST::ConstantFieldFunction E_f("E", E);
MAST::ConstantFieldFunction nu_f("nu", nu);
MAST::ConstantFieldFunction rho_f("rho", rho);
MAST::ConstantFieldFunction kappa_f("kappa", kappa);
MAST::ConstantFieldFunction alpha_f("alpha_expansion", alpha);
MAST::ConstantFieldFunction off_f("off", zero);
MAST::ConstantFieldFunction pressure_f("pressure", pressure);
MAST::ConstantFieldFunction temperature_f("temperature", temperature);
MAST::ConstantFieldFunction ref_temp_f("ref_temperature", zero);

Initialize load.

MAST::BoundaryConditionBase surface_pressure(MAST::SURFACE_PRESSURE);
surface_pressure.add(pressure_f);
discipline.add_volume_load(0, surface_pressure);

MAST::Examples::ThermalJacobianScaling jac_scaling;
jac_scaling.set_enable(true);
MAST::BoundaryConditionBase temperature_load(MAST::TEMPERATURE);
temperature_load.add(temperature_f);
temperature_load.add(ref_temp_f);

temperature_load.add(jac_scaling);

discipline.add_volume_load(0, temperature_load);

Create the material property card ("card" is NASTRAN lingo) and the relevant parameters to it. An isotropic material needs elastic modulus (E) and Poisson ratio (nu) to describe its behavior.

MAST::IsotropicMaterialPropertyCard material;
material.add(E_f);
material.add(nu_f);
material.add(kappa_f);
material.add(alpha_f);
material.add(rho_f);

Create the section property card. Attach all property values.

MAST::Solid2DSectionElementPropertyCard section;
section.add(th_f);
section.add(off_f);
section.set_strain(MAST::NONLINEAR_STRAIN);

Attach material to the card.

section.set_material(material);

Initialize the specify the subdomain in the mesh that it applies to.

discipline.set_property_for_subdomain(0, section);

Create nonlinear assembly object and set the discipline and structural_system. Create reference to system.

MAST::NonlinearImplicitAssembly assembly;
MAST::StructuralNonlinearAssemblyElemOperations elem_ops;
assembly.set_discipline_and_system(discipline, structural_system);
elem_ops.set_discipline_and_system(discipline, structural_system);
MAST::NonlinearSystem& nonlinear_system = assembly.system();
jac_scaling.set_assembly(assembly);

Zero the solution before solving.

nonlinear_system.solution->zero();

output the solution to exo file.

libMesh::ExodusII_IO output(mesh);

write the header to the load.txt file

std::ofstream out;
if (mesh.comm().rank() == 0) {
    out.open("load.txt", std::ofstream::out);
    out
    << std::setw(10) << "iter"
    << std::setw(25) << "temperature"
    << std::setw(25) << "pressure"
    << std::setw(25) << "displ" << std::endl;
}

Solve the system

unsigned int
n_steps  = 800,
dof_num  = nd->dof_number(0, 2, 0);
std::vector<Real> vec1;
std::vector<unsigned int> vec2 = {dof_num};

Real
temp_max = 10.;

for (unsigned int i=0; i<=n_steps; i++) {
    
    temperature = temp_max * (i*1.)/(1.*n_steps);
    libMesh::out << "load step: " << i << "  temp: " << temperature() << std::endl;
    
    nonlinear_system.solve(elem_ops, assembly);

get the value of the node at the center of the plate for output

system.solution->localize(vec1, vec2);

write the value to the load.txt file

    if (mesh.comm().rank() == 0) {
        out
        << std::setw(10) << i
        << std::setw(25) << temperature()
        << std::setw(25) << pressure()
        << std::setw(25) << vec1[0] << std::endl;
    }

    assembly.reset_residual_norm_history();
    output.write_timestep("panel.exo", equation_systems, i+1, system.time);
    system.time += 1.;
}