The azimuthally Fourier-decomposed equations of 3D time-harmonic linear elasticity

The aim of this tutorial is to demonstrate the solution of the time-harmonic equations of linear elasticity in cylindrical polar coordinates, using a Fourier decomposition of the solution in the azimuthal direction. These equations are useful to describe forced, time-harmonic, non-axisymmetric oscillations of axisymmetric elastic bodies.

 Acknowledgement: This implementation of the equations and the documentation were developed jointly with Robert Harter (Thales Underwater Systems Ltd) with financial support from a KTA Secondment grant from University of Manchester's EPSRC-funded Knowledge Transfer Account.

# Theory

Consider a three-dimensional, axisymmetric body (of density , Young's modulus , and Poisson's ratio ), occupying the region whose boundary is . Assuming that the body performs time-harmonic oscillations of frequency of , we use cylindrical coordinates . The equations of time-harmonic linear elasticity can then be written as where , and the stresses, body force and displacements are given by , and respectively. Note that here and henceforth, the superscript asterisk notation is used to distinguish dimensional quantities from their non-dimensional counterparts where required. (The coordinate is by definition dimensionless, and so we do not use an asterisk when referencing this parameter).

The body is subject to imposed time-harmonic displacements along , and subject to an imposed traction along where so that where is the outer unit normal on the boundary.

The stresses and displacements are related by the constitutive equations where represents the transpose of . Note that in cylindrical coordinates, the second-order tensor is given by and is equal to the trace of this matrix.

We non-dimensionalise the equations, using a problem specific reference length, , and a timescale , and use Young's modulus to non-dimensionalise the body force and the stress tensor:  The non-dimensional form of the linear elasticity equations is then given by where , and the non-dimensional parameter is the ratio of the elastic body's intrinsic timescale, , to the problem-specific timescale, , that we used to non-dimensionalise time. The boundary conditions are Given the assumed axisymmetry of the body we expand all quantities in a Fourier series in the azimuthal coordinate by writing, This decomposition allows us to remove the -dependence from the equations by writing , where represents any physical parameter in the problem. Furthermore, since the governing equations are linear, we can solve for each Fourier component separately and simply specify the Fourier wavenumber as a parameter.

# Implementation

Within oomph-lib, the non-dimensional version of the two-dimensional Fourier-decomposed equations (1) with the constitutive equations (2) are implemented in the TimeHarmonicFourierDecomposedLinearElasticityEquations equations class. Following our usual approach, discussed in the (Not-So-)Quick Guide, this equation class is then combined with a geometric finite element to form a fully-functional finite element. For instance, the combination of the TimeHarmonicFourierDecomposedLinearElasticityEquations class with the geometric finite element QElement<2,3> yields a nine-node quadrilateral element. As usual, the mapping between local and global (Eulerian) coordinates within an element is given by, where the coordinates are enumerated as . is the number of nodes in the element, is the -th global (Eulerian) coordinate (enumerated as above) of the -th Node in the element, and the are the element's shape functions, defined in the geometric finite element.

We allow for the presence of damping by allowing the constitutive parameters and forcing frequency to be complex-valued. The three components of the displacement field therefore have real and imaginary parts and we store the six real-valued nodal unknowns in the order and use the shape functions to interpolate the displacements as where is the -th displacement component (enumerated as indicated above) at the -th Node in the element.

# The test problem

The governing equations are fairly complicated and it is difficult to come up with non-trivial analytical solutions that could be used to validate the implementation. We therefore construct an analytical solution by postulating a displacement field and providing a body force that makes this a solution of the equations.

Specifically we consider the time-harmonic non-axisymmetric deformation of an annular elastic body that occupies the region .

The displacement field is an exact solution of the governing equations if the body is subject to a body force where and are the non-dimensional Lamé parameters (non-dimensionalised on ). The body is subject to a non-zero traction on all four boundaries; for example, on the inner boundary (where ) the traction is We choose to set this traction as a boundary condition, whilst pinning the displacements on the remaining boundaries where we impose a prescribed displacement according to (3).

