Plane strain pull-in instability#
Code for plane strain large deformation electro-elasticity with u-p formulation
Electro-elastic pull-in instability of a plane strain VHB block
Units#
Length: mm
Mass: kg
Time: s
Charge: nC
Force: mN
Stress: kPa
Electric potential: kV
Software:#
Dolfinx v0.8.0
In the collection “Example Codes for Coupled Theories in Solid Mechanics,”
By Eric M. Stewart, Shawn A. Chester, and Lallit Anand.
Import modules#
# Import FEnicSx/dolfinx
import dolfinx
# For numerical arrays
import numpy as np
# For MPI-based parallelization
from mpi4py import MPI
comm = MPI.COMM_WORLD
rank = comm.Get_rank()
# PETSc solvers
from petsc4py import PETSc
# specific functions from dolfinx modules
from dolfinx import fem, mesh, io, plot, log
from dolfinx.fem import (Constant, dirichletbc, Function, functionspace, Expression )
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from dolfinx.io import VTXWriter, XDMFFile
# specific functions from ufl modules
import ufl
from ufl import (TestFunctions, TrialFunction, Identity, grad, det, div, dev, inv, tr, sqrt, conditional ,\
gt, dx, inner, derivative, dot, ln, split, exp, eq, cos, acos, ge, le, outer)
# basix finite elements
import basix
from basix.ufl import element, mixed_element, quadrature_element
# Matplotlib for plotting
import matplotlib.pyplot as plt
plt.close('all')
# For timing the code
from datetime import datetime
# Set level of detail for log messages (integer)
# Guide:
# CRITICAL = 50, // errors that may lead to data corruption
# ERROR = 40, // things that HAVE gone wrong
# WARNING = 30, // things that MAY go wrong later
# INFO = 20, // information of general interest (includes solver info)
# PROGRESS = 16, // what's happening (broadly)
# TRACE = 13, // what's happening (in detail)
# DBG = 10 // sundry
#
log.set_log_level(log.LogLevel.WARNING)
Define geometry#
# An axisymetric box
length = 10.0 # mm
# domain = mesh.create_rectangle(MPI.COMM_WORLD, [[0.0,0.0], [length,length]],\
# [10,10], mesh.CellType.quadrilateral)
domain = mesh.create_rectangle(MPI.COMM_WORLD, [[0.0,0.0], [length,length]],\
[10,10], mesh.CellType.triangle, diagonal=mesh.DiagonalType.crossed)
x = ufl.SpatialCoordinate(domain)
# Identify the planar boundaries of the box mesh
#
def Left(x):
return np.isclose(x[0], 0)
def Right(x):
return np.isclose(x[0], length)
def Bottom(x):
return np.isclose(x[1], 0)
def Top(x):
return np.isclose(x[1], length)
# Mark the sub-domains
#
boundaries = [(1,Left),(2,Bottom),(3,Right),(4,Top)]
# build collections of facets on each subdomain and mark them appropriately.
#
facet_indices, facet_markers = [], [] # initalize empty collections of indices and markers.
fdim = domain.topology.dim - 1 # geometric dimension of the facet (mesh dimension - 1)
for (marker, locator) in boundaries:
facets = mesh.locate_entities(domain, fdim, locator) # an array of all the facets in a
# given subdomain ("locator")
facet_indices.append(facets) # add these facets to the collection.
facet_markers.append(np.full_like(facets, marker)) # mark them with the appropriate index.
# Format the facet indices and markers as required for use in dolfinx.
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
#
# Add these marked facets as "mesh tags" for later use in BCs.
facet_tags = mesh.meshtags(domain, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
Print out the unique facet index numbers
# top_imap = domain.topology.index_map(2) # index map of 2D entities in domain (facets)
# values = np.zeros(top_imap.size_global) # an array of zeros of the same size as number of 2D entities
# values[facet_tags.indices]=facet_tags.values # populating the array with facet tag index numbers
print(np.unique(facet_tags.values)) # printing the unique indices
[1 2 3 4]
Visualize reference configuration and boundary facets
import pyvista
pyvista.set_jupyter_backend('html')
from dolfinx.plot import vtk_mesh
pyvista.start_xvfb()
# initialize a plotter
plotter = pyvista.Plotter()
# Add the mesh.
topology, cell_types, geometry = plot.vtk_mesh(domain, domain.topology.dim)
grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)
plotter.add_mesh(grid, show_edges=True)
plotter.view_xy()
labels = dict( xlabel='X', ylabel='Y')
plotter.add_axes(**labels)
plotter.screenshot("results/pe_pullin_mesh.png")
from IPython.display import Image
Image(filename='results/pe_pullin_mesh.png')
Un-comment this cell to see an interactive visualization of the mesh#
# plotter.show()
Define boundary and volume integration measure#
# Define the boundary integration measure "ds" using the facet tags,
# also specify the number of surface quadrature points.
ds = ufl.Measure('ds', domain=domain, subdomain_data=facet_tags, metadata={'quadrature_degree':2})
# Define the volume integration measure "dx"
# also specify the number of volume quadrature points.
dx = ufl.Measure('dx', domain=domain, metadata={'quadrature_degree': 2})
# Create facet to cell connectivity required to determine boundary facets.
domain.topology.create_connectivity(domain.topology.dim, domain.topology.dim)
domain.topology.create_connectivity(domain.topology.dim, domain.topology.dim-1)
domain.topology.create_connectivity(domain.topology.dim-1, domain.topology.dim)
# Define facet normal
n2D = ufl.FacetNormal(domain)
n = ufl.as_vector([n2D[0], n2D[1], 0.0]) # define n as a 3D vector for later use
Material parameters#
# Mechanical parameters
Geq_0 = Constant(domain, 15.0) # Shear modulus, kPa
Kbulk = Constant(domain, PETSc.ScalarType(1.0e3*Geq_0)) # Bulk modulus, kPa
I_m = Constant(domain, 175.0) # Gent locking paramter
# Electrostatic parameters
vareps_0 = Constant(domain, 8.85E-3) # permittivity of free space pF/mm
vareps_r = Constant(domain, 5.0) # relative permittivity, dimensionless
vareps = vareps_r*vareps_0 # permittivity of the material
Function spaces#
U2 = element("Lagrange", domain.basix_cell(), 2, shape=(2,)) # For displacement
P1 = element("Lagrange", domain.basix_cell(), 1) # For pressure and electric potential
#
TH = mixed_element([U2, P1, P1]) # Taylor-Hood style mixed element
ME = functionspace(domain, TH) # Total space for all DOFs
#
V1 = functionspace(domain, P1) # Scalar function space.
V2 = functionspace(domain, U2) # Vector function space
#
# Define actual functions with the required DOFs
w = Function(ME)
u, p, phi = split(w) # displacement u, presssure p, and electric potential phi
# A copy of functions to store values in the previous step
w_old = Function(ME)
u_old, p_old, phi_old = split(w_old)
# Define test functions
u_test, p_test, phi_test = TestFunctions(ME)
# Define trial functions needed for automatic differentiation
dw = TrialFunction(ME)
Initial conditions#
The initial conditions for degrees of freedom \(\mathbf{u}\), \(p\), and \(\phi\) are zero everywhere
These are imposed automatically, since we have not specified any non-zero initial conditions.
Subroutines for kinematics and constitutive equations#
#-------------------------------------------------------------
# Utility subroutines
#-------------------------------------------------------------
# Subroutine for a "safer" sqrt() function which avoids a divide by zero
# when differentiated.
def safe_sqrt(x):
return sqrt(x + 1.0e-16)
# Gradient of vector field u for plane strain
def pe_grad_vector(u):
grad_u = grad(u)
pe_grad_u = ufl.as_tensor([ [grad_u[0,0], grad_u[0,1], 0.0],
[grad_u[1,0], grad_u[1,1], 0.0],
[ 0.0, 0.0, 0.0] ])
return pe_grad_u
# Gradient of scalar field y
# (just need an extra zero for dimensions to work out)
def pe_grad_scalar(y):
grad_y = grad(y)
pe_grad_y = ufl.as_vector([grad_y[0], grad_y[1], 0.])
return pe_grad_y
#-------------------------------------------------------------
# Subroutines for kinematics
#-------------------------------------------------------------
# Axisymmetric deformation gradient
def F_pe_calc(u):
dim = len(u) # dimension of problem (2)
Id = Identity(dim) # 2D Identity tensor
F = Id + grad(u) # 2D Deformation gradient
F_pe = ufl.as_tensor([ [F[0,0], F[0,1], 0.0],
[F[1,0], F[1,1], 0.0],
[ 0.0, 0.0, 1.0]]) # Full plane strain F
return F_pe
#-------------------------------------------------------------
# Subroutines for calculating the electric field and displacement
#-------------------------------------------------------------
# Referential electric displacement
def Dmat_calc(u, phi):
F = F_pe_calc(u)
J = det(F)
C = F.T*F
e_R = - pe_grad_scalar(phi) # referential electric field
Dmat = vareps * J* inv(C)*e_R
return Dmat
#-------------------------------------------------------------
# Subroutines for calculating the stress
#-------------------------------------------------------------
# Generalized shear modulus for Gent model
def Geq_Gent_calc(u):
F = F_pe_calc(u)
J = det(F)
C = F.T*F
Cdis = J**(-2/3)*C
I1 = tr(Cdis)
z = I1-3
z = conditional( gt(z, I_m), 0.95*I_m, z ) # Keep from blowing up
Geq_Gent = Geq_0 * (I_m/(I_m - z))
return Geq_Gent
# Mechanical Cauchy stress for Gent material
def T_mech_calc(u,p):
Id = Identity(3)
F = F_pe_calc(u)
J = det(F)
B = F*F.T
Bdis = J**(-2/3)*B
Geq = Geq_Gent_calc(u)
T_mech = (1/J)* Geq * dev(Bdis) - p * Id
return T_mech
# Maxwell contribution to the Cauchy stress
def T_maxw_calc(u,phi):
F = F_pe_calc(u)
e_R = - pe_grad_scalar(phi) # referential electric field
e_sp = inv(F.T)*e_R # spatial electric field
# Spatial Maxwel stress
T_maxw = vareps*(outer(e_sp,e_sp) - 1/2*(inner(e_sp,e_sp))*Identity(3))
return T_maxw
# Piola stress
def T_mat_calc(u, p, phi):
Id = Identity(3)
F = F_pe_calc(u)
J = det(F)
#
T_mech = T_mech_calc(u,p)
#
T_maxw = T_maxw_calc(u,phi)
#
T = T_mech + T_maxw
#
Tmat = J * T * inv(F.T)
return Tmat
Evaluate kinematics and constitutive relations#
# Kinematical quantities
F = F_pe_calc(u)
J = det(F)
C = F.T*F
Fdis = J**(-1/3)*F
Cdis = J**(-2/3)*C
I1 = tr(Cdis)
# Mechanical Cauchy stress
T_mech = T_mech_calc(u, p)
# Electrostatic Cauchy stress
T_maxw =T_maxw_calc(u, phi)
# Piola stress
Piola = T_mat_calc(u, p, phi)
# Referential electric displacement
Dmat = Dmat_calc(u, phi)
Weak forms#
# The weak form for the equilibrium equation
#
Res_1 = inner( Piola, pe_grad_vector(u_test))*dx
# The auxiliary equation for the pressure
#
Res_2 = dot((p/Kbulk + ln(J)/J) , p_test)*dx
# The weak form for Gauss's equation
Res_3 = inner(Dmat, pe_grad_scalar(phi_test))*dx
# The total residual
Res = Res_1 + Res_2 + Res_3
# Automatic differentiation tangent:
a = derivative(Res, w, dw)
Set-up output files#
# results file name
results_name = "pe_pullin"
# Function space for projection of results
P1 = element("Lagrange", domain.basix_cell(), 1)
VV1 = fem.functionspace(domain, P1) # linear scalar function space
#
U1 = element("Lagrange", domain.basix_cell(), 1, shape=(2,))
VV2 = fem.functionspace(domain, U1) # linear Vector function space
#
T1 = element("Lagrange", domain.basix_cell(), 1, shape=(3,3))
VV3 = fem.functionspace(domain, T1) # linear tensor function space
# For visualization purposes, we need to re-project the stress tensor onto a linear function space before
# we write it (and its components and the von Mises stress, etc) to the VTX file.
#
# This is because the stress is a complicated "mixed" function of the (quadratic Lagrangian) displacements
# and the (quadrature representation) plastic strain tensor and scalar equivalent plastic strain.
#
# First, define a function for setting up this kind of projection problem for visualization purposes:
def setup_projection(u, V):
trial = ufl.TrialFunction(V)
test = ufl.TestFunction(V)
a = ufl.inner(trial, test)*dx
L = ufl.inner(u, test)*dx
projection_problem = dolfinx.fem.petsc.LinearProblem(a, L, [], \
petsc_options={"ksp_type": "cg", "ksp_rtol": 1e-16, "ksp_atol": 1e-16, "ksp_max_it": 1000})
return projection_problem
# Create a linear problem for projecting the stress tensor onto the linear tensor function space VV3.
#
tensor_projection_problem = setup_projection(Piola, VV3)
Piola_temp = tensor_projection_problem.solve()
# primary fields to write to output file
u_vis = Function(VV2, name="disp")
p_vis = Function(VV1, name="p")
phi_vis = Function(VV1, name="phi")
# Mises stress
T = Piola_temp*F.T/J
T0 = T - (1/3)*tr(T)*Identity(3)
Mises = sqrt((3/2)*inner(T0, T0))
Mises_vis= Function(VV1,name="Mises")
Mises_expr = Expression(Mises,VV1.element.interpolation_points())
# Cauchy stress components
T11 = Function(VV1)
T11.name = "T11"
T11_expr = Expression(T[0,0],VV1.element.interpolation_points())
T22 = Function(VV1)
T22.name = "T22"
T22_expr = Expression(T[1,1],VV1.element.interpolation_points())
T33 = Function(VV1)
T33.name = "T33"
T33_expr = Expression(T[2,2],VV1.element.interpolation_points())
# Stretch measure
I1_vis = Function(VV1)
I1_vis.name = "I1"
I1_expr = Expression(I1, VV1.element.interpolation_points())
# Volumetric deformation
J_vis = Function(VV1)
J_vis.name = "J"
J_expr = Expression(J, VV1.element.interpolation_points())
# set up the output VTX files.
file_results = VTXWriter(
MPI.COMM_WORLD,
"results/" + results_name + ".bp",
[ # put the functions here you wish to write to output
u_vis, p_vis, phi_vis, # DOF outputs
Mises_vis, T11, T22, T33, # stress outputs
I1_vis, J_vis, # Kinematical outputs
],
engine="BP4",
)
def writeResults(t):
# Update the output fields before writing to VTX.
#
u_vis.interpolate(w.sub(0))
p_vis.interpolate(w.sub(1))
phi_vis.interpolate(w.sub(2))
#
# re-project to smooth visualization of quadrature functions
# before interpolating.
Piola_temp = tensor_projection_problem.solve()
Mises_vis.interpolate(Mises_expr)
T11.interpolate(T11_expr)
T22.interpolate(T22_expr)
T33.interpolate(T33_expr)
#
I1_vis.interpolate(I1_expr)
J_vis.interpolate(J_expr)
# Finally, write output fields to VTX.
#
file_results.write(t)
Infrastructure for pulling out time history data (force, displacement, etc.)#
# # computing the reaction force using the stress field
# traction = dot(Piola_temp, n)
# Force = dot(traction, n)*ds(4)
# rxnForce = fem.form(Force)
# # infrastructure for evaluating functions at a certain point efficiently
pointForEval = np.array([length, length, 0.0])
bb_tree = dolfinx.geometry.bb_tree(domain,domain.topology.dim)
cell_candidates = dolfinx.geometry.compute_collisions_points(bb_tree, pointForEval)
colliding_cells = dolfinx.geometry.compute_colliding_cells(domain, cell_candidates, pointForEval).array
Boundary condtions#
# Constant for applied electric potential
phi_cons = Constant(domain,PETSc.ScalarType(phiRamp(0)))
# boundaries = [(1,Left),(2,Bottom),(3,Right),(4,Top)]
# Find the specific DOFs which will be constrained.
Left_u1_dofs = fem.locate_dofs_topological(ME.sub(0).sub(0), facet_tags.dim, facet_tags.find(1))
Bot_u2_dofs = fem.locate_dofs_topological(ME.sub(0).sub(1), facet_tags.dim, facet_tags.find(2))
Top_phi_dofs = fem.locate_dofs_topological(ME.sub(2), facet_tags.dim, facet_tags.find(4))
Bot_phi_dofs = fem.locate_dofs_topological(ME.sub(2), facet_tags.dim, facet_tags.find(2))
# building Dirichlet BCs
bcs_1 = dirichletbc(0.0, Left_u1_dofs, ME.sub(0).sub(0)) # u1 fix - Left
bcs_2 = dirichletbc(0.0, Bot_u2_dofs, ME.sub(0).sub(1)) # u2 fix - Bottom
#
bcs_3 = dirichletbc(phi_cons, Top_phi_dofs, ME.sub(2)) # phi ramp - Top
bcs_4 = dirichletbc(0.0, Bot_phi_dofs, ME.sub(2)) # phi ground - Bottom
bcs = [bcs_1, bcs_2, bcs_3, bcs_4]
Define the nonlinear variational problem#
# Set up nonlinear problem
problem = NonlinearProblem(Res, w, bcs, a)
# the global newton solver and params
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-8
solver.atol = 1e-8
solver.max_it = 50
solver.report = True
# The Krylov solver parameters.
ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly" # "preonly" works equally well
opts[f"{option_prefix}pc_type"] = "lu" # do not use 'gamg' pre-conditioner
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
opts[f"{option_prefix}ksp_max_it"] = 30
ksp.setFromOptions()
Start calculation loop#
# Give the step a descriptive name
step = "Actuate"
# Variables for storing time history
totSteps = 1000000
timeHist0 = np.zeros(shape=[totSteps])
timeHist1 = np.zeros(shape=[totSteps])
# Iinitialize a counter for reporting data
ii=0
# Write initial state to file
writeResults(t=0.0)
# print a message for simulation startup
print("------------------------------------")
print("Simulation Start")
print("------------------------------------")
# Store start time
startTime = datetime.now()
# Time-stepping solution procedure loop
while (round(t + dt, 9) <= Ttot):
# increment time
t += dt
# increment counter
ii += 1
# update time variables in time-dependent BCs
phi_cons.value = phiRamp(t)
# Solve the problem
try:
(iter, converged) = solver.solve(w)
except: # Break the loop if solver fails
break
# Collect results from MPI ghost processes
w.x.scatter_forward()
# Print progress of calculation
if ii%1 == 0:
now = datetime.now()
current_time = now.strftime("%H:%M:%S")
print("Step: {} | Increment: {} | Iterations: {}".format(step, ii, iter))
print("dt: {} | Simulation Time: {} s | Percent of total time: {}%".format(round(dt,4), round(t,4), round(100*t/Ttot,4)))
print()
# Write output to file
writeResults(t)
# Store time history variables at this time
timeHist0[ii] = w.sub(0).sub(1).eval([length, length, 0.0],colliding_cells[0])[0] # time history of displacement at a point
timeHist1[ii] = w.sub(2).eval([length, length, 0.0],colliding_cells[0])[0] # time history of voltage phi at a point
# Update DOFs for next step
w_old.x.array[:] = w.x.array
# close the output file.
file_results.close()
# End analysis
print("-----------------------------------------")
print("End computation")
# Report elapsed real time for the analysis
endTime = datetime.now()
elapseTime = endTime - startTime
print("------------------------------------------")
print("Elapsed real time: {}".format(elapseTime))
print("------------------------------------------")
------------------------------------
Simulation Start
------------------------------------
Step: Actuate | Increment: 1 | Iterations: 4
dt: 10.0 | Simulation Time: 10.0 s | Percent of total time: 1.0%
Step: Actuate | Increment: 2 | Iterations: 4
dt: 10.0 | Simulation Time: 20.0 s | Percent of total time: 2.0%
Step: Actuate | Increment: 3 | Iterations: 4
dt: 10.0 | Simulation Time: 30.0 s | Percent of total time: 3.0%
Step: Actuate | Increment: 4 | Iterations: 4
dt: 10.0 | Simulation Time: 40.0 s | Percent of total time: 4.0%
Step: Actuate | Increment: 5 | Iterations: 4
dt: 10.0 | Simulation Time: 50.0 s | Percent of total time: 5.0%
Step: Actuate | Increment: 6 | Iterations: 4
dt: 10.0 | Simulation Time: 60.0 s | Percent of total time: 6.0%
Step: Actuate | Increment: 7 | Iterations: 4
dt: 10.0 | Simulation Time: 70.0 s | Percent of total time: 7.0%
Step: Actuate | Increment: 8 | Iterations: 4
dt: 10.0 | Simulation Time: 80.0 s | Percent of total time: 8.0%
Step: Actuate | Increment: 9 | Iterations: 4
dt: 10.0 | Simulation Time: 90.0 s | Percent of total time: 9.0%
Step: Actuate | Increment: 10 | Iterations: 4
dt: 10.0 | Simulation Time: 100.0 s | Percent of total time: 10.0%
Step: Actuate | Increment: 11 | Iterations: 4
dt: 10.0 | Simulation Time: 110.0 s | Percent of total time: 11.0%
Step: Actuate | Increment: 12 | Iterations: 4
dt: 10.0 | Simulation Time: 120.0 s | Percent of total time: 12.0%
Step: Actuate | Increment: 13 | Iterations: 4
dt: 10.0 | Simulation Time: 130.0 s | Percent of total time: 13.0%
Step: Actuate | Increment: 14 | Iterations: 4
dt: 10.0 | Simulation Time: 140.0 s | Percent of total time: 14.0%
Step: Actuate | Increment: 15 | Iterations: 4
dt: 10.0 | Simulation Time: 150.0 s | Percent of total time: 15.0%
Step: Actuate | Increment: 16 | Iterations: 4
dt: 10.0 | Simulation Time: 160.0 s | Percent of total time: 16.0%
Step: Actuate | Increment: 17 | Iterations: 4
dt: 10.0 | Simulation Time: 170.0 s | Percent of total time: 17.0%
Step: Actuate | Increment: 18 | Iterations: 4
dt: 10.0 | Simulation Time: 180.0 s | Percent of total time: 18.0%
Step: Actuate | Increment: 19 | Iterations: 4
dt: 10.0 | Simulation Time: 190.0 s | Percent of total time: 19.0%
Step: Actuate | Increment: 20 | Iterations: 4
dt: 10.0 | Simulation Time: 200.0 s | Percent of total time: 20.0%
Step: Actuate | Increment: 21 | Iterations: 4
dt: 10.0 | Simulation Time: 210.0 s | Percent of total time: 21.0%
Step: Actuate | Increment: 22 | Iterations: 4
dt: 10.0 | Simulation Time: 220.0 s | Percent of total time: 22.0%
Step: Actuate | Increment: 23 | Iterations: 4
dt: 10.0 | Simulation Time: 230.0 s | Percent of total time: 23.0%
Step: Actuate | Increment: 24 | Iterations: 4
dt: 10.0 | Simulation Time: 240.0 s | Percent of total time: 24.0%
Step: Actuate | Increment: 25 | Iterations: 4
dt: 10.0 | Simulation Time: 250.0 s | Percent of total time: 25.0%
Step: Actuate | Increment: 26 | Iterations: 4
dt: 10.0 | Simulation Time: 260.0 s | Percent of total time: 26.0%
Step: Actuate | Increment: 27 | Iterations: 4
dt: 10.0 | Simulation Time: 270.0 s | Percent of total time: 27.0%
Step: Actuate | Increment: 28 | Iterations: 4
dt: 10.0 | Simulation Time: 280.0 s | Percent of total time: 28.0%
Step: Actuate | Increment: 29 | Iterations: 4
dt: 10.0 | Simulation Time: 290.0 s | Percent of total time: 29.0%
Step: Actuate | Increment: 30 | Iterations: 4
dt: 10.0 | Simulation Time: 300.0 s | Percent of total time: 30.0%
Step: Actuate | Increment: 31 | Iterations: 4
dt: 10.0 | Simulation Time: 310.0 s | Percent of total time: 31.0%
Step: Actuate | Increment: 32 | Iterations: 4
dt: 10.0 | Simulation Time: 320.0 s | Percent of total time: 32.0%
Step: Actuate | Increment: 33 | Iterations: 4
dt: 10.0 | Simulation Time: 330.0 s | Percent of total time: 33.0%
Step: Actuate | Increment: 34 | Iterations: 4
dt: 10.0 | Simulation Time: 340.0 s | Percent of total time: 34.0%
Step: Actuate | Increment: 35 | Iterations: 4
dt: 10.0 | Simulation Time: 350.0 s | Percent of total time: 35.0%
Step: Actuate | Increment: 36 | Iterations: 4
dt: 10.0 | Simulation Time: 360.0 s | Percent of total time: 36.0%
Step: Actuate | Increment: 37 | Iterations: 4
dt: 10.0 | Simulation Time: 370.0 s | Percent of total time: 37.0%
Step: Actuate | Increment: 38 | Iterations: 4
dt: 10.0 | Simulation Time: 380.0 s | Percent of total time: 38.0%
Step: Actuate | Increment: 39 | Iterations: 4
dt: 10.0 | Simulation Time: 390.0 s | Percent of total time: 39.0%
Step: Actuate | Increment: 40 | Iterations: 4
dt: 10.0 | Simulation Time: 400.0 s | Percent of total time: 40.0%
Step: Actuate | Increment: 41 | Iterations: 4
dt: 10.0 | Simulation Time: 410.0 s | Percent of total time: 41.0%
Step: Actuate | Increment: 42 | Iterations: 4
dt: 10.0 | Simulation Time: 420.0 s | Percent of total time: 42.0%
Step: Actuate | Increment: 43 | Iterations: 4
dt: 10.0 | Simulation Time: 430.0 s | Percent of total time: 43.0%
Step: Actuate | Increment: 44 | Iterations: 4
dt: 10.0 | Simulation Time: 440.0 s | Percent of total time: 44.0%
Step: Actuate | Increment: 45 | Iterations: 4
dt: 10.0 | Simulation Time: 450.0 s | Percent of total time: 45.0%
Step: Actuate | Increment: 46 | Iterations: 4
dt: 10.0 | Simulation Time: 460.0 s | Percent of total time: 46.0%
Step: Actuate | Increment: 47 | Iterations: 4
dt: 10.0 | Simulation Time: 470.0 s | Percent of total time: 47.0%
Step: Actuate | Increment: 48 | Iterations: 4
dt: 10.0 | Simulation Time: 480.0 s | Percent of total time: 48.0%
Step: Actuate | Increment: 49 | Iterations: 4
dt: 10.0 | Simulation Time: 490.0 s | Percent of total time: 49.0%
Step: Actuate | Increment: 50 | Iterations: 4
dt: 10.0 | Simulation Time: 500.0 s | Percent of total time: 50.0%
Step: Actuate | Increment: 51 | Iterations: 4
dt: 10.0 | Simulation Time: 510.0 s | Percent of total time: 51.0%
Step: Actuate | Increment: 52 | Iterations: 4
dt: 10.0 | Simulation Time: 520.0 s | Percent of total time: 52.0%
Step: Actuate | Increment: 53 | Iterations: 4
dt: 10.0 | Simulation Time: 530.0 s | Percent of total time: 53.0%
Step: Actuate | Increment: 54 | Iterations: 4
dt: 10.0 | Simulation Time: 540.0 s | Percent of total time: 54.0%
Step: Actuate | Increment: 55 | Iterations: 4
dt: 10.0 | Simulation Time: 550.0 s | Percent of total time: 55.0%
Step: Actuate | Increment: 56 | Iterations: 4
dt: 10.0 | Simulation Time: 560.0 s | Percent of total time: 56.0%
Step: Actuate | Increment: 57 | Iterations: 4
dt: 10.0 | Simulation Time: 570.0 s | Percent of total time: 57.0%
Step: Actuate | Increment: 58 | Iterations: 4
dt: 10.0 | Simulation Time: 580.0 s | Percent of total time: 58.0%
Step: Actuate | Increment: 59 | Iterations: 5
dt: 10.0 | Simulation Time: 590.0 s | Percent of total time: 59.0%
Step: Actuate | Increment: 60 | Iterations: 5
dt: 10.0 | Simulation Time: 600.0 s | Percent of total time: 60.0%
Step: Actuate | Increment: 61 | Iterations: 5
dt: 10.0 | Simulation Time: 610.0 s | Percent of total time: 61.0%
Step: Actuate | Increment: 62 | Iterations: 5
dt: 10.0 | Simulation Time: 620.0 s | Percent of total time: 62.0%
Step: Actuate | Increment: 63 | Iterations: 6
dt: 10.0 | Simulation Time: 630.0 s | Percent of total time: 63.0%
Step: Actuate | Increment: 64 | Iterations: 5
dt: 10.0 | Simulation Time: 640.0 s | Percent of total time: 64.0%
Step: Actuate | Increment: 65 | Iterations: 6
dt: 10.0 | Simulation Time: 650.0 s | Percent of total time: 65.0%
Step: Actuate | Increment: 66 | Iterations: 5
dt: 10.0 | Simulation Time: 660.0 s | Percent of total time: 66.0%
Step: Actuate | Increment: 67 | Iterations: 5
dt: 10.0 | Simulation Time: 670.0 s | Percent of total time: 67.0%
Step: Actuate | Increment: 68 | Iterations: 5
dt: 10.0 | Simulation Time: 680.0 s | Percent of total time: 68.0%
Step: Actuate | Increment: 69 | Iterations: 5
dt: 10.0 | Simulation Time: 690.0 s | Percent of total time: 69.0%
Step: Actuate | Increment: 70 | Iterations: 5
dt: 10.0 | Simulation Time: 700.0 s | Percent of total time: 70.0%
Step: Actuate | Increment: 71 | Iterations: 4
dt: 10.0 | Simulation Time: 710.0 s | Percent of total time: 71.0%
Step: Actuate | Increment: 72 | Iterations: 4
dt: 10.0 | Simulation Time: 720.0 s | Percent of total time: 72.0%
Step: Actuate | Increment: 73 | Iterations: 4
dt: 10.0 | Simulation Time: 730.0 s | Percent of total time: 73.0%
Step: Actuate | Increment: 74 | Iterations: 4
dt: 10.0 | Simulation Time: 740.0 s | Percent of total time: 74.0%
Step: Actuate | Increment: 75 | Iterations: 4
dt: 10.0 | Simulation Time: 750.0 s | Percent of total time: 75.0%
Step: Actuate | Increment: 76 | Iterations: 4
dt: 10.0 | Simulation Time: 760.0 s | Percent of total time: 76.0%
Step: Actuate | Increment: 77 | Iterations: 4
dt: 10.0 | Simulation Time: 770.0 s | Percent of total time: 77.0%
Step: Actuate | Increment: 78 | Iterations: 4
dt: 10.0 | Simulation Time: 780.0 s | Percent of total time: 78.0%
Step: Actuate | Increment: 79 | Iterations: 4
dt: 10.0 | Simulation Time: 790.0 s | Percent of total time: 79.0%
Step: Actuate | Increment: 80 | Iterations: 4
dt: 10.0 | Simulation Time: 800.0 s | Percent of total time: 80.0%
Step: Actuate | Increment: 81 | Iterations: 4
dt: 10.0 | Simulation Time: 810.0 s | Percent of total time: 81.0%
Step: Actuate | Increment: 82 | Iterations: 4
dt: 10.0 | Simulation Time: 820.0 s | Percent of total time: 82.0%
Step: Actuate | Increment: 83 | Iterations: 4
dt: 10.0 | Simulation Time: 830.0 s | Percent of total time: 83.0%
Step: Actuate | Increment: 84 | Iterations: 4
dt: 10.0 | Simulation Time: 840.0 s | Percent of total time: 84.0%
Step: Actuate | Increment: 85 | Iterations: 4
dt: 10.0 | Simulation Time: 850.0 s | Percent of total time: 85.0%
Step: Actuate | Increment: 86 | Iterations: 4
dt: 10.0 | Simulation Time: 860.0 s | Percent of total time: 86.0%
Step: Actuate | Increment: 87 | Iterations: 4
dt: 10.0 | Simulation Time: 870.0 s | Percent of total time: 87.0%
Step: Actuate | Increment: 88 | Iterations: 4
dt: 10.0 | Simulation Time: 880.0 s | Percent of total time: 88.0%
Step: Actuate | Increment: 89 | Iterations: 4
dt: 10.0 | Simulation Time: 890.0 s | Percent of total time: 89.0%
Step: Actuate | Increment: 90 | Iterations: 4
dt: 10.0 | Simulation Time: 900.0 s | Percent of total time: 90.0%
Step: Actuate | Increment: 91 | Iterations: 4
dt: 10.0 | Simulation Time: 910.0 s | Percent of total time: 91.0%
Step: Actuate | Increment: 92 | Iterations: 4
dt: 10.0 | Simulation Time: 920.0 s | Percent of total time: 92.0%
Step: Actuate | Increment: 93 | Iterations: 4
dt: 10.0 | Simulation Time: 930.0 s | Percent of total time: 93.0%
Step: Actuate | Increment: 94 | Iterations: 4
dt: 10.0 | Simulation Time: 940.0 s | Percent of total time: 94.0%
Step: Actuate | Increment: 95 | Iterations: 4
dt: 10.0 | Simulation Time: 950.0 s | Percent of total time: 95.0%
Step: Actuate | Increment: 96 | Iterations: 4
dt: 10.0 | Simulation Time: 960.0 s | Percent of total time: 96.0%
Step: Actuate | Increment: 97 | Iterations: 4
dt: 10.0 | Simulation Time: 970.0 s | Percent of total time: 97.0%
Step: Actuate | Increment: 98 | Iterations: 4
dt: 10.0 | Simulation Time: 980.0 s | Percent of total time: 98.0%
Step: Actuate | Increment: 99 | Iterations: 4
dt: 10.0 | Simulation Time: 990.0 s | Percent of total time: 99.0%
Step: Actuate | Increment: 100 | Iterations: 4
dt: 10.0 | Simulation Time: 1000.0 s | Percent of total time: 100.0%
-----------------------------------------
End computation
------------------------------------------
Elapsed real time: 0:00:02.400832
------------------------------------------
Plot results#
# set plot font to size 18
font = {'size' : 18}
plt.rc('font', **font)
# Get array of default plot colors
prop_cycle = plt.rcParams['axes.prop_cycle']
colors = prop_cycle.by_key()['color']
# Plot the normalized dimensionless quantity for $\phi$ used in Wang et al. 2016
# versus stretch in the vertical direction.
#
normVolts = timeHist1/(length * np.sqrt(float(Geq_0)/float(vareps)))
normVolts = normVolts[0:ii]
#
stretch = timeHist0/length + 1.0
stretch = stretch[0:ii]
#
plt.plot(normVolts, stretch, c=colors[0], linewidth=2.0)
plt.grid(linestyle="--", linewidth=0.5, color='b')
ax = plt.gca()
#
ax.set_ylabel(r'$\lambda$')
# ax.set_ylim([0.6,1.05])
# ax.set_yticks([0.2, 0.4, 0.6, 0.8, 1.0])
#
# ax.set_xlabel(r'$(\phi/\ell_0)/\sqrt{G_0/\varepsilon} $')
ax.set_xlabel(r'Normalized potential, $\frac{\phi/\ell_0}{\sqrt{G_0/\varepsilon}}$',size=18)
# ax.set_xlim([0.0,0.8])
# ax.set_xticks([0.0, 0.25, 0.5, 0.75, 1.0])
#
# from matplotlib.ticker import AutoMinorLocator,FormatStrFormatter
# ax.xaxis.set_minor_locator(AutoMinorLocator())
# ax.yaxis.set_minor_locator(AutoMinorLocator())
# plt.show()
# plt.text(normVolts[ii-1], stretch[ii-1], r'X', color='k',\
# horizontalalignment='center',verticalalignment='center')
fig = plt.gcf()
fig.set_size_inches(7,5)
plt.tight_layout()
plt.savefig("results/pe_pullin_Im175.png", dpi=600)
