Transcript Constant Gravity on a 1D neo

Verification of GIMP with
Manufactured Solutions
Philip Wallstedt – [email protected]
Jim Guilkey – [email protected]
Verification: Necessary Versus Sufficient
Eyeball Norms – no obvious error
• Not predictive: you already know the answer
Symmetry – some coding mistakes exposed
• Many mistakes are symmetric
Compare to existing code (Finite Element)
• Existing code solves different problems
• Existing code has unverified accuracy
• When differences are found, are they errors or not?
Experimental results – scattered data shows same trends
• Lack of data
• Differences don’t show what’s wrong with code
Known Solutions to PDE’s
• No solutions for large deformation
We need a better way: the Method of
Manufactured Solutions (MMS)
Recently proposed as ASME standard
“V&V 10 - 2006 Guide for Verification and Validation
in Computational Solid Mechanics”
Sufficient, not just necessary, if we test all modes:
• Boundary conditions
• Non-square cells and particles
• Time integration algorithms
• Shape functions
Each mode must be tested, but not all in the same test.
Once a mode has “passed”, then further testing not needed.
Rate of convergence is very sensitive to errors and
can be applied to individual pieces of a method
Displacement error compares current config to reference.
 u  x p  X p  u( X ) EXACT
Average error
Worst Error
  p 
 p

L1  max

 N 


L  max p 
Body Force on a 1D Bar
Body Force b(x,t)
Zero
Displacement
u(0,t) = 0
Zero Stress
(L,t) = 0
Body Force on a 1D Bar
Given
Momentum
Neo-Hookean Constitutive
Model
Constitutive Model with
assumptions: 1D, Poisson = 0
  σ  b  a
σ

J
ln JI 


FF
J
T
I
E
1
σ  F  
2
F
Find displacement u(x) – in general this cannot be done.

Start with the answer and reformulate backwards
Given Displacement
1D Neo-Hookean with
Poisson’s ratio = 0
Momentum
Solve for Gravity
u(X)
E
1
P  F  
2
F
  P  0b  0a
b a
Now we just take derivatives . . .
1
0
P
What answer (displacement field) do we start with?
The chosen displacement field(s) must:
• exercise all features of the code (large deformation,
translation, rotation, Dirichlet and Neumann boundaries
• be “smooth enough” – sufficiently differentiable in
time and space
• Conform to assumptions made by the method. For
GIMP this means zero normal stress at material
boundaries.
For the 1D rod a parabolic form should work:


u  c0  c1 X  c2 X 2 A(t )
Constants for the 1D bar
Zero displacement at X = 0

0.2

0  c0  c1 0  c2 02 A(t )
Scale displacement at L


A(t )  c0  c1L  c2 L A(t )
2
u( X )
0
0.2
0.1
0
0
0
0.5
1
X
L
Zero stress at X = L

E
1 E
1

P( L)  0   F    1  c1  2c2 L A(t ) 
2
F 2
1  c1  2c2 L A(t ) 
X (2 L  X )
u
A(t )
2
L
Choose a convenient time function A(t)
Trig function:
• Easy to differentiate
• Stays close to un-deformed shape
• Tests ability to preserve energy
• Can be made self-similar in time – same number of time
steps per period, regardless of problem stiffness.
But other functions work just fine, provided:
• problem always has sufficient particles per cell
• displacement field is well-behaved (for us A(t) > -1/2)
 E

X (2 L  X )


u
0
.
2
cos
π
t
2
 

L
 0

A detour and a review:
reference versus current configuration
Particles stationary in reference configuration
Grid stationary in current configuration
X=5
x=5
Reference
X = 4.75
x = 4.75
X = 4.75
x = 5.2
Current
X=?
x=5
Why manufacture solutions in the reference configuration?
– Because boundaries move in the current configuration.
How find the current length and apply boundary?
x(2 L0  x)
u ( x) 
A(t )
2
L0
( L0  L)( 2 L0  ( L0  L))
L  u ( L0  L) 
A(t )
2
L0
This is icky. We can avoid recursive / implicit definitions like
the above by using the reference configuration.
Reference Configuration vs Current Configuration
Momentum
Deformation
Gradient
Neo-Hookean
1D,
Poisson = 0
Reference Configuration
“Total Lagrange”
Current Configuration
“Updated Lagrange”
  P  0b  0a
  σ  b  a
F(x)  I  u(x)
1
F(X)  I  u(X)

P   ln JF  F FF  I
1
1
T
E
1 

P   F(X) 
2
F(X) 
Stress Transformation:

σ

J
ln JI 


FF
J
T
I
E
1 

σ   F(x) 
2
F(x) 
1
P  JF  σ

Return to the 1D Bar: Take Derivatives
Given
Displacement
Deformation Gradient
Divergence of Stress
Solve for b(X)
X (2 L  X )
u
A(t )
2
L
2( L  X )
F  1
A
2
L
2
E   2( L  X )  
  P   2 1  1 
A A
2

L  
L
 
1
b 2
L

  2( L  X )  2  
E
  1  1 
A  A
 X (2 L  X ) A
2

0  
L
  

1D Bar: Restate the Problem
Body Force
1
b 2
L
Zero
Displacement
u(0,t) = 0

  2( L  X )  2  
E
  1  1 
A  A
 X (2 L  X ) A
2

0  
L
  

The answer is:
X (2 L  X )
u( X ) 
A(t )
2
L
Zero Stress
(L,t) = 0
 E 
Solve with GIMP where A(t )  0.2 cos
π t 
 0 
Now we can measure convergence under large
deformation – the kind of problem MPM/GIMP
is designed to solve
max error
0.1
0.01
0.001
0.0001
GIMP
UGIMP
0.00001
0.000001
10
100 mesh size
1000
Skeptical Questions
1st Piola-Kirchoff is neither objective nor fully
Lagrangian – doesn’t that cause problems?
MPM is a first-order, fully non-linear method. It can’t
be expected to agree with your manufactured solution
due to its non-linearity.