# Results

The figures below show plots of and for a Fourier wavenumber of and geometric parameters . We set , corresponding to an exponentially growing time-periodic forcing; , corresponding to a slightly dissipative material (see Comments ); and . The imaginary part of the solution is small (though not identically equal to zero) but it converges to zero under mesh refinement; see Exercises . Computed (red) and exact (green) solution for real part of the radial displacement component. Computed (red) and exact (green) solution for real part of the axial displacement component. Computed (red) and exact (green) solution for real part of the azimuthal displacement component.

# Global parameters and functions

As usual, we define all non-dimensional parameters in a namespace. In this namespace, we also define the (Fourier-decomposed) body force, the traction to be applied on boundary 3, and the exact solution. Note that, consistent with the enumeration of the unknowns, discussed above, the order of the components in the functions that specify the body force and the surface traction is .

//===start_of_namespace=================================================
/// Namespace for global parameters
//======================================================================
{
/// Define Poisson's ratio Nu
std::complex<double> Nu(0.3,0.05);
/// Define the non-dimensional Young's modulus
std::complex<double> E(1.0,0.01);
// Lame parameters
std::complex<double> lambda = E*Nu/(1.0+Nu)/(1.0-2.0*Nu);
std::complex<double> mu = E/2.0/(1.0+Nu);
/// Define Fourier wavenumber
/// Define the non-dimensional square angular frequency of
/// time-harmonic motion
std::complex<double> Omega_sq (10.0,5.0);
/// Length of domain in r direction
double Lr = 1.0;
/// Length of domain in z-direction
double Lz = 2.0;
// Set up min & max (r,z) coordinates
double rmin = 0.1;
double zmin = 0.3;
double rmax = rmin+Lr;
double zmax = zmin+Lz;
/// Define the imaginary unit
const std::complex<double> I(0.0,1.0);
/// The traction function at r=rmin: (t_r, t_z, t_theta)
void boundary_traction(const Vector<double> &x,
const Vector<double> &n,
Vector<std::complex<double> > &result)
{
result = -6.0*pow(x,2)*mu*cos(x)-
lambda*(I*double(Fourier_wavenumber)*pow(x,2)*pow(x,3)+
(4.0*pow(x,2)+pow(x,3))*cos(x));
result = -mu*(3.0*pow(x,2)-pow(x,3))*sin(x);
result = -mu*pow(x,2)*(2*pow(x,3)+I*double(Fourier_wavenumber)*
cos(x));
}
/// The body force function; returns vector of complex doubles
/// in the order (b_r, b_z, b_theta)
void body_force(const Vector<double> &x,
Vector<std::complex<double> > &result)
{
result =
x*(-2.0*I*lambda*double(Fourier_wavenumber)*pow(x,3)-cos(x)*
(lambda*(8.0+3.0*x)-
mu*(pow(double(Fourier_wavenumber),2)
-16.0+x*(x-3.0))+pow(x,2)*Omega_sq));
result =
x*sin(x)*(mu*(pow(double(Fourier_wavenumber),2)-9.0)+
4.0*x*(lambda+mu)+pow(x,2)*
(lambda+2.0*mu-Omega_sq))-
3.0*I*double(Fourier_wavenumber)*pow(x,2)*pow(x,2)*(lambda+mu);
result =
-x*(8.0*mu*pow(x,3)-pow(double(Fourier_wavenumber),2)*pow(x,3)*
(lambda+2.0*mu)+pow(x,2)*(pow(x,3)*Omega_sq+6.0*mu*x)+
I*cos(x)*double(Fourier_wavenumber)*
(lambda*(4.0+x)+mu*(6.0+x)));
}
/// The exact solution in a flat-packed vector:
// 0: u_r[real], 1: u_z[real],..., 5: u_theta[imag]
void exact_solution(const Vector<double> &x,
Vector<double> &u)
{
u = pow(x,3)*cos(x);
u = pow(x,3)*sin(x);
u = pow(x,3)*pow(x,3);
u = 0.0;
u = 0.0;
u = 0.0;
}
} // end_of_namespace
Namespace for global parameters.
Definition: cylinder.cc:43
double Lz
Length of domain in z-direction.
Definition: cylinder.cc:65
const std::complex< double > I(0.0, 1.0)
Define the imaginary unit.
std::complex< double > lambda
Definition: cylinder.cc:51
double Lr
Length of domain in r direction.
Definition: cylinder.cc:62
std::complex< double > mu
Definition: cylinder.cc:52
void boundary_traction(const Vector< double > &x, const Vector< double > &n, Vector< std::complex< double > > &result)
The traction function at r=rmin: (t_r, t_z, t_theta)
Definition: cylinder.cc:77
std::complex< double > Nu(0.3, 0.05)
Define Poisson's ratio Nu.
void exact_solution(const Vector< double > &x, Vector< double > &u)
The exact solution in a flat-packed vector:
Definition: cylinder.cc:114
void body_force(const Vector< double > &x, Vector< std::complex< double > > &result)
The body force function; returns vector of complex doubles in the order (b_r, b_z,...
Definition: cylinder.cc:92
std::complex< double > Omega_sq(10.0, 5.0)
Define the non-dimensional square angular frequency of time-harmonic motion.
std::complex< double > E(1.0, 0.01)
Define the non-dimensional Young's modulus.
int Fourier_wavenumber
Define Fourier wavenumber.
Definition: cylinder.cc:55

# The driver code

We start by setting the number of elements in each of the two coordinate directions before creating a DocInfo object to store the output directory.

//===start_of_main======================================================
/// Driver code
//======================================================================
int main(int argc, char* argv[])
{
// Number of elements in r-direction
unsigned nr=5;
// Number of elements in z-direction (for (approximately) square elements)
unsigned nz=unsigned(double(nr)*Global_Parameters::Lz/Global_Parameters::Lr);
// Set up doc info
DocInfo doc_info;
// Set output directory
doc_info.set_directory("RESLT");
int main(int argc, char *argv[])
Driver code.
Definition: cylinder.cc:360

We build the problem using two-dimensional QTimeHarmonicFourierDecomposedLinearElasticityElements, solve using the Problem::newton_solve() function, and document the results.

// Set up problem
<QTimeHarmonicFourierDecomposedLinearElasticityElement<3> >
// Solve
problem.newton_solve();
// Output the solution
problem.doc_solution(doc_info);
} // end_of_main
Class to validate time harmonic linear elasticity (Fourier decomposed)
Definition: cylinder.cc:134

# The problem class

The Problem class is very simple. As in other problems with Neumann boundary conditions, we provide separate meshes for the "bulk" elements and the face elements that apply the traction boundary conditions. The latter are attached to the relevant faces of the bulk elements by the function assign_traction_elements().

//===start_of_problem_class=============================================
/// Class to validate time harmonic linear elasticity (Fourier
/// decomposed)
//======================================================================
template<class ELEMENT>
{
public:
/// Constructor: Pass number of elements in r and z directions
/// and boundary locations
const unsigned &nr, const unsigned &nz,
const double &rmin, const double& rmax,
const double &zmin, const double& zmax);
/// Update before solve is empty
/// Update after solve is empty
/// Doc the solution
void doc_solution(DocInfo& doc_info);
private:
/// Allocate traction elements on the bottom surface
/// Pointer to the bulk mesh
Mesh* Bulk_mesh_pt;
/// Pointer to the mesh of traction elements
}; // end_of_problem_class
FourierDecomposedTimeHarmonicLinearElasticityProblem(const unsigned &nr, const unsigned &nz, const double &rmin, const double &rmax, const double &zmin, const double &zmax)
Constructor: Pass number of elements in r and z directions and boundary locations.
Definition: cylinder.cc:174
Mesh * Surface_mesh_pt
Pointer to the mesh of traction elements.
Definition: cylinder.cc:163
void assign_traction_elements()
Allocate traction elements on the bottom surface.
Definition: cylinder.cc:289
Mesh * Bulk_mesh_pt
Pointer to the bulk mesh.
Definition: cylinder.cc:160
void actions_after_newton_solve()
Update after solve is empty.
Definition: cylinder.cc:149
void actions_before_newton_solve()
Update before solve is empty.
Definition: cylinder.cc:146
void doc_solution(DocInfo &doc_info)
Doc the solution.
Definition: cylinder.cc:318

