Surface coupling with oomph-lib

Tutorial

Before we can do coupled simulations, see How to run oomph-lib codes in MercuryDPM to learn how to run oomph-lib simulations in the framework of MercuryDPM.

Now you can include oomph-lib source files into your Driver codes. A good example of this is SolidBeamDemo.cpp:

#include "Species/LinearViscoelasticSlidingFrictionSpecies.h"
#include "Oomph/SCoupling/SCoupledSolidProblem.h"
 
int main()
{
    SCoupledSolidProblem<RefineableQDPVDElement<3, 2>> problem;
    //set name
    problem.setName("CoupledBeamUnitTest");
 
    //setup mercury
    problem.setDomain({0,0,0},{0.16,0.04,0.01});
    // add species
    auto species = problem.speciesHandler.copyAndAddObject(LinearViscoelasticSpecies());
    species->setDensity(1000);
    species->setStiffness(5.2);
    // add particle
    SphericalParticle p(species);
    p.setRadius(5e-3);
    p.setPosition({problem.getXCenter(),problem.getYCenter(),problem.getZMax()+p.getRadius()});
    p.setVelocity({0,0,-0.05});
    auto particle = problem.particleHandler.copyAndAddObject(p);
    // set time step
    double tc = species->getCollisionTime(species->getMassFromRadius(p.getRadius()));
    problem.setTimeStep(0.02*tc);
 
    // setup oomph
    problem.setElasticModulus(1e6);
    problem.setDensity(2500);
    problem.setSolidCubicMesh(7, 5, 5, 0, 0.16, 0, 0.04, 0, 0.01);
    problem.pinBoundaries({problem.Boundary::X_MIN, problem.Boundary::X_MAX});
    problem.setNewtonSolverTolerance(3e-8);
    problem.prepareForSolve();
    problem.coupleBoundary(problem.Boundary::Z_MAX);
    // sync time steps
    problem.setOomphTimeStep(problem.getTimeStep());
 
    // setup time
    problem.setTimeMax(100*problem.getTimeStep());
 
    // Solve the  problem and save mesh
    problem.solveSurfaceCoupling();
    problem.saveSolidMesh();
    return 0;
}

This code includes the class SCoupledSolidProblem, which is derived from the SolidProblem class. This class is defined in the MercuryDPM kernel, and is specifically designed to solve problems with coupled elastic bodies. As you can see, the driver code defines a class of type SCoupledSolidProblem, with the element type RefineableQDPVDElement<3,2>, sets the properties of the elastic solid, then calls solveSteady to solve the linear algebra problem.

Further test cases

We have implemented the following test cases:

SolidOnParticleBed.cpp

Illustrates the surface coupling technique implemented in MercuryDPM and oomph-lib. A solid cube (10cm side length) rests on four particles (2cm diameter), which rests on a fixed bottom wall. Gravity acts on both particles and the solid cube. The particle stiffness is chosen very low such that the compression due to gravitational load can be seen: the equilibrium overlap of the particle-wall and particle-cube contacts is 1mm, and you can see in the simulation that the overlaps oscillate between 0 and 2 mm. Both solid cube and particle contacts are non-dissipative, so oscillations are not damped out.

CoupledBeam.cpp

This is the demo from the coupling paper by Hongyang: A 20m x 80cm x 80cm beam fixed on the left end, with two particles of 20 cm diameter rolling down the beam at 0.1 m/s velocity, under gravity. Both solid cube and particle contacts are non-dissipative, so oscillations are not damped out. (Does not fully work yet; the smooth wall implementation is not there)

CoupledBeamUnitTest.cpp

A short simulation showing the unsteady interaction of a particle impacting a beam that is fixed on both sides. Colors indicate the displacement of the beam, due to gravity and the particle-beam interaction.