Solid mechanics: Simulation Setup

We provide here an example of data and main files which are used to setup a simulation.

$\newcommand{\bm}[1]{\boldsymbol{\mathbf{#1}}} \newcommand{\bvec}[1]{\mathbf{#1}} \newcommand{\bx}{\mathbf{x}} \newcommand{\bd}{\mathbf{d}} \newcommand{\bv}{\dot{\mathbf{d}}} \newcommand{\ba}{\ddot{\mathbf{d}}} \newcommand{\bmw}{\mathbf{w}} \newcommand{\bmv}{\mathbf{v}} \newcommand{\Nh}[1][]{N_{h,#1}} \newcommand{\bU}{\mathbf{U}} \newcommand{\bp}{\bvec{p}} \newcommand{\VV}{\mathbf{V}} \newcommand{\WW}{\mathbf{W}} \newcommand{\QQ}{\mathbf{Q}} \renewcommand{\AA}{\mathbf{A}} \newcommand{\BB}{\mathbf{B}} \newcommand{\CC}{\mathbf{B}} \newcommand{\EE}{\mathbf{E}} \newcommand{\FF}{\mathbf{F}} \newcommand{\FFi}{\mathbf{f}^i} \newcommand{\FFe}{\mathbf{f}^e} \newcommand{\GG}{\mathbf{G}} \newcommand{\CC}{\mathbf{C}} \newcommand{\MM}{\mathbf{M}} \newcommand{\XX}{\mathbf{X}} \newcommand{\KK}{\mathbf{K}} \newcommand{\UU}{\mathbf{U}} \newcommand{\II}{\mathbf{I}} \renewcommand{\SS}{\mathbf{S}} \newcommand{\HH}{\mathbf{H}} \newcommand{\DD}{\mathbf{D}} \newcommand{\NN}{\mathbf{N}} \newcommand{\bmg}{\mathbf{g}} \newcommand{\R}{\mathbb{R}} \newcommand{\JJ}{\mathbf{J}} \newcommand{\RR}{\mathbf{R}} \newcommand{\PP}{\mathbf{P}}$

Datafile options

Define the 3 components of the volumetric forces (i.e. the right-hand side term $\mathbf{f}$ )

% Source term
data.force{1} = @(x, y, z, t, param)(0.*x.*y.*z);
data.force{2} = @(x, y, z, t, param)(0.*x.*y.*z);
data.force{3} = @(x, y, z, t, param)(0.*x.*y.*z);

Set the flags for the boundary conditions: on boundaries [2 3 4 5 6] we impose Dirichlet (essential) conditions on all the 3 components of the displacement. On boundary 200 we impose a normal pressure, i.e.

$\PP(\bd) \mathbf{n} = p \mathbf{n}.$

On boundary 210 we impose a Robin condition plus an external normal pressure, i.e.

$\PP(\bd) \mathbf{n} + \alpha \bd = p \mathbf{n}.$

data.flag_dirichlet{1} = [2 3 4 5 6];
data.flag_neumann{1}   = [];
data.flag_pressure{1}  = [200 210];
data.flag_robin{1}     = [210];

data.flag_dirichlet{2} = [2 3 4 5 6];
data.flag_neumann{2}   = [];
data.flag_pressure{2}  = [200 210];
data.flag_robin{2}     = [210];

data.flag_dirichlet{3} = [2 3 4 5 6];
data.flag_neumann{3}   = [];
data.flag_pressure{3}  = [200 210];
data.flag_robin{3}     = [210];

Specify the functions for Dirichlet, Neumann, Robin, and normal pressure boundary conditions

% Dirichlet
data.bcDir{1} = @(x, y, z, t, param)(0.*x.*y.*z);
data.bcDir{2} = @(x, y, z, t, param)(0.*x.*y.*z);
data.bcDir{3} = @(x, y, z, t, param)(0.*x.*y.*z);

% Neumann
data.bcNeu{1} = @(x, y, z, t, param)(0.*x.*y.*z);
data.bcNeu{2} = @(x, y, z, t, param)(0.*x.*y.*z);
data.bcNeu{3} = @(x, y, z, t, param)(0 + 0.*x.*y.*z);

% Normal Pressure
data.bcPrex{200}   = @(x, y, z, t, param)(-80*1333 + 0*x.*y.*z);%dyn/cm^2
data.bcPrex{210}   = @(x, y, z, t, param)(-3*1333 + 0*x.*y.*z);%dyn/cm^2

data.ElasticCoefRobin = 10^6;

Set material model and corresponding parameters. In particular for Linear, StVenantKirchhoff and NeoHookean material models, Young modulus and Poisson ratio should be provided. For the RaghavanVorp material model, the Alpha, Beta and Bulk parameters should be provided instead. Density is needed only for time-dependent simulations.

% material parameters
data.Material_Model = 'StVenantKirchhoff';% 'Linear','NeoHookean','RaghavanVorp'
data.Young   = 10^7;%dyn/cm^2
data.Poisson = 0.45;
%data.Alpha  = ... ;
%data.Beta   = ... ;
%data.Bulk   = ... ;

data.Density = 1;

Set tolerance and maximum number of iterations of the nonlinear solver (Newton)

% NonLinear Solver
data.NonLinearSolver.tol               = 1e-6;
data.NonLinearSolver.maxit             = 20;
data.NonLinearSolver.backtrackIter     = 4;
data.NonLinearSolver.backtrackFactor   = 0.7;

Set parameters for the linear solver and preconditioner. In this case, if a parallel pool is open we use gmres with a one-level Additive Schwarz preconditioner; otherwise we use a direct solver

% Solver and Preconditioner
%   If parallel pool available, use gmres with AdditiveSchwarz preconditioner
%   Otherwise use direct solver
poolobj = gcp('nocreate');
if isempty(poolobj)
    poolsize = 0;
else
    poolsize = poolobj.NumWorkers;
end

if poolsize > 0
    % Linear Solver
    data.LinearSolver.type              = 'gmres'; % MUMPS, backslash, gmres
    data.LinearSolver.tol               = 1e-7;
    data.LinearSolver.maxit             = 800;
    data.LinearSolver.gmres_verbosity   = 10;
    data.LinearSolver.mumps_reordering  = 7;
    
    % Preconditioner
    data.Preconditioner.type         = 'AdditiveSchwarz'; % AdditiveSchwarz, None, ILU
    data.Preconditioner.local_solver = 'MUMPS'; % matlab_lu, MUMPS
    data.Preconditioner.overlap_level     = 4;
    data.Preconditioner.mumps_reordering  = 7;
    data.Preconditioner.num_subdomains    = poolsize; %poolsize, number of subdomains
    
else
    
    % Linear Solver
    data.LinearSolver.type           = 'MUMPS'; % MUMPS, backslash, gmres
    % Preconditioner
    data.Preconditioner.type         = 'None'; % AdditiveSchwarz, None, ILU
end

Set time parameters for the generalized-$\alpha$ method

% Time options
data.time.t0         = 0;
data.time.dt         = 0.005; 
data.time.tf         = 5;
data.time.gamma      = 1/2;
data.time.beta       = 1/4;
data.time.alpha_m    = 0;
data.time.alpha_f    = 0;

Set output options. In this case we require to compute and export the Von Mises stress of the final solution

% OutPut Options
data.Output.ComputeVonMisesStress = true;

The Von Mises stress is defined as (in 3D)

$\sigma_{VM} = \sqrt{ \frac{1}{2} \left( (\sigma_{11} - \sigma_{22})^2 + (\sigma_{22} - \sigma_{33})^2 + (\sigma_{11} - \sigma_{33})^2 + 6 (\sigma_{23}^2 + \sigma_{31}^2 + \sigma_{12}^2)\right) }$

where $\bm{\sigma} = J^{-1} \PP \FF^T$ is the Cauchy stress tensor.