Using SciPy and Other Third-party Python Packages
This example demonstrates using an optimization routine from scipy
to solve an inverse elastic problem.
A constant force, f, is applied to the tip of a cantilever beam. From a given beam geometry and displacement we use PFC to find the modulus of the beam. The beam is modeled with parallel bonded PFC particles.
The deflection is known analytically,
where f is the applied force, E is the Youngs Modulus, I is the moment of inertia and x is the length along the beam and L is the total length of the beam.
This is a simplification of a PFC verification problem see Cantilever Beam for more information.
import matplotlib
matplotlib.use("QT5Agg")
matplotlib.rcParams['backend.qt5']='PySide2'
import pylab as plt
import numpy as np
from scipy.constants import inch
import itasca as it
it.command("python-reset-state false")
from itasca import ballarray as ba
r = 8.0 # beam radius [inches]
l = 200.0 # beam length [inches]
x = np.linspace(0, l, 100)
E = 1e4 # [ksi]
I = 0.25*np.pi*r**4 # moment of inertia.
f = 1.0 # [kip] force applied in the positive y direction
First we import the required libraries and setup some parameters for the beam.
def analytical_deflection(x):
"Return the analytical solution for y deflection at the x point(s)"
return f * x**2 / 2.0 / E / I * ( l - x / 3.0)
def pfc_deflection(E_trial):
"""For a given parallel bond stiffness, return a tuple containing the
ball x positions and the y displacements.
calling this function cycles the PFC engine
"""
n = 11 # number of pfc particles
pb_kn = E_trial/l*(n-1)
print ("calculating deflection with pb_kn = {} ... ".format(pb_kn),)
it.command("""
model new
model large-strain on
model random
model domain ext -50 250 -50 50
contact cmat default model linearpbond
""")
for i in range(n):
it.ball.create(10.0, (i*20.0, 0.0, 0.0))
tip_ball = it.ball.find(n)
tip_ball.set_force_app_y(f)
it.command("""
ball attribute dens=1000 damp 0.7
ball property 'kn'=0 'ks'=0
model clean
contact method bond
contact property pb_ten 1e5 pb_coh 1e5 pb_kn {kn} pb_ks 188 pb_rmul 0.8
ball fix velocity-x velocity-y velocity-z ...
spin-x spin-y spin-z range id=1
model solve ratio-maximum 1e-4
""".format(kn=pb_kn))
pfc_x = ba.pos()[:,0]
pfc_dy = ba.disp()[:,1]
print ("done")
return pfc_x, pfc_dy
Two functions are defined one to give the analytical solution and one to give the PFC solution.
We can test everything is working by calling the pfc_deflection function with the known correct modulus. The pfc and analytical solution should match,
pfc_x, pfc_dy = pfc_deflection(E)
plt.plot(x, analytical_deflection(x))
plt.plot(pfc_x, pfc_dy, "o")
plt.xlabel("Length along beam [inches]")
plt.ylabel("Beam vertical displacement [inches]")
plt.draw()
plt.savefig("p3d-opt_comp.png")
- output:
calculating deflection with pb_kn = 500.0 ... done
The PFC beam deflection is shown as points. The PFC solution matches well against the analytical solution which is shown as a solid blue line.
Next we do the actual optimization calculation. The function
fminbound is imported from the scipy.optimize module. This
function takes three arguments, a callable function, an upper bound
and a lower bound. We define the function residual to return the
sum of the square of the difference between the PFC and analytical
solutions. We bound the optimization between 1e2 and 5.5e4.
- For more about
fminboundsee: http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.fminbound.html
from scipy.optimize import fminbound
def residual(pb_kn):
"""Return the sum of squares of the residual. This is a scalar measure
of the difference between the PFC solution and the analytical solution.
"""
pfc_x, pfc_dy = pfc_deflection(pb_kn)
a_dy = analytical_deflection(pfc_x)
plt.close('all')
plt.plot(x, analytical_deflection(x))
plt.plot(pfc_x, pfc_dy, "o")
plt.xlabel("Length along beam [inches]")
plt.ylabel("Beam vertical displacement [inches]")
plt.draw()
return ((a_dy - pfc_dy)**2).sum()
E_opt = fminbound(residual, 1e2, 5.5e4)
print ("""
E specified: {}
E calculated: {}""".format(E, E_opt))
plt.close('all')
- output:
calculating deflection with pb_kn = 1053.4967008815386 ... done calculating deflection with pb_kn = 1701.5032991184614 ... done calculating deflection with pb_kn = 653.0065982369227 ... done calculating deflection with pb_kn = 405.490102644616 ... done calculating deflection with pb_kn = 530.7543627593386 ... done calculating deflection with pb_kn = 529.6311596321268 ... done calculating deflection with pb_kn = 482.21349526223526 ... done calculating deflection with pb_kn = 501.8190835557665 ... done calculating deflection with pb_kn = 501.70523214550855 ... done calculating deflection with pb_kn = 499.7989842280888 ... done calculating deflection with pb_kn = 493.08192515191905 ... done calculating deflection with pb_kn = 500.68789623953967 ... done calculating deflection with pb_kn = 500.5925746896699 ... done calculating deflection with pb_kn = 500.66760871604833 ... done calculating deflection with pb_kn = 500.6732746302266 ... done calculating deflection with pb_kn = 500.6776216232916 ... done calculating deflection with pb_kn = 500.6759566385157 ... done calculating deflection with pb_kn = 500.681546177477 ... done calculating deflection with pb_kn = 500.67693683977063 ... done calculating deflection with pb_kn = 500.67762921620755 ... done calculating deflection with pb_kn = 500.6791253622799 ... done calculating deflection with pb_kn = 500.678200693155 ... done calculating deflection with pb_kn = 500.6778475009777 ... done calculating deflection with pb_kn = 500.6777125935705 ... done calculating deflection with pb_kn = 500.6776610635263 ... done calculating deflection with pb_kn = 500.6776413808009 ... done E specified: 10000.0 E calculated: 10013.55258432415
The calculated value for E and the known value are close.
Example File
The source code for the examples here is using_scipy.py.
| Was this helpful? ... | 3DEC © 2019, Itasca | Updated: Feb 25, 2024 |