Buckling of a column with square cross-section#
Units#
Length: mm
Mass: kg
Time: s
Force: milliNewtons
Stress: kPa
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)
# basix finite elements (necessary for dolfinx v0.8.0)
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#
# A 3-D beam
width = 1.0 # mm
length = 20 # mm
domain = mesh.create_box(MPI.COMM_WORLD, [[0.0,0.0,0.0], [width, width,length]],\
[4,4,50], mesh.CellType.tetrahedron)
# Add an initial imperfection to control buckling mode
imperf = 0.005*width # mm
# Map the coordinates of the uniform box mesh to the biased spacing
xOrig = domain.geometry.x
xMap1 = np.zeros((len(xOrig),3))
# Mapping functions
for i in range(0,len(xMap1)):
xMap1[i,0] = xOrig[i,0] - (imperf/2.0)*(1.0 - np.cos(2*np.pi*xOrig[i,2]/length))
xMap1[i,1] = xOrig[i,1] - (imperf/2.0)*(1.0 - np.cos(2*np.pi*xOrig[i,2]/length))
xMap1[i,2] = xOrig[i,2]
domain.geometry.x[:] = xMap1
# This says "spatial coordinates" but is really the referential coordinates,
# since the mesh does not convect in FEniCS.
x = ufl.SpatialCoordinate(domain)
Identify boundaries of the domain
# Identify the planar boundaries of the box mesh
#
def xBot(x):
return np.isclose(x[0], 0)
def xTop(x):
return np.isclose(x[0], width)
def yBot(x):
return np.isclose(x[1], 0)
def yTop(x):
return np.isclose(x[1], width)
def zBot(x):
return np.isclose(x[2], 0)
def zTop(x):
return np.isclose(x[2], length)
# Mark the sub-domains
boundaries = [(1, xBot),(2,xTop),(3,yBot),(4,yTop),(5,zBot),(6,zTop)]
# 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
# Surface numbering:
# boundaries = [(1, xBot),(2,xTop),(3,yBot),(4,yTop),(5,zBot),(6,zTop)]
[5 6]
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 -- I make the 3D mesh opaque, so that 2D surfaces stand out.
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, opacity=0.5)
# Add colored 2D surfaces for the named surfaces
zBot_surf = pyvista.UnstructuredGrid(*vtk_mesh(domain, domain.topology.dim-1,facet_tags.indices[facet_tags.values==5]) )
zTop_surf = pyvista.UnstructuredGrid(*vtk_mesh(domain, domain.topology.dim-1,facet_tags.indices[facet_tags.values==6]) )
#
actor = plotter.add_mesh(zBot_surf, show_edges=True,color="blue") # zBot face is blue
actor2 = plotter.add_mesh(zTop_surf, show_edges=True,color="red") # zTop is red
labels = dict(zlabel='Z', xlabel='X', ylabel='Y')
plotter.add_axes(**labels)
plotter.screenshot("mesh.png")
from IPython.display import Image
Image(filename='mesh.png')
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':4})
# Define the volume integration measure "dx"
# also specify the number of volume quadrature points.
dx = ufl.Measure('dx', domain=domain, metadata={'quadrature_degree': 4})
# Define facet normal
n = ufl.FacetNormal(domain)
Material parameters#
-Arruda-Boyce model
Gshear_0 = Constant(domain,PETSc.ScalarType(280.0)) # Ground state shear modulus
lambdaL = Constant(domain,PETSc.ScalarType(5.12)) # Locking stretch
Kbulk = Constant(domain,PETSc.ScalarType(1000.0*Gshear_0))
Function spaces#
# Define function space, both vectorial and scalar
# dolfinx v0.8.0 syntax:
U2 = element("Lagrange", domain.basix_cell(), 2, shape=(3,)) # For displacement
P1 = element("Lagrange", domain.basix_cell(), 1) # For pressure
#
TH = mixed_element([U2, P1]) # Taylor-Hood style mixed element
ME = functionspace(domain, TH) # Total space for all DOFs
# Define actual functions with the required DOFs
w = Function(ME)
u, p = split(w) # displacement u, pressure p
# A copy of functions to store values in the previous step
w_old = Function(ME)
u_old, p_old = split(w_old)
# Define test functions
u_test, p_test = TestFunctions(ME)
# Define trial functions needed for automatic differentiation
dw = TrialFunction(ME)
Initial conditions#
The initial conditions for degrees of freedom u and p are zero everywhere
These are imposed automatically, since we have not specified any non-zero initial conditions.
Subroutines for kinematics and constitutive equations#
# Deformation gradient
def F_calc(u):
Id = Identity(3)
F = Id + grad(u)
return F
def lambdaBar_calc(u):
F = F_calc(u)
C = F.T*F
Cdis = J**(-2/3)*C
I1 = tr(Cdis)
lambdaBar = sqrt(I1/3.0)
return lambdaBar
def zeta_calc(u):
lambdaBar = lambdaBar_calc(u)
# Use Pade approximation of Langevin inverse
z = lambdaBar/lambdaL
z = conditional(gt(z,0.95), 0.95, z) # Keep simulation from blowing up
beta = z*(3.0 - z**2.0)/(1.0 - z**2.0)
zeta = (lambdaL/(3*lambdaBar))*beta
return zeta
# Generalized shear modulus for Arruda-Boyce model
def Gshear_AB_calc(u):
zeta = zeta_calc(u)
Gshear = Gshear_0 * zeta
return Gshear
#---------------------------------------------
# Subroutine for calculating the Cauchy stress
#---------------------------------------------
def T_calc(u,p):
Id = Identity(3)
F = F_calc(u)
J = det(F)
B = F*F.T
Bdis = J**(-2/3)*B
Gshear = Gshear_AB_calc(u)
T = (1/J)* Gshear * dev(Bdis) - p * Id
return T
#----------------------------------------------
# Subroutine for calculating the Piola stress
#----------------------------------------------
def Piola_calc(u, p):
Id = Identity(3)
F = F_calc(u)
J = det(F)
#
T = T_calc(u,p)
#
Tmat = J * T * inv(F.T)
return Tmat
Evaluate kinematics and constitutive relations#
F = F_calc(u)
J = det(F)
lambdaBar = lambdaBar_calc(u)
# Piola stress
Tmat = Piola_calc(u, p)
Weak forms#
# Residuals:
# Res_0: Balance of forces (test fxn: u)
# Res_1: Coupling pressure (test fxn: p)
# The weak form for the balance of forces
Res_0 = inner(Tmat, grad(u_test) )*dx
# The weak form for the pressure
fac_p = ln(J)/J
#
Res_1 = dot( (p/Kbulk + fac_p), p_test)*dx
# Total weak form
Res = Res_0 + Res_1
# Automatic differentiation tangent:
a = derivative(Res, w, dw)
Set-up output files#
# results file name
results_name = "3D_column_buckling"
# Function space for projection of results
# v0.8.0 syntax:
U1 = element("DG", domain.basix_cell(), 1, shape=(3,)) # For displacement
P0 = element("DG", domain.basix_cell(), 1) # For pressure
V2 = fem.functionspace(domain, U1) #Vector function space
V1 = fem.functionspace(domain, P0) #Scalar function space
# fields to write to output file
u_vis = Function(V2)
u_vis.name = "disp"
p_vis = Function(V1)
p_vis.name = "p"
J_vis = Function(V1)
J_vis.name = "J"
J_expr = Expression(J,V1.element.interpolation_points())
lambdaBar_vis = Function(V1)
lambdaBar_vis.name = "lambdaBar"
lambdaBar_expr = Expression(lambdaBar,V1.element.interpolation_points())
P11 = Function(V1)
P11.name = "P11"
P11_expr = Expression(Tmat[0,0],V1.element.interpolation_points())
P22 = Function(V1)
P22.name = "P22"
P22_expr = Expression(Tmat[1,1],V1.element.interpolation_points())
P33 = Function(V1)
P33.name = "P33"
P33_expr = Expression(Tmat[2,2],V1.element.interpolation_points())
T = Tmat*F.T/J
T0 = T - (1/3)*tr(T)*Identity(3)
Mises = sqrt((3/2)*inner(T0, T0))
Mises_vis= Function(V1,name="Mises")
Mises_expr = Expression(Mises,V1.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, J_vis, P11, P22, P33, lambdaBar_vis,
Mises_vis,
],
engine="BP4",
)
def writeResults(t):
# Output field interpolation
u_vis.interpolate(w.sub(0))
p_vis.interpolate(w.sub(1))
J_vis.interpolate(J_expr)
P11.interpolate(P11_expr)
P22.interpolate(P22_expr)
P33.interpolate(P33_expr)
lambdaBar_vis.interpolate(lambdaBar_expr)
Mises_vis.interpolate(Mises_expr)
# Write output fields
file_results.write(t)
Infrastructure for pulling out time history data (force, displacement, etc.)#
# infrastructure for evaluating functions at a certain point efficiently
pointForDisp = np.array([width, width, length])
bb_tree = dolfinx.geometry.bb_tree(domain,domain.topology.dim)
cell_candidates = dolfinx.geometry.compute_collisions_points(bb_tree, pointForDisp)
colliding_cells = dolfinx.geometry.compute_colliding_cells(domain, cell_candidates, pointForDisp).array
# Form for evaluating the engineering stress at each step
engineeringStress = fem.form(P33*ds(6))
# Surface numbering:
# boundaries = [(1, xBot),(2,xTop),(3,yBot),(4,yTop),(5,zBot),(6,zTop)]
Name the analysis step#
# Give the step a descriptive name
step = "Buckle"
Boundary condtions#
# Surface numbering:
# boundaries = [(1, xBot),(2,xTop),(3,yBot),(4,yTop),(5,zBot),(6,zTop)]
# Constant for applied displacement
dispCons = Constant(domain,PETSc.ScalarType(dispRamp(0)))
# Find the specific DOFs which will be constrained.
#
# Bottom surface displacement degrees of freedom
Btm_dofs_u1 = fem.locate_dofs_topological(ME.sub(0).sub(0), facet_tags.dim, facet_tags.find(5))
Btm_dofs_u2 = fem.locate_dofs_topological(ME.sub(0).sub(1), facet_tags.dim, facet_tags.find(5))
Btm_dofs_u3 = fem.locate_dofs_topological(ME.sub(0).sub(2), facet_tags.dim, facet_tags.find(5))
# Top surface displacement degrees of freedom
Top_dofs_u1 = fem.locate_dofs_topological(ME.sub(0).sub(0), facet_tags.dim, facet_tags.find(6))
Top_dofs_u2 = fem.locate_dofs_topological(ME.sub(0).sub(1), facet_tags.dim, facet_tags.find(6))
Top_dofs_u3 = fem.locate_dofs_topological(ME.sub(0).sub(2), facet_tags.dim, facet_tags.find(6))
# Build the Dirichlet BCs
bcs_0 = dirichletbc(0.0, Btm_dofs_u1, ME.sub(0).sub(0)) # u1 fix - zBtm
bcs_1 = dirichletbc(0.0, Btm_dofs_u2, ME.sub(0).sub(1)) # u2 fix - zBtm
bcs_2 = dirichletbc(0.0, Btm_dofs_u3, ME.sub(0).sub(2)) # u3 fix - zBtm
#
bcs_3 = dirichletbc(0.0, Top_dofs_u1, ME.sub(0).sub(0)) # u1 fix - zTop
bcs_4 = dirichletbc(0.0, Top_dofs_u2, ME.sub(0).sub(1)) # u2 fix - zTop
bcs_5 = dirichletbc(dispCons, Top_dofs_u3, ME.sub(0).sub(2)) # u3 move - zTop
# collect all BCs in one object.
bcs = [bcs_0, bcs_1, bcs_2, bcs_3, bcs_4, bcs_5]
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#
# Variables for storing time history
totSteps = numSteps+1
timeHist0 = np.zeros(shape=[totSteps])
timeHist1 = np.zeros(shape=[totSteps])
timeHist2 = np.zeros(shape=[totSteps])
#Iinitialize a counter for reporting data
ii=0
# Write initial state to file
writeResults(t=0.0)
# Print out message for simulation start
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
dispCons.value = dispRamp(t)
# Solve the problem
try:
(iter, converged) = solver.solve(w)
except: # Break the loop if solver fails
print("Ended Early")
break
# Collect results from MPI ghost processes
w.x.scatter_forward()
# Write output to file
writeResults(t)
# Update DOFs for next step
w_old.x.array[:] = w.x.array
# Store time history variables at this time
timeHist0[ii] = t # current time
#
timeHist1[ii] = dispRamp(t) # time history of applied pressure
#
timeHist2[ii] = domain.comm.gather(fem.assemble_scalar(engineeringStress))[0] # time history of engineering stress
# 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(" Simulation Time: {} s of {} s".\
format(round(t,4), Ttot))
print()
# 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: Buckle | Increment: 1, Iterations: 4
Simulation Time: 0.01 s of 1 s
Step: Buckle | Increment: 2, Iterations: 4
Simulation Time: 0.02 s of 1 s
Step: Buckle | Increment: 3, Iterations: 4
Simulation Time: 0.03 s of 1 s
Step: Buckle | Increment: 4, Iterations: 4
Simulation Time: 0.04 s of 1 s
Step: Buckle | Increment: 5, Iterations: 4
Simulation Time: 0.05 s of 1 s
Step: Buckle | Increment: 6, Iterations: 5
Simulation Time: 0.06 s of 1 s
Step: Buckle | Increment: 7, Iterations: 11
Simulation Time: 0.07 s of 1 s
Step: Buckle | Increment: 8, Iterations: 6
Simulation Time: 0.08 s of 1 s
Step: Buckle | Increment: 9, Iterations: 5
Simulation Time: 0.09 s of 1 s
Step: Buckle | Increment: 10, Iterations: 5
Simulation Time: 0.1 s of 1 s
Step: Buckle | Increment: 11, Iterations: 5
Simulation Time: 0.11 s of 1 s
Step: Buckle | Increment: 12, Iterations: 5
Simulation Time: 0.12 s of 1 s
Step: Buckle | Increment: 13, Iterations: 4
Simulation Time: 0.13 s of 1 s
Step: Buckle | Increment: 14, Iterations: 4
Simulation Time: 0.14 s of 1 s
Step: Buckle | Increment: 15, Iterations: 4
Simulation Time: 0.15 s of 1 s
Step: Buckle | Increment: 16, Iterations: 4
Simulation Time: 0.16 s of 1 s
Step: Buckle | Increment: 17, Iterations: 4
Simulation Time: 0.17 s of 1 s
Step: Buckle | Increment: 18, Iterations: 4
Simulation Time: 0.18 s of 1 s
Step: Buckle | Increment: 19, Iterations: 4
Simulation Time: 0.19 s of 1 s
Step: Buckle | Increment: 20, Iterations: 4
Simulation Time: 0.2 s of 1 s
Step: Buckle | Increment: 21, Iterations: 4
Simulation Time: 0.21 s of 1 s
Step: Buckle | Increment: 22, Iterations: 4
Simulation Time: 0.22 s of 1 s
Step: Buckle | Increment: 23, Iterations: 4
Simulation Time: 0.23 s of 1 s
Step: Buckle | Increment: 24, Iterations: 4
Simulation Time: 0.24 s of 1 s
Step: Buckle | Increment: 25, Iterations: 4
Simulation Time: 0.25 s of 1 s
Step: Buckle | Increment: 26, Iterations: 4
Simulation Time: 0.26 s of 1 s
Step: Buckle | Increment: 27, Iterations: 4
Simulation Time: 0.27 s of 1 s
Step: Buckle | Increment: 28, Iterations: 4
Simulation Time: 0.28 s of 1 s
Step: Buckle | Increment: 29, Iterations: 4
Simulation Time: 0.29 s of 1 s
Step: Buckle | Increment: 30, Iterations: 4
Simulation Time: 0.3 s of 1 s
Step: Buckle | Increment: 31, Iterations: 4
Simulation Time: 0.31 s of 1 s
Step: Buckle | Increment: 32, Iterations: 4
Simulation Time: 0.32 s of 1 s
Step: Buckle | Increment: 33, Iterations: 4
Simulation Time: 0.33 s of 1 s
Step: Buckle | Increment: 34, Iterations: 4
Simulation Time: 0.34 s of 1 s
Step: Buckle | Increment: 35, Iterations: 4
Simulation Time: 0.35 s of 1 s
Step: Buckle | Increment: 36, Iterations: 4
Simulation Time: 0.36 s of 1 s
Step: Buckle | Increment: 37, Iterations: 4
Simulation Time: 0.37 s of 1 s
Step: Buckle | Increment: 38, Iterations: 4
Simulation Time: 0.38 s of 1 s
Step: Buckle | Increment: 39, Iterations: 4
Simulation Time: 0.39 s of 1 s
Step: Buckle | Increment: 40, Iterations: 4
Simulation Time: 0.4 s of 1 s
Step: Buckle | Increment: 41, Iterations: 4
Simulation Time: 0.41 s of 1 s
Step: Buckle | Increment: 42, Iterations: 4
Simulation Time: 0.42 s of 1 s
Step: Buckle | Increment: 43, Iterations: 4
Simulation Time: 0.43 s of 1 s
Step: Buckle | Increment: 44, Iterations: 4
Simulation Time: 0.44 s of 1 s
Step: Buckle | Increment: 45, Iterations: 4
Simulation Time: 0.45 s of 1 s
Step: Buckle | Increment: 46, Iterations: 4
Simulation Time: 0.46 s of 1 s
Step: Buckle | Increment: 47, Iterations: 4
Simulation Time: 0.47 s of 1 s
Step: Buckle | Increment: 48, Iterations: 4
Simulation Time: 0.48 s of 1 s
Step: Buckle | Increment: 49, Iterations: 4
Simulation Time: 0.49 s of 1 s
Step: Buckle | Increment: 50, Iterations: 4
Simulation Time: 0.5 s of 1 s
Step: Buckle | Increment: 51, Iterations: 4
Simulation Time: 0.51 s of 1 s
Step: Buckle | Increment: 52, Iterations: 4
Simulation Time: 0.52 s of 1 s
Step: Buckle | Increment: 53, Iterations: 4
Simulation Time: 0.53 s of 1 s
Step: Buckle | Increment: 54, Iterations: 4
Simulation Time: 0.54 s of 1 s
Step: Buckle | Increment: 55, Iterations: 4
Simulation Time: 0.55 s of 1 s
Step: Buckle | Increment: 56, Iterations: 4
Simulation Time: 0.56 s of 1 s
Step: Buckle | Increment: 57, Iterations: 4
Simulation Time: 0.57 s of 1 s
Step: Buckle | Increment: 58, Iterations: 4
Simulation Time: 0.58 s of 1 s
Step: Buckle | Increment: 59, Iterations: 4
Simulation Time: 0.59 s of 1 s
Step: Buckle | Increment: 60, Iterations: 4
Simulation Time: 0.6 s of 1 s
Step: Buckle | Increment: 61, Iterations: 4
Simulation Time: 0.61 s of 1 s
Step: Buckle | Increment: 62, Iterations: 4
Simulation Time: 0.62 s of 1 s
Step: Buckle | Increment: 63, Iterations: 4
Simulation Time: 0.63 s of 1 s
Step: Buckle | Increment: 64, Iterations: 4
Simulation Time: 0.64 s of 1 s
Step: Buckle | Increment: 65, Iterations: 4
Simulation Time: 0.65 s of 1 s
Step: Buckle | Increment: 66, Iterations: 4
Simulation Time: 0.66 s of 1 s
Step: Buckle | Increment: 67, Iterations: 4
Simulation Time: 0.67 s of 1 s
Step: Buckle | Increment: 68, Iterations: 4
Simulation Time: 0.68 s of 1 s
Step: Buckle | Increment: 69, Iterations: 4
Simulation Time: 0.69 s of 1 s
Step: Buckle | Increment: 70, Iterations: 4
Simulation Time: 0.7 s of 1 s
Step: Buckle | Increment: 71, Iterations: 4
Simulation Time: 0.71 s of 1 s
Step: Buckle | Increment: 72, Iterations: 4
Simulation Time: 0.72 s of 1 s
Step: Buckle | Increment: 73, Iterations: 4
Simulation Time: 0.73 s of 1 s
Step: Buckle | Increment: 74, Iterations: 4
Simulation Time: 0.74 s of 1 s
Step: Buckle | Increment: 75, Iterations: 4
Simulation Time: 0.75 s of 1 s
Step: Buckle | Increment: 76, Iterations: 4
Simulation Time: 0.76 s of 1 s
Step: Buckle | Increment: 77, Iterations: 4
Simulation Time: 0.77 s of 1 s
Step: Buckle | Increment: 78, Iterations: 4
Simulation Time: 0.78 s of 1 s
Step: Buckle | Increment: 79, Iterations: 4
Simulation Time: 0.79 s of 1 s
Step: Buckle | Increment: 80, Iterations: 4
Simulation Time: 0.8 s of 1 s
Step: Buckle | Increment: 81, Iterations: 4
Simulation Time: 0.81 s of 1 s
Step: Buckle | Increment: 82, Iterations: 4
Simulation Time: 0.82 s of 1 s
Step: Buckle | Increment: 83, Iterations: 4
Simulation Time: 0.83 s of 1 s
Step: Buckle | Increment: 84, Iterations: 4
Simulation Time: 0.84 s of 1 s
Step: Buckle | Increment: 85, Iterations: 4
Simulation Time: 0.85 s of 1 s
Step: Buckle | Increment: 86, Iterations: 4
Simulation Time: 0.86 s of 1 s
Step: Buckle | Increment: 87, Iterations: 4
Simulation Time: 0.87 s of 1 s
Step: Buckle | Increment: 88, Iterations: 4
Simulation Time: 0.88 s of 1 s
Step: Buckle | Increment: 89, Iterations: 4
Simulation Time: 0.89 s of 1 s
Step: Buckle | Increment: 90, Iterations: 4
Simulation Time: 0.9 s of 1 s
Step: Buckle | Increment: 91, Iterations: 4
Simulation Time: 0.91 s of 1 s
Step: Buckle | Increment: 92, Iterations: 4
Simulation Time: 0.92 s of 1 s
Step: Buckle | Increment: 93, Iterations: 4
Simulation Time: 0.93 s of 1 s
Step: Buckle | Increment: 94, Iterations: 4
Simulation Time: 0.94 s of 1 s
Step: Buckle | Increment: 95, Iterations: 4
Simulation Time: 0.95 s of 1 s
Step: Buckle | Increment: 96, Iterations: 4
Simulation Time: 0.96 s of 1 s
Step: Buckle | Increment: 97, Iterations: 4
Simulation Time: 0.97 s of 1 s
Step: Buckle | Increment: 98, Iterations: 4
Simulation Time: 0.98 s of 1 s
Step: Buckle | Increment: 99, Iterations: 4
Simulation Time: 0.99 s of 1 s
Step: Buckle | Increment: 100, Iterations: 4
Simulation Time: 1.0 s of 1 s
-----------------------------------------
End computation
------------------------------------------
Elapsed real time: 0:03:25.355426
------------------------------------------
Plot results#
# set plot font to size 14
font = {'size' : 14}
plt.rc('font', **font)
# Get array of default plot colors
prop_cycle = plt.rcParams['axes.prop_cycle']
colors = prop_cycle.by_key()['color']
# Only plot as far as we have time history data
ind = np.argmax(timeHist0)
#
fig = plt.figure()
#fig.set_size_inches(7,4)
ax=fig.gca()
#--------------------------------------------------------------
plt.plot(-timeHist1[0:ind], -timeHist2[0:ind], linewidth=2.0,\
color=colors[0], marker='.')
#-------------------------------------------------------------
# plt.xlim([1,8])
# plt.ylim([0,8])
plt.axis('tight')
plt.ylabel(r'Force, mN')
plt.xlabel(r'$u_3$ (mm)')
plt.grid(linestyle="--", linewidth=0.5, color='b')
#-------------------------------------------------------------
#ax.set_title("Stress dispalcement curve", size=14, weight='normal')
from matplotlib.ticker import AutoMinorLocator,FormatStrFormatter
ax.xaxis.set_minor_locator(AutoMinorLocator())
ax.yaxis.set_minor_locator(AutoMinorLocator())
plt.show()
fig = plt.gcf()
fig.set_size_inches(7,5)
plt.tight_layout()
plt.savefig("results/3D_column_buckling.png", dpi=600)
<Figure size 700x500 with 0 Axes>