# The problem constructor

We begin by building the meshes and pin the displacements on the appropriate boundaries. Recall that the order of the six real unknowns stored at the nodes is //===start_of_constructor=============================================
/// Problem constructor: Pass number of elements in coordinate
/// directions and size of domain.
//====================================================================
template<class ELEMENT>
(const unsigned &nr, const unsigned &nz,
const double &rmin, const double& rmax,
const double &zmin, const double& zmax)
{
//Now create the mesh
//Create the surface mesh of traction elements
Surface_mesh_pt=new Mesh;
assign_traction_elements();
// Set the boundary conditions for this problem: All nodes are
// free by default -- just pin & set the ones that have Dirichlet
// conditions here
// storage for nodal position
Vector<double> x(2);
// Storage for prescribed displacements
Vector<double> u(6);
// Now set displacements on boundaries 0 (z=zmin),
//------------------------------------------------
// 1 (r=rmax) and 2 (z=zmax)
//--------------------------
for (unsigned ibound=0;ibound<=2;ibound++)
{
unsigned num_nod=Bulk_mesh_pt->nboundary_node(ibound);
for (unsigned inod=0;inod<num_nod;inod++)
{
// Get pointer to node
Node* nod_pt=Bulk_mesh_pt->boundary_node_pt(ibound,inod);
// get r and z coordinates
x=nod_pt->x(0);
x=nod_pt->x(1);
// Pinned in r, z and theta
nod_pt->pin(0);nod_pt->pin(1);nod_pt->pin(2);
nod_pt->pin(3);nod_pt->pin(4);nod_pt->pin(5);
// Compute the value of the exact solution at the nodal point
Vector<double> u(6);
// Set the displacements
nod_pt->set_value(0,u);
nod_pt->set_value(1,u);
nod_pt->set_value(2,u);
nod_pt->set_value(3,u);
nod_pt->set_value(4,u);
nod_pt->set_value(5,u);
}
} // end_of_loop_over_boundary_nodes

Next we loop over the bulk mesh elements and assign the constitutive parameters, the body force, the Fourier wavenumber and the non-dimensional frequency to each element.

// Complete the problem setup to make the elements fully functional
// Loop over the elements
unsigned n_el = Bulk_mesh_pt->nelement();
for(unsigned e=0;e<n_el;e++)
{
// Cast to a bulk element
ELEMENT *el_pt = dynamic_cast<ELEMENT*>(Bulk_mesh_pt->element_pt(e));
// Set the body force
el_pt->body_force_fct_pt() = &Global_Parameters::body_force;
// Set the pointer to Poisson's ratio
el_pt->nu_pt() = &Global_Parameters::Nu;
// Set the pointer to Fourier wavenumber
el_pt->fourier_wavenumber_pt() = &Global_Parameters::Fourier_wavenumber;
// Set the pointer to non-dim Young's modulus
el_pt->youngs_modulus_pt() = &Global_Parameters::E;
// Set the pointer to square of the angular frequency
el_pt->omega_sq_pt() = &Global_Parameters::Omega_sq;
}// end loop over elements

We then loop over the traction elements and specify the applied traction.

// Loop over the traction elements
unsigned n_traction = Surface_mesh_pt->nelement();
for(unsigned e=0;e<n_traction;e++)
{
// Cast to a surface element
TimeHarmonicFourierDecomposedLinearElasticityTractionElement<ELEMENT>*
el_pt =
dynamic_cast<TimeHarmonicFourierDecomposedLinearElasticityTractionElement
<ELEMENT>* >(Surface_mesh_pt->element_pt(e));
// Set the applied traction
el_pt->traction_fct_pt() = &Global_Parameters::boundary_traction;
}// end loop over traction elements

The two sub-meshes are now added to the problem and a global mesh is constructed before the equation numbering scheme is set up, using the function assign_eqn_numbers().

// Add the submeshes to the problem
// Now build the global mesh
build_global_mesh();
// Assign equation numbers
cout << assign_eqn_numbers() << " equations assigned" << std::endl;
} // end of constructor

