Calculate total potential energy?

1

Hello

I am working on a simple 2-D linear elastic problem, after solving for displacement u, I need to calculated the total potential energy.

I wanted to obtain this scalar Pi from the Fenics solution u and edge traction T

Pi = 0.5 * inner(epsilon(u), sigma(u)) - inner(T, u)

where sigma and epsilon are defined as :

def epsilon(u): return sym(grad(u))
def sigma(u): return lam*tr(epsilon(u)) * Identity(2) + 2.0*mu*epsilon(u)

I have tried a few ways to get this Pi but not successful yet, could anyone offer some solutions? My Fenics version is 2019.1.0.

Thank you!

2

Hi,
you need to add the integration measures dx and ds in your formulation and then use assemble to obtain the corresponding total energy.

3

Thank you very much for your advice!

I did this:

pi = assemble( 0.5 * inner(epsilon(u), sigma(u)) *dx - inner(T, u)*dss(1) )

where dss(1) is previously define that the one edge that has applied traction.

I did obtain a scalar value for pi, however it is extremely small, down to machine epsilon. I wonder if this value is not the total potential energy, but the variation of the total potential energy ( =0 ). What’s your thought on this?

4

This is my full code, it is just the classic plate with a hole and traction, using quarter symmetry and plane stress. The for loop is just to vary the grid size.

``
import matplotlib.pyplot as plt
from dolfin import *
import numpy as np
from ufl import nabla_div
from mshr import *

def epsilon(u):
return sym(grad(u))
def sigma(u):
return lamtr(epsilon(u)) * Identity(2) + 2.0mu*epsilon(u)

def left(x): return x[0] < DOLFIN_EPS
def bottom(x): return x[1] < DOLFIN_EPS

def right(self,x,on_boundary): return x[0] > 2-DOLFIN_EPS
class Trac(SubDomain):
def inside(self,x,on_boundary):
return near(x[0], 2.0)

p=0.1
E=2.2
nu=1.0/10.0
mu=E/(2*(1+nu))

lam = Enu/((1+nu)(1-2nu)) ## plane strain
lam = 2mulam/(lam+2mu)

geo = Rectangle(dolfin.Point(0., 0.), dolfin.Point(2., 1.))
- Circle(dolfin.Point(0.0, 0.0), 4.0/7.0)

for grid in range(2,50):

mesh=generate_mesh(geo, grid )

V = VectorFunctionSpace(mesh, 'P', 1)
zero = Constant(0.0)
bc1 = DirichletBC(V.sub(0), zero, left)    
bc2 = DirichletBC(V.sub(1), zero, bottom)
bc = [bc1, bc2]


du = TrialFunction(V)
v = TestFunction(V)

a = inner(sigma(du), epsilon(v))*dx
boundary_subdomains = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundary_subdomains.set_all(3)  
trac=Trac()
trac.mark(boundary_subdomains, 1)
dss=Measure('ds', domain=mesh, subdomain_data=boundary_subdomains)

T = Constant( (p, 0) ) # surface traction
L = inner(T, v)*dss(1)

# Compute solution
u = Function(V)
solve(a == L, u, bc)

a = inner(sigma(u), epsilon(u))*dx
L = inner(T, u)*dss(1)

pi = assemble(a-L)

print(pi)

``

5

Hi,
you’re missing the 1/2 factor in front of a for the potential energy