Revision of Method of Manufactured Solutions from Sun, 2007-05-06 02:12

The revisions let you track differences between multiple versions of a post.

Nowadays, we often use the tools of computational analysis in engineering design. For the results of such analyses to be believable, the tools that are used have to pass rigorous tests. There are two categories of tests involved:

  1. Verification tests where one checks whether the mathematical model has been implemented according to design.
  2. Validation tests where one checks how closely the correctly implemented mathematical model mimics reality.

A method of verification of a numerical code that has received some attention is is the method of manufactured solutions. A systematic application of the method for complex nonlinear constitutive models in solid mechanics is yet to be found in the literature. The fluid mechanics community, however, uses the method on a regular basis.

Last year, at the World Conference on Computational Mechanics in Los Angeles, Prof. R. C. Batra told me that he had been using the method to verify his codes for a while. However, he had never felt the need to publish the approach because of its seemingly trivial nature. In this post I will give you an example of the method of fictitious body forces used by Prof. Batra [1,2].

The method of fictitious body forces involves assuming a displacement field over the body; finding the body force, initial conditions, and boundary conditions that fit that solution; and then using that body force and the computed boundary conditions and initial conditions in the numerical algorithm to arrive at a solution. If the computed solution matches the assumed displacement, then we're good to go.

The example below was meant to provide a starting point for a student. Please let me know if you find it useful or whether the equations are confusing rather that clarifying.

A 1-D example:

Consider a body in which the material point $ \ensuremath{\mathbf{X}}$ moves to $ \ensuremath{\mathbf{x}}$ during a motion. The motion is governed by the equations

$\displaystyle \begin{aligned} &\ensuremath{\ensuremath{\ensuremath{\boldsymbol{... ...ensuremath{\boldsymbol{\nabla}}_0 \ensuremath{\mathbf{u}}})^T \\ \end{aligned}$

where

$\displaystyle \ensuremath{\mathbf{x}}= \ensuremath{\boldsymbol{\varphi}}(\ensur... ...{\sigma}}= J^{-1} \ensuremath{\boldsymbol{P}} \ensuremath{\boldsymbol{F}}^T . $

If we now consider a one-dimensional bar of the material, then

$\displaystyle x_1 = \varphi(X_1, t) ;   F_{11} = 1 + \ensuremath{\frac{\partia... ...{\partial X_1}} ;   J = 1 + \ensuremath{\frac{\partial u_1}{\partial X_1}} . $

The left Cauchy-Green deformation tensor is

$\displaystyle B_{11} = 1 + 2\ensuremath{\frac{\partial u_1}{\partial X_1}} + \left(\ensuremath{\frac{\partial u_1}{\partial X_1}}\right)^2 . $

The Cauchy stress is given by (using a simplified constitutive model):

$\displaystyle \sigma_{11} = \ensuremath{\frac{1}{2}}  C (B_{11} - 1) = \ensure... ... X_1}} + \left(\ensuremath{\frac{\partial u_1}{\partial X_1}}\right)^2\right] $

where $ C$ is an elastic constant. Also, the Cauchy stress and the 1st Piola-Kirchhoff stress are related by

$\displaystyle \sigma_{11} = \frac{1}{F_{11}} P_{11} F_{11} = P_{11} . $

The one-dimensional momentum equation is

$\displaystyle \ensuremath{\frac{\partial P_{11}}{\partial X_1}} + \rho_0 b_1 = \rho_0 \ddot{u_1} . $

 

The method of manufactured solutions:

Let us try to apply the method of fictitious body forces to the bar problem.

The first step is to assume that the displacement field is given by (you can choose any other appropriate functional form)

$\displaystyle u_1(X_1, t) = X_1^2 \sin(\omega t) . $

Therefore,

$\displaystyle \ensuremath{\frac{\partial u_1}{\partial X_1}} = 2 X_1 \sin(\omeg... ...uremath{\frac{\partial^2 u_1}{\partial t^2}} = -X_1^2 \omega^2 \cos(\omega t) $

and the deformation gradient is

$\displaystyle F_{11} = 1 + 2 X_1 \sin(\omega t) . $

The left Cauchy-Green deformation is

$\displaystyle B_{11} = 1 + 4 X_1 \sin(\omega t) + 4 X_1^2 \sin^2(\omega t) . $

The first Piola-Kirchhoff stress is

$\displaystyle P_{11} = \ensuremath{\frac{1}{2}} C [4 X_1 \sin(\omega t) + 4 X_1^2 \sin^2(\omega t)] . $

Therefore,

$\displaystyle \ensuremath{\frac{\partial P_{11}}{\partial X_1}} = 2 C \sin(\omega t)[1 + 2 X_1 \sin(\omega t)] . $

Plugging $ \ensuremath{\frac{\partial P_{11}}{\partial X_1}}$ and $ \ensuremath{\frac{\partial^2 u_1}{\partial t^2}}$ into the momentum equation, we get

$\displaystyle 2 C \sin(\omega t)[1 + 2 X_1 \sin(\omega t)] + \rho_o b_1 = -\rho_0 X_1^2 \omega^2 \cos(\omega t) . $

Therefore, the body force is

$\displaystyle \boxed{ b_1 = - X_1^2 \omega^2 \cos(\omega t) - \cfrac{2 C}{\rho_0}  \sin(\omega t)[1 + 2 X_1 \sin(\omega t)] . } $

If the bar is of length $ L$ , the boundary conditions are

$\displaystyle \boxed{ u_1(0, t) = 0 \qquad\text{and}\qquad u_1(L, t) = L^2 \sin(\omega t) . } $

 

The initial conditions are

$\displaystyle \boxed{ u_1(X_1, 0) = X_1^2 \sin(0) = 0 . } $

When you apply the body force $ b_1(X_1, t)$ and the initial and boundary conditions, you should get a solution that matches the chosen function.

Note that all these have been done in a Lagrangian configuration. You could alternatively do the same assuming a function $ \ensuremath{\mathbf{u}}(\ensuremath{\mathbf{x}}, t)$ and transforming the equations accordingly.

 

References

 

1
R. C. Batra and X. Q. Liang.
Finite dynamic deformations of smart structures.
Computational Mechanics, 20:427-438, 1997.

 

2
R. C. Batra and B. M. Love.
Multiscale analysis of adiabatic shear bands in tungsten heavy alloy particulate composites.
Int. J. Multiscale Computational Engineering, 4(1):95-114, 2006.