“Ideal” sphere magnetostriction#
Code for large deformation soft-magentoelasticity of elastomeric materials
Magnetic actuation of a spherical specimen using an Axisymmetric formulation.
Units#
Basic:
Length: mm
Time: s
Mass: kg
Charge: kC
Derived:
Force: mN
Pressure: kPa
Current: kA
Mag. flux density: mT
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, sin, acos, ge, le, outer, tanh, cosh, atan, atan2)
# 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#
# Geometry parameters
radius = 5.0 # radius of sphere, mm
scale = 50 # side length of air domain, mm
# Read in the 2D mesh and cell tags
with XDMFFile(MPI.COMM_WORLD,"meshes/sphere_air_5mm.xdmf",'r') as infile:
domain = infile.read_mesh(name="Grid",xpath="/Xdmf/Domain")
cell_tags = infile.read_meshtags(domain,name="Grid")
domain.topology.create_connectivity(domain.topology.dim, domain.topology.dim-1)
# Also read in 1D facets for applying BCs
with XDMFFile(MPI.COMM_WORLD,"meshes/facet_sphere_air_5mm.xdmf",'r') as infile:
facet_tags = infile.read_meshtags(domain,name="Grid")
# A single point for "grounding" the displacement and phi
def ground(x):
# return np.logical_and(np.isclose(x[0], 0), np.isclose(x[1], 0))
return np.logical_and(np.isclose(x[1], 0), np.less_equal(x[0], radius))
x = ufl.SpatialCoordinate(domain)
Print out the unique cell index numbers
top_imap = domain.topology.index_map(2) # index map of 2D entities in domain
values = np.zeros(top_imap.size_global) # an array of zeros of the same size as number of 2D entities
values[cell_tags.indices]=cell_tags.values # populating the array with facet tag index numbers
print(np.unique(cell_tags.values)) # printing the unique indices
[11 12]
Print out the unique facet index numbers
top_imap = domain.topology.index_map(1) # index map of 2D entities in domain
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
# Physical Curve("Top", 9) = {8};
# Physical Curve("Bottom", 10) = {6};
# Physical Curve("MRE_surf", 13) = {1};
# Physical Curve("Right", 14) = {7};
# Physical Curve("Left", 15) = {2, 3, 4, 5};
[ 9 10 13 14 15]
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, opacity=0.25)
# Add colored 1D surfaces for the named surfaces
mre_surf = pyvista.UnstructuredGrid(*vtk_mesh(domain, domain.topology.dim,cell_tags.indices[cell_tags.values==12]) )
air_surf = pyvista.UnstructuredGrid(*vtk_mesh(domain, domain.topology.dim,cell_tags.indices[cell_tags.values==11]) )
#
actor = plotter.add_mesh(mre_surf, show_edges=False,color="red") # mre is red
actor2 = plotter.add_mesh(air_surf, show_edges=False,color="blue") # air is blue
plotter.view_xy()
labels = dict(zlabel='Z', xlabel='X', ylabel='Y')
plotter.add_axes(**labels)
plotter.screenshot("results/sphere_mesh.png")
from IPython.display import Image
Image(filename='results/sphere_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':2})
# Define the volume integration measure "dx"
# also specify the number of volume quadrature points.
dx = ufl.Measure('dx', domain=domain, subdomain_data=cell_tags, 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])
Material parameters#
# A function for constructing spatially varying (piecewise-constant) material parameters
# Need some extra infrastructure for the spatially-discontinuous material property fields
Vmat = functionspace(domain, ("DG", 0)) # create a DG0 function space on the domain
def mat(prop_val_mre, prop_val_air):
# Define an empty "prop" material parameter function,
# which lives on the DG0 function space.
prop = Function(Vmat)
# Now, actualy assign the desired values of shear moduli to the new field.
#
coords = Vmat.tabulate_dof_coordinates()
#
# loop over the coordinates and assign the relevant material property,
# based on the local cell tag number.
for i in range(coords.shape[0]):
if cell_tags.values[i] == 12:
prop.vector.setValueLocal(i, prop_val_mre)
else:
prop.vector.setValueLocal(i, prop_val_air)
return prop
# Elasticity parameters
Gshear = mat(411.0, 0.001) # Shear modulus, kPa
Kbulk = mat(411.0e3, 101.0) # Nearly-incompressible
# Mass density
rho = mat(2000e-6, 1.293e-6) # 1e3 kg/m^3 = 1e-6 kg/mm^3
# Magnetization parameters
#
# Vacuum permeability
mu0 = Constant(domain, 1.256e-6*1e9) # Vacuum permeability, mN / mA^2
#
# material permeability (paramagnetic response)
chi = Constant(domain, 1.75) #mat(2.0, 0.0)# unitless magnetic susceptibility
mu = mu0*(1.0 + chi) # magnetic permeability
ms = Constant(domain, 0.495) # magnetic saturation value, kA/mm = MA/m
# Switch for whether to include the Maxwell stress
TmaxInd = mat(1, 1)
#
# Switch for whether the domain is the air or not
TairInd = mat( 0, 1)
Function spaces#
U2 = element("Lagrange", domain.basix_cell(), 2, shape=(2,)) # For displacement
P1 = element("DG", domain.basix_cell(), 0) # For pressure (must be DG to accomodate different materials)
P2 = element("Lagrange", domain.basix_cell(), 2) # For phi
#
TH = mixed_element([U2, P1, P2]) # Taylor-Hood style mixed element
ME = functionspace(domain, TH) # Total space for all DOFs
#
V1 = functionspace(domain, P1) # Scalar function space.
V1_2 = functionspace(domain, P2) # Scalar (quadratic) function space.
V2 = functionspace(domain, U2) # Vector function space
# V3 = functionspace(domain, T0) # Tensor function space
#
# Define actual functions with the required DOFs
w = Function(ME)
u, p, phi = split(w)
# 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 automatically differentiated.
def safe_sqrt(x):
return sqrt(x + 1.0e-16)
#-------------------------------------------------------
# Subroutines for axisymmetry
#-------------------------------------------------------
# Special gradient operators for axisymmetric functions
#
# Gradient of vector field u
def axi_grad_vector(u):
grad_u = grad(u)
axi_grad_33_exp = conditional(eq(x[0], 0), 0.0, u[0]/x[0])
axi_grad_u = ufl.as_tensor([[grad_u[0,0], grad_u[0,1], 0],
[grad_u[1,0], grad_u[1,1], 0],
[0, 0, axi_grad_33_exp]])
return axi_grad_u
# Gradient of scalar field y
# (just need an extra zero for dimensions to work out)
def axi_grad_scalar(y):
grad_y = grad(y)
axi_grad_y = ufl.as_vector([grad_y[0], grad_y[1], 0.])
return axi_grad_y
# Axisymmetric deformation gradient
def F_axi_calc(u):
dim = len(u) # dimension of problem (2)
Id = Identity(dim) # 2D Identity tensor
F = Id + grad(u) # 2D Deformation gradient
F33_exp = 1.0 + u[0]/x[0] # axisymmetric F33, R/R0
F33 = conditional(eq(x[0], 0.0), 1.0, F33_exp) # avoid divide by zero at r=0
F_axi = ufl.as_tensor([[F[0,0], F[0,1], 0],
[F[1,0], F[1,1], 0],
[0, 0, F33]]) # Full axisymmetric F
return F_axi
def h_calc(u, phi):
F = F_axi_calc(u)
h_R = -axi_grad_scalar(phi)
h_sp = inv(F.T)*h_R
return h_sp
#-------------------------------------------------------------
# Subroutines for calculating the equilibrium Piola stress
#-------------------------------------------------------------
# Subrountine for computing the effective stretch
def lambdaBar_calc(u):
F = F_axi_calc(u)
J = det(F)
Fbar = J**(-1/3)*F
Cbar = Fbar.T*Fbar
I1 = tr(Cbar)
lambdaBar = safe_sqrt(I1/3.0)
return lambdaBar
# Subroutine for calculating the equilibrium Cauchy stress
def T_eq_calc(u,p):
F = F_axi_calc(u)
J = det(F)
Fbar = J**(-1/3)*F
Bbar = Fbar*Fbar.T
T_eq = (1/J)* Gshear * dev(Bbar) - p * Identity(3)
return T_eq
# Subroutine for calculating the equilibrium Cauchy stress
def T_eq_vis_calc(u,p):
F = F_axi_calc(u)
J = det(F)
Fbar = J**(-1/3)*F
Bbar = Fbar*Fbar.T
T_eq = (1/J)* Gshear * dev(Bbar) - p * Identity(3)
return T_eq
# Subroutine for calculating the Magnetic Cauchy stress
def T_mag_calc(u, phi):
h_sp = h_calc(u, phi)
h_norm = safe_sqrt( dot(h_sp, h_sp) )
T_mag_maxw = mu0 * (outer(h_sp, h_sp) - (1./2.) * dot(h_sp, h_sp) * Identity(3) )
T_mag_para = mu * (outer(h_sp, h_sp) - (1./2.) * dot(h_sp, h_sp) * Identity(3) )
# - mu0 * ( (ms**2) / chi ) * ln( cosh( ( chi / ms ) * h_norm ) ) * Identity(3) \
# + mu0 * ( ms / h_norm ) * tanh( ( chi / ms ) * h_norm ) * outer(h_sp, h_sp)
# Exclude the paramagnetic response in the air domain.
# T_mag = T_mag_maxw + conditional(eq(TairInd,1), 0.0*Identity(3), T_mag_para)
T_mag = conditional(eq(TairInd,1), T_mag_maxw, T_mag_para)
return T_mag
def Piola_calc(u, p, phi):
F = F_axi_calc(u)
J = det(F)
T_eq = T_eq_calc(u,p)
T_mag = T_mag_calc(u, phi)
T = T_eq + conditional(eq(TairInd,1), 0*Identity(3), T_mag)
Piola = J * T * inv(F.T)
return Piola
#-------------------------------------------------------------
# Subroutines for calculating the magnetic flux density
#-------------------------------------------------------------
def b_R_calc(u, phi):
F = F_axi_calc(u)
J = det(F)
h_sp = h_calc(u, phi)
h_norm = safe_sqrt( dot(h_sp, h_sp) )
b_maxw = mu0 * h_sp
b_para = mu * h_sp # mu0 * ( ms / h_norm ) * tanh( ( chi / ms ) * h_norm ) * h_sp
# Exclude the paramagnetic response in the air domain.
# b_sp = b_maxw + conditional(eq(TairInd,1), ufl.as_vector([0.,0.,0.]), b_para)
b_sp = conditional(eq(TairInd,1), b_maxw, b_para)
b_R = J * inv(F) * b_sp
return b_R
def Piola_out(u, phi):
F = F_axi_calc(u)
J = det(F)
# Calculate the Maxwell stress in the air
h_sp = h_calc(u, phi)
T_max = mu0 * (outer(h_sp, h_sp) - (1./2.) * dot(h_sp, h_sp) * Identity(3) )
# map it to the Piola stress
Piola_out = J*T_max*inv(F.T)
return Piola_out
def Cauchy_out(u, phi):
# Calculate the Maxwell stress in the air
h_sp = h_calc(u, phi)
T_max = mu0 * (outer(h_sp, h_sp) - (1./2.) * dot(h_sp, h_sp) * Identity(3) )
return T_max
Evaluate kinematics and constitutive relations#
# kinematical quantities
F = F_axi_calc(u)
J = det(F)
lambdaBar = lambdaBar_calc(u)
# Evaulate the total Piola stress
Piola = Piola_calc(u, p, phi)
# Evaluate the magnetic flux density
b_R = b_R_calc(u, phi)
Weak forms#
# Where TairInd on the '+' side is equal to 1, the '-' side indicates the MRE.
# Thus we use the following form of the maxwell traction, where:
# -> J, F, and n are those of the MRE on the '-' side.
# -> Cauchy_out is that of the Air on the '+' side.
#
# The leading factor of TairInd('+') acts as a conditional,
# which makes sure that this traction is zero where the MRE is on the '+' side.
#
maxw_trac_plus = TairInd('+')*dot(J('-')*Cauchy_out(u, phi)('+')*inv(F.T)('-'),n('-'))
# Where TairInd on the '-' side is equal to 1, the '+' side indicates the MRE.
# Thus we use the following form of the maxwell traction, where:
# -> J, F, and n are those of the MRE on the '+' side.
# -> Cauchy_out is that of the Air on the '-' side.
#
# The leading factor of TairInd('-') acts as a conditional,
# which makes sure that this traction is zero where the MRE is on the '-' side.
#
maxw_trac_minus = TairInd('-')*dot(J('+')*Cauchy_out(u, phi)('-')*inv(F.T)('+'),n('+'))
# The total maxwell traction on ds_out will be the contribution of both of these terms, where:
# -> maxw_trac_plus applies on portions of ds_out where the air is on the '+' side
# -> maxw_trac_minus applies on portions of ds_out where the air is on the '-' side
# -> the whole boundary must have one of these two conditions, so together they cover the whole domain.
#
maxw_trac = maxw_trac_plus + maxw_trac_minus
# The weak form for the equilibrium equation
#
Res_1 = inner( Piola, axi_grad_vector(u_test))*x[0]*dx \
- inner(maxw_trac, ufl.avg(ufl.as_vector([u_test[0], u_test[1], 0.0])) )*x[0]*ds(13)
# The auxiliary equation for the pressure
Res_2 = inner((J-1) + p/Kbulk, p_test)*x[0]*dx
# Weak form of Gauss's Law for magnetism
Res_3 = inner( b_R, axi_grad_scalar(phi_test) )*x[0]*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 = "Axi_sphere_air_5mm_ideal"
# Function space for projection of results
P1 = element("DG", domain.basix_cell(), 1)
VV1 = fem.functionspace(domain, P1) # linear scalar function space
#
U1 = element("DG", domain.basix_cell(), 1, shape=(2,))
VV2 = fem.functionspace(domain, U1) # linear Vector function space
#
T1 = element("DG", 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.
#
# Create a linear problem for projecting the stress tensor onto the linear tensor function space VV3.
def setup_projection(u, V):
trial = ufl.TrialFunction(V)
test = ufl.TestFunction(V)
a = ufl.inner(trial, test)*x[0]*dx
L = ufl.inner(u, test)*x[0]*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
Piola_projection = setup_projection(Piola, VV3)
Piola_temp = Piola_projection.solve()
T_mag = T_mag_calc(u, phi)
T_eq = T_eq_vis_calc(u, p)
T = Piola*F.T/J
T22_projection = setup_projection(T[1,1], VV1)
T22_temp = T22_projection.solve()
T22_temp.name = "T22 projection"
T22_projection1 = setup_projection(T_mag[1,1], VV1)
T22_temp1 = T22_projection1.solve()
T22_temp1.name = "T22 mag projection"
p_proj = setup_projection(p, VV1)
p_vis = p_proj.solve()
p_vis.name = "p"
# 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())
T12 = Function(VV1)
T12.name = "T12"
T12_expr = Expression(T[0,1],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())
F22 = Function(VV1)
F22.name = "F22"
F22_expr = Expression(F[1,1],VV1.element.interpolation_points())
# Effective stretch
lambdaBar_vis = Function(VV1)
lambdaBar_vis.name = "LambdaBar"
lambdaBar_expr = Expression(lambdaBar, VV1.element.interpolation_points())
# Volumetric deformation
J_vis = Function(VV1)
J_vis.name = "J"
J_expr = Expression(J, VV1.element.interpolation_points())
# air flag
TairInd_vis = Function(VV1)
TairInd_vis.name = "Air flag"
TairInd_expr = Expression(TairInd, VV1.element.interpolation_points())
# h-, b-, and m-fields
h_vis = Function(VV2)
h_vis.name = "h-field"
h_sp = h_calc(u, phi)
h_2D = ufl.as_vector([h_sp[0], h_sp[1]])
h_vis_expr = Expression(h_2D, VV2.element.interpolation_points())
#
b_vis = Function(VV2)
b_vis.name = "b-field, T"
b_sp = (F * b_R_calc(u, phi) / J)/1e3
b_2D = ufl.as_vector([b_sp[0], b_sp[1]])
b_vis_expr = Expression(b_2D, VV2.element.interpolation_points())
#
m_vis = Function(VV2)
m_vis.name = "m-field"
m_sp = b_sp*1e3/mu0 - h_sp
m_2D = ufl.as_vector([m_sp[0], m_sp[1]])
m_vis_expr = Expression(m_2D, VV2.element.interpolation_points())
# Magnitude of magnetic fields
h_mag_vis = Function(VV1)
h_mag_vis.name = "h-field magnitude"
h_mag = safe_sqrt(dot(h_sp, h_sp))
h_mag_expr = Expression(h_mag, VV1.element.interpolation_points())
#
b_mag_vis = Function(VV1)
b_mag_vis.name = "b-field magnitude, T"
b_mag = safe_sqrt(dot(b_sp, b_sp))
b_mag_expr = Expression(b_mag, VV1.element.interpolation_points())
#
m_mag_vis = Function(VV1)
m_mag_vis.name = "m-field magnitude"
m_mag = safe_sqrt(dot(m_sp, m_sp))
m_mag_expr = Expression(m_mag, VV1.element.interpolation_points())
#
m_mag_proj = setup_projection(safe_sqrt(dot(m_sp, m_sp)), VV1)
m_mag_proj_vis = m_mag_proj.solve()
m_mag_proj_vis.name = "m proj."
# induced h- and b-fields
#
# Constant for applied h-field
h_cons = Constant(domain,PETSc.ScalarType(hRamp(0)))
#
h_ind_vis = Function(VV2)
h_ind_vis.name = "Induced h-field"
h_ind_sp = h_calc(u, phi) - ufl.as_vector([0, h_cons, 0])
h_ind_2D = ufl.as_vector([h_ind_sp[0], h_ind_sp[1]])
h_ind_vis_expr = Expression(h_ind_2D, VV2.element.interpolation_points())
#
b_ind_vis = Function(VV2)
b_ind_vis.name = "Induced b-field, T"
b_ind_sp = (F * b_R_calc(u, phi) / J - mu0 * ufl.as_vector([0, h_cons, 0]))/1e3
b_ind_2D = ufl.as_vector([b_ind_sp[0], b_ind_sp[1]])
b_ind_vis_expr = Expression(b_ind_2D, VV2.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, T12, T22, T22_temp, T22_temp1, T33, F22, # stress outputs
lambdaBar_vis, J_vis, # Kinematical outputs
TairInd_vis, # material params
h_vis, b_vis, m_vis, m_mag_proj_vis, h_ind_vis, b_ind_vis, # magnetic 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))
p_vis = p_proj.solve()
phi_vis.interpolate(w.sub(2))
#
# re-project to smooth visualization of quadrature functions
# before interpolating.
Piola_temp = Piola_projection.solve()
T22_temp = T22_projection.solve()
T22_temp1 = T22_projection1.solve()
Mises_vis.interpolate(Mises_expr)
T11.interpolate(T11_expr)
T12.interpolate(T12_expr)
T22.interpolate(T22_expr)
T33.interpolate(T33_expr)
F22.interpolate(F22_expr)
#
lambdaBar_vis.interpolate(lambdaBar_expr)
J_vis.interpolate(J_expr)
#
TairInd_vis.interpolate(TairInd_expr)
#
h_vis.interpolate(h_vis_expr)
b_vis.interpolate(b_vis_expr)
m_vis.interpolate(m_vis_expr)
#
h_mag_vis.interpolate(h_mag_expr)
b_mag_vis.interpolate(b_mag_expr)
m_mag_vis.interpolate(m_mag_expr)
#
m_mag_proj_vis = m_mag_proj.solve()
#
h_ind_vis.interpolate(h_ind_vis_expr)
b_ind_vis.interpolate(b_ind_vis_expr)
# Finally, write output fields to VTX.
#
file_results.write(t)
P1_Lg = element("Lagrange", domain.basix_cell(), 1)
VV1_Lg = fem.functionspace(domain, P1_Lg) # linear lagrange scalar function space
TairInd_lg_vis = Function(VV1_Lg, name="Air flag")
U1_Lg = element("Lagrange", domain.basix_cell(), 1, shape=(2,))
VV2_Lg = fem.functionspace(domain, U1_Lg) # linear lagrange Vector function space
u_lg_vis = Function(VV2_Lg, name="disp")
# set up the output VTX files.
file_results_bound = VTXWriter(
MPI.COMM_WORLD,
"results/" + results_name + "_boundary.bp",
[ # put the functions here you wish to write to output
u_lg_vis, TairInd_lg_vis,
],
engine="BP4",
)
def writeResults_bound(t):
# Update the output fields before writing to VTX.
#
u_lg_vis.interpolate(w.sub(0))
TairInd_lg_vis.interpolate(TairInd_expr)
# Finally, write output fields to VTX.
#
file_results_bound.write(t)
Infrastructure for pulling out time history data (force, displacement, etc.)#
# infrastructure for evaluating functions at a certain point efficiently
pointForEval = np.array([0.0, radius, 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
#
# infrastructure for evaluating functions at a certain point efficiently
pointForEval2 = np.array([0.0, 0.0, 0.0])
cell_candidates2 = dolfinx.geometry.compute_collisions_points(bb_tree, pointForEval2)
colliding_cells2 = dolfinx.geometry.compute_colliding_cells(domain, cell_candidates2, pointForEval2).array
#
pointForEval3 = np.array([0.0, -radius, 0.0])
cell_candidates3 = dolfinx.geometry.compute_collisions_points(bb_tree, pointForEval3)
colliding_cells3 = dolfinx.geometry.compute_colliding_cells(domain, cell_candidates3, pointForEval3).array
Boundary condtions#
# Physical Curve("Top", 9) = {8};
# Physical Curve("Bottom", 10) = {6};
# Physical Curve("MRE_surf", 13) = {1};
# Physical Curve("Right", 14) = {7};
# Physical Curve("Left", 15) = {2, 3, 4, 5};
# Constant for applied displacement
phi_cons = Constant(domain,PETSc.ScalarType(phiRamp(0)))
# 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(15))
#
right_u1_dofs = fem.locate_dofs_topological(ME.sub(0).sub(0), facet_tags.dim, facet_tags.find(14))
top_u1_dofs = fem.locate_dofs_topological(ME.sub(0).sub(0), facet_tags.dim, facet_tags.find(9))
btm_u1_dofs = fem.locate_dofs_topological(ME.sub(0).sub(0), facet_tags.dim, facet_tags.find(10))
#
right_u2_dofs = fem.locate_dofs_topological(ME.sub(0).sub(1), facet_tags.dim, facet_tags.find(14))
top_u2_dofs = fem.locate_dofs_topological(ME.sub(0).sub(1), facet_tags.dim, facet_tags.find(9))
btm_u2_dofs = fem.locate_dofs_topological(ME.sub(0).sub(1), facet_tags.dim, facet_tags.find(10))
#
top_phi_dofs = fem.locate_dofs_topological(ME.sub(2), facet_tags.dim, facet_tags.find(9))
btm_phi_dofs = fem.locate_dofs_topological(ME.sub(2), facet_tags.dim, facet_tags.find(10))
# Zero bc quantity
V0, submap = ME.sub(0).sub(1).collapse()
fixed_u_dofs = fem.locate_dofs_geometrical((ME.sub(0).sub(1), V0), ground)
fixed_disp = Function(V0)
fixed_disp.interpolate(lambda x: np.zeros(x.shape[1]) )
# building Dirichlet BCs
bcs_0 = dirichletbc(0.0, left_u1_dofs, ME.sub(0).sub(0)) # u1 fix - left
#
bcs_1 = dirichletbc(0.0, right_u1_dofs, ME.sub(0).sub(0)) # u1 fix - right
bcs_2 = dirichletbc(0.0, top_u1_dofs, ME.sub(0).sub(0)) # u1 fix - top
bcs_3 = dirichletbc(0.0, btm_u1_dofs, ME.sub(0).sub(0)) # u1 fix - btm
#
bcs_4 = dirichletbc(0.0, right_u2_dofs, ME.sub(0).sub(1)) # u2 fix - right
bcs_5 = dirichletbc(0.0, top_u2_dofs, ME.sub(0).sub(1)) # u2 fix - top
bcs_6 = dirichletbc(0.0, btm_u2_dofs, ME.sub(0).sub(1)) # u2 fix - btm
#
bcs_7 = dirichletbc(0.0, top_phi_dofs, ME.sub(2)) # phi fix - Top
bcs_8 = dirichletbc(phi_cons, btm_phi_dofs, ME.sub(2)) # phi ramp - Btm
#
bcs_9 = dirichletbc(fixed_disp, fixed_u_dofs, ME.sub(0).sub(1)) # u fix - center node
bcs = [bcs_0, bcs_1, bcs_2, bcs_3, bcs_4, bcs_5, bcs_6, bcs_7, bcs_8, bcs_9]
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])
timeHist2a = np.zeros(shape=[totSteps])
timeHist2b = np.zeros(shape=[totSteps])
timeHist3 = np.zeros(shape=[totSteps])
timeHist4 = np.zeros(shape=[totSteps])
timeHist5 = np.zeros(shape=[totSteps])
#Iinitialize a counter for reporting data
ii=0
# Write initial state to file
writeResults(t=0.0)
writeResults_bound(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 = float(phiRamp(t))
h_cons.value = float(hRamp(t))
# Solve the problem
(iter, converged) = solver.solve(w)
# Collect results from MPI ghost processes
w.x.scatter_forward()
# Print progress of calculation periodically
if ii%1 == 0:
now = datetime.now()
current_time = now.strftime("%H:%M:%S")
print("Step: {} | Increment: {} | Iterations: {}".format(step, ii, iter))
print("dt: {} s | Simulation Time: {} s of {} s".format(round(dt, 4), round(t,4), Ttot))
# print("Maximum viscous strain increment: {} ".format(round(gbarV_delta, 6)))
print()
# Write output to file
writeResults(t)
writeResults_bound(t)
# Store time history variables at this time
timeHist0[ii] = t # current time
#
timeHist1[ii] = hRamp(t) # time history of applied h-field
#
timeHist2a[ii] = w.sub(0).sub(1).eval([0, radius, 0],colliding_cells[0])[0] # time history of vertical displacement
#
timeHist2b[ii] = w.sub(0).sub(1).eval([0, -radius, 0],colliding_cells3[0])[0] # time history of vertical displacement
#
timeHist3[ii] = h_mag_vis.eval([0, 0, 0],colliding_cells2[0])[0] # time history of h-field
#
timeHist4[ii] = b_mag_vis.eval([0, 0, 0],colliding_cells2[0])[0] # time history of b-field
#
timeHist5[ii] = m_mag_vis.eval([0, 0, 0],colliding_cells2[0])[0] # time history of m-field
# Update DOFs for next step
w_old.x.array[:] = w.x.array
# close the output file.
file_results.close()
file_results_bound.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: 25.0 s | Simulation Time: 25.0 s of 1000.0 s
Step: Actuate | Increment: 2 | Iterations: 4
dt: 25.0 s | Simulation Time: 50.0 s of 1000.0 s
Step: Actuate | Increment: 3 | Iterations: 4
dt: 25.0 s | Simulation Time: 75.0 s of 1000.0 s
Step: Actuate | Increment: 4 | Iterations: 4
dt: 25.0 s | Simulation Time: 100.0 s of 1000.0 s
Step: Actuate | Increment: 5 | Iterations: 4
dt: 25.0 s | Simulation Time: 125.0 s of 1000.0 s
Step: Actuate | Increment: 6 | Iterations: 4
dt: 25.0 s | Simulation Time: 150.0 s of 1000.0 s
Step: Actuate | Increment: 7 | Iterations: 4
dt: 25.0 s | Simulation Time: 175.0 s of 1000.0 s
Step: Actuate | Increment: 8 | Iterations: 4
dt: 25.0 s | Simulation Time: 200.0 s of 1000.0 s
Step: Actuate | Increment: 9 | Iterations: 4
dt: 25.0 s | Simulation Time: 225.0 s of 1000.0 s
Step: Actuate | Increment: 10 | Iterations: 4
dt: 25.0 s | Simulation Time: 250.0 s of 1000.0 s
Step: Actuate | Increment: 11 | Iterations: 4
dt: 25.0 s | Simulation Time: 275.0 s of 1000.0 s
Step: Actuate | Increment: 12 | Iterations: 4
dt: 25.0 s | Simulation Time: 300.0 s of 1000.0 s
Step: Actuate | Increment: 13 | Iterations: 4
dt: 25.0 s | Simulation Time: 325.0 s of 1000.0 s
Step: Actuate | Increment: 14 | Iterations: 4
dt: 25.0 s | Simulation Time: 350.0 s of 1000.0 s
Step: Actuate | Increment: 15 | Iterations: 4
dt: 25.0 s | Simulation Time: 375.0 s of 1000.0 s
Step: Actuate | Increment: 16 | Iterations: 4
dt: 25.0 s | Simulation Time: 400.0 s of 1000.0 s
Step: Actuate | Increment: 17 | Iterations: 4
dt: 25.0 s | Simulation Time: 425.0 s of 1000.0 s
Step: Actuate | Increment: 18 | Iterations: 4
dt: 25.0 s | Simulation Time: 450.0 s of 1000.0 s
Step: Actuate | Increment: 19 | Iterations: 4
dt: 25.0 s | Simulation Time: 475.0 s of 1000.0 s
Step: Actuate | Increment: 20 | Iterations: 4
dt: 25.0 s | Simulation Time: 500.0 s of 1000.0 s
Step: Actuate | Increment: 21 | Iterations: 4
dt: 25.0 s | Simulation Time: 525.0 s of 1000.0 s
Step: Actuate | Increment: 22 | Iterations: 4
dt: 25.0 s | Simulation Time: 550.0 s of 1000.0 s
Step: Actuate | Increment: 23 | Iterations: 4
dt: 25.0 s | Simulation Time: 575.0 s of 1000.0 s
Step: Actuate | Increment: 24 | Iterations: 4
dt: 25.0 s | Simulation Time: 600.0 s of 1000.0 s
Step: Actuate | Increment: 25 | Iterations: 4
dt: 25.0 s | Simulation Time: 625.0 s of 1000.0 s
Step: Actuate | Increment: 26 | Iterations: 4
dt: 25.0 s | Simulation Time: 650.0 s of 1000.0 s
Step: Actuate | Increment: 27 | Iterations: 4
dt: 25.0 s | Simulation Time: 675.0 s of 1000.0 s
Step: Actuate | Increment: 28 | Iterations: 4
dt: 25.0 s | Simulation Time: 700.0 s of 1000.0 s
Step: Actuate | Increment: 29 | Iterations: 4
dt: 25.0 s | Simulation Time: 725.0 s of 1000.0 s
Step: Actuate | Increment: 30 | Iterations: 4
dt: 25.0 s | Simulation Time: 750.0 s of 1000.0 s
Step: Actuate | Increment: 31 | Iterations: 4
dt: 25.0 s | Simulation Time: 775.0 s of 1000.0 s
Step: Actuate | Increment: 32 | Iterations: 4
dt: 25.0 s | Simulation Time: 800.0 s of 1000.0 s
Step: Actuate | Increment: 33 | Iterations: 4
dt: 25.0 s | Simulation Time: 825.0 s of 1000.0 s
Step: Actuate | Increment: 34 | Iterations: 4
dt: 25.0 s | Simulation Time: 850.0 s of 1000.0 s
Step: Actuate | Increment: 35 | Iterations: 4
dt: 25.0 s | Simulation Time: 875.0 s of 1000.0 s
Step: Actuate | Increment: 36 | Iterations: 4
dt: 25.0 s | Simulation Time: 900.0 s of 1000.0 s
Step: Actuate | Increment: 37 | Iterations: 4
dt: 25.0 s | Simulation Time: 925.0 s of 1000.0 s
Step: Actuate | Increment: 38 | Iterations: 4
dt: 25.0 s | Simulation Time: 950.0 s of 1000.0 s
Step: Actuate | Increment: 39 | Iterations: 4
dt: 25.0 s | Simulation Time: 975.0 s of 1000.0 s
Step: Actuate | Increment: 40 | Iterations: 4
dt: 25.0 s | Simulation Time: 1000.0 s of 1000.0 s
-----------------------------------------
End computation
------------------------------------------
Elapsed real time: 0:00:12.543437
------------------------------------------
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) +1
plt.figure()
plt.axvline(0, c='k', linewidth=1.)
plt.axhline(0, c='k', linewidth=1.)
plt.plot(float(mu0)*timeHist1[0:ind]/1e3, (timeHist2a[0:ind] - timeHist2b[0:ind])/(2*radius), linewidth=2.0, color=colors[2])
plt.axis('tight')
plt.ylabel(r"$\lambda$ - 1")
plt.xlabel(r"$\mu_0 \, |\mathbf{h}^{app}|$ (Teslas)")
plt.grid(linestyle="--", linewidth=0.5, color='b')
plt.ylim(0,0.25)
plt.xlim(0,1.0)
# plt.legend()
fig = plt.gcf()
fig.set_size_inches(7,5)
plt.tight_layout()
plt.savefig("results/2D_ax_sphere_air_stretch.png", dpi=600)
stretchData = np.array([float(mu0)*timeHist1[0:ind]/1e3, (timeHist2a[0:ind] - timeHist2b[0:ind])/(2*radius)])
np.savetxt("results/5mm_stretch_data_linear.csv", stretchData, delimiter=",")
plt.figure()
plt.axvline(0, c='k', linewidth=1.)
plt.axhline(0, c='k', linewidth=1.)
plt.plot(float(mu0)*timeHist1[0:ind]/1e3, timeHist4[0:ind], linewidth=2.0,\
color=colors[0], label=r'$|\mathbf{b}|$')
plt.plot(float(mu0)*timeHist1[0:ind]/1e3, timeHist3[0:ind]*float(mu0)/1e3, linewidth=2.0,\
color=colors[1], label=r'$\mu_0\, |\mathbf{h}|$')
plt.plot(float(mu0)*timeHist1[0:ind]/1e3, timeHist5[0:ind]*float(mu0)/1e3, linewidth=2.0,\
color=colors[2], label=r'$\mu_0\, |\mathbf{m}|$')
plt.axis('tight')
# plt.ylabel(r"$|\mathbf{b}| , \ \mu_0\,|\mathbf{h}|, \ \mu_0\,|\mathbf{m}|$, Teslas")
plt.ylabel(r"Internal field strength (Teslas)")
plt.xlabel(r"$\mu_0 \, |\mathbf{h}^{app}|$ (Teslas)")
plt.grid(linestyle="--", linewidth=0.5, color='b')
plt.ylim(0,2.0)
plt.xlim(0,1.0)
plt.legend()
fig = plt.gcf()
fig.set_size_inches(7,5)
plt.tight_layout()
plt.savefig("results/2D_ax_sphere_air_fields.png", dpi=600)