# The traction elements

We create the face elements that apply the traction to the boundary .

//===start_of_traction===============================================
/// Make traction elements along the boundary r=rmin
//===================================================================
template<class ELEMENT>
{
unsigned bound, n_neigh;
// How many bulk elements are next to boundary 3
bound=3;
n_neigh = Bulk_mesh_pt->nboundary_element(bound);
// Now loop over bulk elements and create the face elements
for(unsigned n=0;n<n_neigh;n++)
{
// Create the face element
FiniteElement *traction_element_pt
= new TimeHarmonicFourierDecomposedLinearElasticityTractionElement<ELEMENT>
(Bulk_mesh_pt->boundary_element_pt(bound,n),
Bulk_mesh_pt->face_index_at_boundary(bound,n));
}
} // end of assign_traction_elements

# Post-processing

As expected, this member function documents the computed solution.

//==start_of_doc_solution=================================================
/// Doc the solution
//========================================================================
template<class ELEMENT>
doc_solution(DocInfo& doc_info)
{
ofstream some_file;
char filename;
// Number of plot points
unsigned npts=5;
// Output solution
sprintf(filename,"%s/soln.dat",doc_info.directory().c_str());
some_file.open(filename);
Bulk_mesh_pt->output(some_file,npts);
some_file.close();
// Output exact solution
sprintf(filename,"%s/exact_soln.dat",doc_info.directory().c_str());
some_file.open(filename);
Bulk_mesh_pt->output_fct(some_file,npts,
some_file.close();
// Doc error
double error=0.0;
double norm=0.0;
sprintf(filename,"%s/error.dat",doc_info.directory().c_str());
some_file.open(filename);
Bulk_mesh_pt->compute_error(some_file,
error,norm);
some_file.close();
// Doc error norm:
cout << "\nNorm of error: " << sqrt(error) << std::endl;
cout << "Norm of solution: " << sqrt(norm) << std::endl << std::endl;
cout << std::endl;
} // end_of_doc_solution

• Given that we non-dimensionalised all stresses on Young's modulus it seems odd that we provide the option to specify a non-dimensional Young's modulus via the member function TimeHarmonicFourierDecomposedLinearElasticityEquations::youngs_modulus_pt(). The explanation for this is that this function specifies the ratio of the material's actual Young's modulus to the Young's modulus used in the non-dimensionalisation of the equations. The capability to specify such ratios is important in problems where the elastic body is made of multiple materials with different constitutive properties. If the body is made of a single, homogeneous material, the specification of the non-dimensional Young's modulus is not required – it defaults to 1.0. In the example considered above, the specification of the non-dimensional Young's modulus as creates a small amount of damping in the material whose actual stiffness is still characterised by the (real-valued and dimensional) Young's modulus used to non-dimensionalise the equations.
• Note that we also allow Poisson's ratio (whose specification is required) to be complex-valued. We are not aware of any meaningful physical interpretation of non-real Poisson ratios but provide this option because it appears to allow a better characterisation of some materials.

## Exercises

• Confirm that the specification of Poisson's ratio is required: What happens if you comment out its assignment in the problem constructor?
• Confirm that the small imaginary part of the computed displacement field for the test problem goes to zero under mesh refinement.
• Change the problem setup to the (less contrived) case where the deformation of the cylinder is driven by a time-periodic pressure load acting on the inside while its upper and lower ends are held at a fixed position. (You can cheat – there's another tutorial that shows you how to do it...).